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 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();

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;

}

}