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;
        }
}