using CalcVolServe.Models;
|
using OSGeo.GDAL;
|
using OSGeo.OGR;
|
using OSGeo.OSR;
|
using System;
|
using System.Collections.Generic;
|
using System.IO;
|
using System.Linq;
|
using System.Reflection;
|
using System.Web;
|
|
namespace CalcVolServe.cs
|
{
|
public class GdalOp
|
{
|
private static bool isInited = false;
|
|
private static object obj = new object();
|
|
private static readonly char[] c1 = new char[] { ',' };
|
|
private static readonly char[] c2 = new char[] { ' ' };
|
|
public GdalOp()
|
{
|
lock (obj)
|
{
|
if (!isInited)
|
{
|
RegisterGDal();
|
isInited = true;
|
}
|
}
|
}
|
|
public void RegisterGDal()
|
{
|
string path = Path.Combine(System.AppDomain.CurrentDomain.BaseDirectory, "GDAL");
|
|
string dir = Environment.GetEnvironmentVariable("PATH");
|
//if (dir.IndexOf(path) == -1)
|
//{
|
// string gPath = dir + path + ";" + Path.Combine(path, "plugins") + ";";
|
// Environment.SetEnvironmentVariable("PATH", gPath);
|
//}
|
|
//string gdalData = Path.Combine(path, "gdal-data");
|
//Environment.SetEnvironmentVariable("GDAL_DATA", gdalData);
|
//Gdal.SetConfigOption("GDAL_DATA", gdalData);
|
|
//Environment.SetEnvironmentVariable("GEOTIFF_CSV", gdalData);
|
//Gdal.SetConfigOption("GEOTIFF_CSV", gdalData);
|
|
//string driverPath = Path.Combine(path, "plugins");
|
//Environment.SetEnvironmentVariable("GDAL_DRIVER_PATH", driverPath);
|
//Gdal.SetConfigOption("GDAL_DRIVER_PATH", driverPath);
|
|
//string projSharePath = Path.Combine(path, "share");
|
//Environment.SetEnvironmentVariable("PROJ_LIB", projSharePath);
|
//Gdal.SetConfigOption("PROJ_LIB", projSharePath);
|
|
Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // NO
|
Gdal.SetConfigOption("SHAPE_ENCODING", "CP936"); // 空,gb2312,CP936
|
|
Ogr.RegisterAll();
|
Gdal.AllRegister();
|
}
|
|
#region 读取shp
|
public List<Triangle> ReadShp(string file, string wkt)
|
{
|
OSGeo.OGR.Driver driver = Ogr.GetDriverByName("ESRI Shapefile");
|
DataSource shpDS = driver.Open(file, 0);
|
if (shpDS == null)
|
{
|
driver.Dispose();
|
return null;
|
}
|
|
Layer fLayer = shpDS.GetLayerByIndex(0);
|
SpatialReference sr = fLayer.GetSpatialRef();
|
Geometry ring = CreateRing(wkt, sr);
|
string layerName = fLayer.GetName();
|
fLayer.Dispose();
|
|
//string name = Path.GetFileNameWithoutExtension(shpFile);
|
string sql = string.Format("select FID from {0}", layerName);
|
|
Layer layer = shpDS.ExecuteSQL(sql, ring, null);
|
if (layer == null)
|
{
|
Dispose(null, driver, shpDS, null);
|
return null;
|
}
|
|
List<Triangle> list = new List<Triangle>();
|
Feature f = layer.GetNextFeature();
|
while (f != null)
|
{
|
Geometry g = f.GetGeometryRef();
|
if (g != null && g.GetGeometryName() == "POLYGON" && g.GetGeometryCount() > 0)
|
{
|
Triangle t = new Triangle();
|
Geometry r = g.GetGeometryRef(0);
|
for (int i = 0, c = r.GetPointCount() - 1; i < c; i++)
|
{
|
double[] arr = new double[] { 0, 0, 0 };
|
r.GetPoint(i, arr);
|
|
Point p = new Point(arr[0], arr[1], arr[2]);
|
t.Points.Add(p);
|
}
|
list.Add(t);
|
}
|
|
f = layer.GetNextFeature();
|
}
|
|
Dispose(null, driver, shpDS, layer);
|
|
return list;
|
}
|
|
public Geometry CreateRing(string wkt, SpatialReference sr)
|
{
|
string str = wkt.Replace("LINEARRING Z (", "").Replace(")", "");
|
string[] strs = str.Split(c1);
|
|
Geometry ring = new Geometry(wkbGeometryType.wkbLinearRing);
|
foreach (string s in strs)
|
{
|
string[] arr = s.Trim().Split(c2);
|
ring.AddPoint_2D(double.Parse(arr[0]), double.Parse(arr[1]));
|
}
|
|
Geometry polygon = new Geometry(wkbGeometryType.wkbPolygon);
|
polygon.AssignSpatialReference(sr);
|
polygon.AddGeometry(ring);
|
|
return polygon;
|
}
|
|
private void Dispose(Dataset ds, OSGeo.OGR.Driver driver, DataSource dS, Layer layer)
|
{
|
if (layer != null)
|
{
|
layer.Dispose();
|
}
|
if (dS != null)
|
{
|
dS.Dispose();
|
}
|
if (driver != null)
|
{
|
driver.Dispose();
|
}
|
if (ds != null)
|
{
|
ds.Dispose();
|
}
|
}
|
#endregion
|
|
#region 读取tif
|
public void ReadTif(string file, List<Triangle> list, string prop)
|
{
|
Type type = typeof(Point);
|
PropertyInfo pi = type.GetProperty(prop);
|
if (pi == null) return;
|
|
Dataset ds = Gdal.Open(file, Access.GA_ReadOnly);
|
if (ds == null) return;
|
|
Band band = ds.GetRasterBand(1);
|
foreach (Triangle t in list)
|
{
|
foreach (Point p in t.Points)
|
{
|
int[] xy = Geo2ImageSpace(ds, p.X, p.Y);
|
if (xy[0] < 1 || xy[1] < 1 || xy[0] > ds.RasterXSize || xy[1] > ds.RasterYSize)
|
{
|
continue;
|
}
|
|
double val = ReadBand(band, xy[0], xy[1]);
|
pi.SetValue(p, val);
|
}
|
}
|
|
band.Dispose();
|
ds.Dispose();
|
}
|
|
private double ReadBand(Band band, int row, int col)
|
{
|
double val = 0;
|
switch (band.DataType)
|
{
|
case DataType.GDT_Byte:
|
//public CPLErr ReadRaster(int xOff, int yOff, int xSize, int ySize, byte[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace);
|
//byte[] byteBuffer = new byte[band.XSize * band.YSize];
|
|
short[] byteBuffer = new short[1];
|
band.ReadRaster(row, col, 1, 1, byteBuffer, 1, 1, 1, 1);
|
val = byteBuffer[0];
|
break;
|
case DataType.GDT_Int16:
|
case DataType.GDT_CInt16:
|
case DataType.GDT_UInt16:
|
//short[] shortBuffer = new short[band.XSize * band.YSize];
|
//public CPLErr ReadRaster(int xOff, int yOff, int xSize, int ySize, short[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace);
|
|
//short[] shortBuffer = new short[band.XSize * band.YSize];
|
//band.ReadRaster(0, 0, band.XSize, band.YSize, shortBuffer, band.XSize, band.YSize, row, col);
|
|
short[] shortBuffer = new short[1];
|
band.ReadRaster(row, col, 1, 1, shortBuffer, 1, 1, 1, 1);
|
val = shortBuffer[0];
|
break;
|
case DataType.GDT_Int32:
|
case DataType.GDT_CInt32:
|
case DataType.GDT_UInt32:
|
//public CPLErr ReadRaster(int xOff, int yOff, int xSize, int ySize, int[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace);
|
|
int[] intBuffer = new int[1];
|
band.ReadRaster(row, col, 1, 1, intBuffer, 1, 1, 1, 1);
|
val = intBuffer[0];
|
break;
|
case DataType.GDT_Float32:
|
case DataType.GDT_CFloat32:
|
//public CPLErr ReadRaster(int xOff, int yOff, int xSize, int ySize, float[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace);
|
|
float[] floatBuffer = new float[1];
|
band.ReadRaster(row, col, 1, 1, floatBuffer, 1, 1, 1, 1);
|
val = floatBuffer[0];
|
break;
|
case DataType.GDT_Float64:
|
case DataType.GDT_CFloat64:
|
//public CPLErr ReadRaster(int xOff, int yOff, int xSize, int ySize, double[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace);
|
|
double[] doubleBuffer = new double[1];
|
band.ReadRaster(row, col, 1, 1, doubleBuffer, 1, 1, 1, 1);
|
val = doubleBuffer[0];
|
break;
|
}
|
|
return val;
|
}
|
|
private int[] Geo2ImageSpace(Dataset ds, double x, double y)
|
{
|
double[] transform = new double[6];
|
ds.GetGeoTransform(transform); // 影像坐标变换参数
|
//string prj = ds.GetProjection(); // 影像坐标系信息（WKT格式字符串）
|
|
int col = (int)((y * transform[1] - x * transform[4] + transform[0] * transform[4] - transform[3] * transform[1]) / (transform[5] * transform[1] - transform[2] * transform[4])); // 像素所在列
|
int row = (int)((x - transform[0] - col * transform[2]) / transform[1]); // 像素所在行
|
|
return new int[] { row, col };
|
}
|
#endregion
|
}
|
}
|
