GIS矢量切片算法（轉載）

轉自： https://www.giserdqy.com/database/postgresql/25838/

對於大範圍矢量數據，因爲類型衆多範圍普遍每每數據量極大，加載渲染會形成平臺卡頓。所以對矢量數據進行四叉樹索引切片能夠高效加載當前區域矢量，提升效率。

常見的矢量數據爲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;
    }
}

----------------------完-------------------------------