GIS矢量切片算法
作者：互联网
转自： https://www.giserdqy.com/database/postgresql/25838/
参考：
https://my.oschina.net/u/1464512/blog/1631972
https://github.com/mapbox/tippecanoe
https://github.com/mapbox/tile-cover
https://github.com/mapbox/vector-tile-base
https://github.com/mapbox/vtzero
https://github.com/mapbox/mapnik-vector-tile
对于大范围矢量数据，由于类型众多范围广泛往往数据量极大，加载渲染会造成平台卡顿。因此对矢量数据进行四叉树索引切片可以高效加载当前区域矢量，提高效率。
常见的矢量数据为shapefile，可以通过GDAL读取shp范围进行四叉树划分，构建某一层级瓦块。
以下为C#调用GDAL进行矢量四叉树切片算法：
struct TileStructure
{
public int level;
public int x;
public int y;
public OSGeo.OGR.Geometry extentPolygon;
public string path;
public OSGeo.OGR.DataSource ds;
public OSGeo.OGR.Layer layer;
}
public class VectorTileTool
{
List<TileStructure> tiles;
public VectorTileTool()
{
}
public bool SeprateShpLayer(string sourcePath, string resultFolder, int level)
{
OSGeo.GDAL.Gdal.SetConfigOption("SHAPE_ENCODING", "");
OSGeo.OGR.Ogr.RegisterAll();
OSGeo.OGR.Driver dr = OSGeo.OGR.Ogr.GetDriverByName("ESRI shapefile");
if (dr == null)
{
return false;
}
OSGeo.OGR.DataSource ds = dr.Open(sourcePath, 0);
int layerCount = ds.GetLayerCount();
OSGeo.OGR.Layer layer = ds.GetLayerByIndex(0);
//投影信息
OSGeo.OSR.SpatialReference coord = layer.GetSpatialRef();
string coordString;
coord.ExportToWkt(out coordString);
//地理范围
Envelope layerEx = new Envelope();
layer.GetExtent(layerEx, 0);
////如果瓦块数据存在，全部删除
//if (Directory.Exists(resultFolder))
//{
// Directory.Delete(resultFolder,true);
//}
//创建文件夹
Directory.CreateDirectory(resultFolder + "\\JSON\\");
//针对本项目，划分第16级，根据地理范围求出瓦片
int y0 = Convert.ToInt32((90 - layerEx.MaxY) * Math.Pow(2, level)/180.0);
int x0 = Convert.ToInt32((180 + layerEx.MinX) * Math.Pow(2, level)/180.0);
int y1 = Convert.ToInt32((90 - layerEx.MinY) * Math.Pow(2, level) / 180.0);
int x1 = Convert.ToInt32((180 + layerEx.MaxX) * Math.Pow(2, level) / 180.0);
//20160621 ZXQ 创建层行列配置文件
string filePath = resultFolder + "\\JSON\\" + "\\tile.txt";
FileStream fs = new FileStream(filePath, FileMode.CreateNew);
StreamWriter sw = new StreamWriter(fs);
//写入层行列
sw.Write(level.ToString());
sw.Write(",");
sw.Write(x0.ToString());
sw.Write(",");
sw.Write(x1.ToString());
sw.Write(",");
sw.Write(y0.ToString());
sw.Write(",");
sw.Write(y1.ToString());
sw.Write(",");
sw.Write("json");
sw.Flush();
sw.Close();
fs.Close();
tiles = new List<TileStructure>();
for (int x =x0;x<=x1;x++)
{
for (int y=y0;y<=y1;y++)
{
TileStructure tile;
tile.level = level;
tile.x = x;
tile.y = y;
double lonMin = -180 + 180 / (Math.Pow(2, level)) * x;
double lonMax = -180 + 180 / (Math.Pow(2, level)) * (x + 1);
double latMax = 90 - 180 / (Math.Pow(2, level)) * y;
double latMin = 90 - 180 / (Math.Pow(2, level)) * (y + 1);
tile.extentPolygon = new OSGeo.OGR.Geometry(OSGeo.OGR.wkbGeometryType.wkbPolygon);
OSGeo.OGR.Geometry geo = new OSGeo.OGR.Geometry(OSGeo.OGR.wkbGeometryType.wkbLinearRing);
geo.AddPoint(lonMin,latMax,0);
geo.AddPoint(lonMax, latMax, 0);
geo.AddPoint(lonMin, latMin, 0);
geo.AddPoint(lonMax, latMin, 0);
tile.extentPolygon.AddGeometryDirectly(geo);
tile.extentPolygon.CloseRings();
//创建空shp文件
string tileFolder = resultFolder + "\\SHP\\" + level.ToString() + "\\" + x.ToString();
string fileName = y.ToString() + ".shp";
string tilePath = tileFolder + "\\" + fileName;
if (!Directory.Exists(tileFolder))
{
Directory.CreateDirectory(tileFolder);
}
tile.path = tilePath;
tile.ds = dr.CreateDataSource(tilePath, null);
tile.layer = tile.ds.CreateLayer("house", coord, OSGeo.OGR.wkbGeometryType.wkbPolygon, null);
FieldDefn fd = new FieldDefn("HEIGHT", FieldType.OFTReal);
tile.layer.CreateField(fd,1);
tiles.Add(tile);
Console.WriteLine("创建第{0}层第{1}行第{2}列瓦块空shapefile数据", level, x, y);
}
}
OSGeo.OGR.Feature feat;
//读取shp文件
while ((feat = layer.GetNextFeature()) != null)
{
int id = feat.GetFID();
OSGeo.OGR.Geometry geometry = feat.GetGeometryRef();
OSGeo.OGR.wkbGeometryType goetype = geometry.GetGeometryType();
if (goetype != wkbGeometryType.wkbPolygon)
{
continue;
}
geometry.CloseRings();
//随机楼层3-15层
Random random = new Random();
double height = random.Next(3,15)*3;// feat.GetFieldAsDouble("房屋层数") * 3;
for (int i = 0; i < tiles.Count;i++ )
{
TileStructure tile = tiles[i];
//如果瓦片与要素相交，则将要素放入该瓦片
if (tile.extentPolygon.Intersect(geometry))
{
OSGeo.OGR.Feature poFeature = new Feature(tile.layer.GetLayerDefn());
poFeature.SetField(0, height.ToString());
poFeature.SetGeometry(geometry);
tile.layer.CreateFeature(poFeature);
Console.WriteLine("写入第{0}个要素入shp", id);
}
}
}
return true;
}
}
标签：Write,GIS,int,矢量,sw,OGR,切片,OSGeo,public 来源： https://www.cnblogs.com/gispathfinder/p/12929319.html