GIS矢量切片算法（转载）

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