GIS vector slicing algorithm (Reprint)

Time:2020-11-25

Transferred from: https://www.giserdqy.com/database/postgresql/25838/

For large-scale vector data, due to the large number of types and wide range of data, loading and rendering will cause the platform to be stuck. Therefore, the quadtree index slicing of vector data can efficiently load the current region vector and improve the efficiency.

image

The common vector data is ShapeFile, which can be read from the SHP range by GDAL to divide the quadtree and construct a certain level of tiles.

The following is a vector quadtree slicing algorithm called GDAL by C ා:

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);
        //Projection information
        OSGeo.OSR.SpatialReference coord = layer.GetSpatialRef();
        string coordString;
        coord.ExportToWkt(out coordString);
        //Geographical scope
        Envelope layerEx = new Envelope();
        layer.GetExtent(layerEx, 0);
        ////If tile data exists, delete all tiles
        //if (Directory.Exists(resultFolder))
        //{
        //    Directory.Delete(resultFolder,true);
        //}
        //Create folder
        Directory.CreateDirectory(resultFolder + "\JSON\");
        //For this project, the 16th level is divided and tiles are calculated according to the geographical scope
        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 create layer row and column configuration file
        string filePath = resultFolder + "\JSON\" + "\tile.txt";
        FileStream fs = new FileStream(filePath, FileMode.CreateNew);
        StreamWriter sw = new StreamWriter(fs);
        //Write layer row and column
        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();
                //Create an empty SHP file
                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 (create empty ShapeFile data of tile in {0} layer, line {1} and column {2} ", level, x, y);
            }
        }
        OSGeo.OGR.Feature feat;
        //Read SHP file
        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();
            //Random floor 3-15
            Random random = new Random();
            double height =  random.Next (3,15)*3;//  feat.GetFieldAsDouble (the "number of floors") * 3;
            for (int i = 0; i < tiles.Count;i++ )
            {
                TileStructure tile = tiles[i];
                //If the tile intersects the feature, the feature is placed in the tile
                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 (write the {0} element into SHP, ID);
                }
            }
        }
        return true;
    }
}

———————-Over——————————-