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<string> EPSGS = new List<string>() { "EPSG:4326", "EPSG:4490", "ESRI:104903" };
|
|
/// <summary>
|
/// 构造函数
|
/// </summary>
|
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;
|
}
|
}
|
}
|
|
/// <summary>
|
/// 实例
|
/// </summary>
|
public static GdalHelper Instance
|
{
|
get
|
{
|
lock (obj)
|
{
|
if (null == instance) instance = new GdalHelper();
|
|
return instance;
|
}
|
}
|
}
|
|
/// <summary>
|
/// 注册GDAL
|
/// </summary>
|
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 读取栅格信息
|
/// <summary>
|
/// 读取栅格信息
|
/// </summary>
|
public void ReadRasterInfo(ViewData vd)
|
{
|
Dataset ds = null;
|
try
|
{
|
ds = Gdal.Open(vd.FilePath, Access.GA_ReadOnly);
|
if (null == ds || ds.RasterCount == 0 || null == ds.GetSpatialRef()) return;
|
|
vd.Meta.gridsize = string.Format("{0},{1}", ds.RasterXSize, ds.RasterYSize); // 行列数
|
|
SpatialReference sr = ds.GetSpatialRef();
|
vd.Meta.coor_sys = sr.GetName(); // 坐标系统
|
if (MOON200 == vd.Meta.coor_sys)
|
{
|
vd.Meta.epsg = "ESRI:104903"; // EPSG编码
|
}
|
else
|
{
|
string code = sr.GetAuthorityCode(null); // PROJCS、GEOGCS、GEOGCS|UNIT 或 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; // 高程基准
|
setMinAndMax(ds, vd.Meta); // 设置最值
|
|
double[] tr = new double[6];
|
ds.GetGeoTransform(tr);
|
vd.Meta.resolution = string.Format("{0},{1}", tr[1], Math.Abs(tr[5])); // 分辨率
|
|
if (tr[0] == 0.0 && tr[1] == 1.0 && tr[2] == 0.0 && tr[3] == 0.0 && tr[4] == 0.0 && tr[5] == 1.0) 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 geom = string.Format("ST_GeomFromText('POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{0} {1}))')", xmin, ymax, xmax, ymax, xmax, ymin, xmin, ymin);
|
|
vd.Meta.geom = geom; // 四至范围
|
}
|
catch (Exception ex)
|
{
|
LogOut.Error(ex.StackTrace);
|
}
|
finally
|
{
|
if (ds != null) ds.Dispose();
|
}
|
}
|
|
/// <summary>
|
/// 获取Dataset的最小点
|
/// </summary>
|
private Geometry GetMinPoint(Dataset ds)
|
{
|
double[] transform = new double[6];
|
ds.GetGeoTransform(transform);
|
|
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);
|
}
|
|
/// <summary>
|
/// 获取Dataset的最大点
|
/// </summary>
|
private Geometry GetMaxPoint(Dataset ds)
|
{
|
/**
|
* transform[0] 左上角x坐标
|
* transform[1] 东西方向分辨率
|
* transform[2] 旋转角度, 0表示图像 "北方朝上"
|
*
|
* transform[3] 左上角y坐标
|
* transform[4] 旋转角度, 0表示图像 "北方朝上"
|
* transform[5] 南北方向分辨率
|
*/
|
double[] transform = new double[6];
|
ds.GetGeoTransform(transform);
|
|
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);
|
}
|
|
/// <summary>
|
/// 坐标转换
|
/// </summary>
|
private Geometry Transform(Dataset ds, Geometry point)
|
{
|
point.AssignSpatialReference(ds.GetSpatialRef());
|
if (ds.GetSpatialRef().IsGeographic() > 0)
|
{
|
return point;
|
}
|
|
string srsName = ds.GetSpatialRef().GetName();
|
if (srsName.Contains("CGCS2000"))
|
{
|
point.TransformTo(sr4490);
|
}
|
else
|
{
|
point.TransformTo(sr4326);
|
}
|
point.SwapXY();
|
|
return point;
|
}
|
|
/// <summary>
|
/// 创建多边形
|
/// </summary>
|
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;
|
}
|
|
/// <summary>
|
/// 设置最值
|
/// GDALRasterBand::GetHistogram 统计直方图
|
/// </summary>
|
private void setMinAndMax(Dataset ds, SysMeta meta)
|
{
|
List<double> minList = new List<double>();
|
List<double> maxList = new List<double>();
|
|
for (int i = 1; i <= ds.RasterCount; i++)
|
{
|
int hasval;
|
double min, max;
|
ds.GetRasterBand(i).GetMinimum(out min, out hasval);
|
ds.GetRasterBand(i).GetMaximum(out max, out hasval);
|
|
minList.Add(min);
|
maxList.Add(max);
|
}
|
|
meta.min = string.Join(",", minList.ToArray());
|
meta.max = string.Join(",", maxList.ToArray());
|
}
|
#endregion
|
|
public void CsTransform(double x, double y, int epsg)
|
{
|
SpatialReference srs = new SpatialReference(null);
|
srs.ImportFromEPSG(epsg);
|
|
Geometry point = new Geometry(wkbGeometryType.wkbPoint);
|
point.AddPoint(x, y, 0);
|
point.AssignSpatialReference(srs);
|
|
if (srs.GetName().Contains("CGCS2000"))
|
{
|
point.TransformTo(sr4490);
|
}
|
else
|
{
|
point.TransformTo(sr4326);
|
}
|
point.SwapXY();
|
|
double[] xy = new double[] { point.GetX(0), point.GetY(0) };
|
}
|
}
|
}
|
