using DataLoader.Model; using OSGeo.GDAL; using OSGeo.OGR; using OSGeo.OSR; using System; using System.Collections.Generic; using System.IO; using System.Linq; namespace DataLoader.CS { public class GdalHelper { #region 成员变量+构造函数+初始化 private static bool isInited = false; private static GdalHelper instance = null; private static readonly object obj = new object(); public static SpatialReference sr4326 = null; public static SpatialReference sr4490 = null; public static SpatialReference sr104903 = null; public static string MOON200 = "GCS_Moon_2000"; public static List EPSGS = new List() { "EPSG:4326", "EPSG:4490", "ESRI:104903" }; /// /// 构造函数 /// private GdalHelper() { lock (obj) { if (!isInited) { RegisterGDal(); sr4326 = new SpatialReference(null); sr4326.ImportFromEPSG(4326); sr4490 = new SpatialReference(null); sr4490.ImportFromEPSG(4490); string wkt = "GEOGCS[\"GCS_Moon_2000\",\r\n" + " DATUM[\"D_Moon_2000\",\r\n" + " SPHEROID[\"Moon_2000_IAU_IAG\",1737400,0,\r\n" + " AUTHORITY[\"ESRI\",\"107903\"]],\r\n" + " AUTHORITY[\"ESRI\",\"106903\"]],\r\n" + " PRIMEM[\"Reference_Meridian\",0,\r\n" + " AUTHORITY[\"ESRI\",\"108900\"]],\r\n" + " UNIT[\"degree\",0.0174532925199433,\r\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\r\n" + " AUTHORITY[\"ESRI\",\"104903\"]]"; sr104903 = new SpatialReference(wkt); //sr104903.ImportFromWkt(ref wkt); isInited = true; } } } /// /// 实例 /// public static GdalHelper Instance { get { lock (obj) { if (null == instance) instance = new GdalHelper(); return instance; } } } /// /// 注册GDAL /// public void RegisterGDal() { string gdalData = Path.Combine(Tools.BaseDir, "gdal-data"); Environment.SetEnvironmentVariable("GDAL_DATA", gdalData); string proj7 = Path.Combine(Tools.BaseDir, "proj7\\share"); Environment.SetEnvironmentVariable("PROJ_LIB", proj7); Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // NO,YES Gdal.SetConfigOption("SHAPE_ENCODING", ""); // 空,gb2312,CP936 Ogr.RegisterAll(); Gdal.AllRegister(); } #endregion #region 读取栅格信息 /// /// 读取栅格信息 /// public void ReadRasterInfo(ViewData vd) { Dataset ds = null; try { ds = Gdal.Open(vd.FilePath, Access.GA_ReadOnly); if (null == ds || ds.RasterCount == 0) return; vd.Meta.gridsize = string.Format("{0},{1}", ds.RasterXSize, ds.RasterYSize); // 行列数 SpatialReference sr = ds.GetSpatialRef(); if (sr != null) { vd.Meta.coor_sys = sr.GetName(); // 坐标系统 if (MOON200 == vd.Meta.coor_sys) { vd.Meta.epsg = "ESRI:104903"; // EPSG编码 } else { string code = sr.GetAuthorityCode(null); vd.Meta.epsg = string.IsNullOrEmpty(code) ? null : "EPSG:" + code; // EPSG编码 } } vd.Meta.bands = ds.RasterCount.ToString(); // 波段数 vd.Meta.band_type = Enum.GetName(typeof(DataType), ds.GetRasterBand(1).DataType); // 数据类型 ColorTable ct = ds.GetRasterBand(1).GetRasterColorTable(); vd.Meta.ct = null == ct ? null : ct.ToString(); // 数据颜色表 vd.Meta.h_datum = null; // 高程基准 double[] transform = new double[6]; ds.GetGeoTransform(transform); vd.Meta.resolution = string.Format("{0},{1}", transform[1], transform[5]); // 分辨率 if (!EPSGS.Contains(vd.Meta.epsg)) return; Geometry minPoint = GetMinPoint(ds); Geometry maxPoint = GetMaxPoint(ds); double xmin = minPoint.GetX(0); double ymin = minPoint.GetY(0); double xmax = maxPoint.GetX(0); double ymax = maxPoint.GetY(0); string wkt = string.Format("POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{0} {1}))", xmin, ymax, xmax, ymax, xmax, ymin, xmin, ymin); vd.Meta.geom = wkt; // 四至范围 } catch (Exception ex) { LogOut.Error(ex.StackTrace); } finally { if (ds != null) ds.Dispose(); } } /// /// 获取Dataset的最小点 /// private Geometry GetMinPoint(Dataset ds) { double[] transform = new double[6]; ds.GetGeoTransform(transform); string epsg = ds.GetSpatialRef().GetAuthorityCode(null); double xMin = transform[0]; double yMin = transform[3] - ds.RasterYSize * transform[1]; Geometry point = new Geometry(wkbGeometryType.wkbPoint); point.AddPoint(xMin, yMin, 0); return Transform(ds, point, epsg); } /// /// 获取Dataset的最大点 /// private Geometry GetMaxPoint(Dataset ds) { double[] transform = new double[6]; ds.GetGeoTransform(transform); string epsg = ds.GetSpatialRef().GetAuthorityCode(null); double xMax = transform[0] + (ds.RasterXSize * transform[1]); double yMax = transform[3]; Geometry point = new Geometry(wkbGeometryType.wkbPoint); point.AddPoint(xMax, yMax, 0); return Transform(ds, point, epsg); } /// /// 坐标转换 /// private Geometry Transform(Dataset ds, Geometry point, string epsg) { if (string.IsNullOrEmpty(epsg)) { point.AssignSpatialReference(sr104903); return point; } if ("4326" == epsg) { point.AssignSpatialReference(sr4326); return point; } if ("4490" == epsg) { point.AssignSpatialReference(sr4490); return point; } point.AssignSpatialReference(ds.GetSpatialRef()); if (ds.GetSpatialRef().GetName().Contains("CGCS2000")) { point.TransformTo(sr4490); } else { point.TransformTo(sr4326); } return point; } /// /// 创建多边形 /// public Geometry CreatePolygon(double xmin, double xmax, double ymin, double ymax, SpatialReference sr) { string kwt = string.Format("POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{0} {1}))", xmin, ymax, xmax, ymax, xmax, ymin, xmin, ymin); Geometry geo = Geometry.CreateFromWkt(kwt); geo.AssignSpatialReference(sr); return geo; } #endregion } }