using OSGeo.GDAL;
|
using OSGeo.OGR;
|
using OSGeo.OSR;
|
using System;
|
using System.Collections.Generic;
|
using System.Diagnostics;
|
using System.IO;
|
using System.Linq;
|
using System.Threading.Tasks;
|
using System.Windows;
|
using System.Windows.Markup;
|
|
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 sr104903 = null;
|
|
public static string MOON200 = "GCS_Moon_2000";
|
|
private GdalHelper()
|
{
|
lock (obj)
|
{
|
if (!isInited)
|
{
|
RegisterGDal();
|
|
sr4326 = new SpatialReference(null);
|
sr4326.ImportFromEPSG(4326);
|
|
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;
|
}
|
}
|
}
|
|
public void RegisterGDal()
|
{
|
string gdalData = Path.Combine(StaticData.BasePath, "gdal-data");
|
Environment.SetEnvironmentVariable("GDAL_DATA", gdalData);
|
|
string proj7 = Path.Combine(StaticData.BasePath, "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
|
|
public void ReadTiff()
|
{
|
string outPath = "D:\\xyz\\ce\\xyz";
|
string tif = "D:\\Moon\\data\\dom_tif\\moon.tif";
|
|
Dataset ds = Gdal.Open(tif, Access.GA_ReadOnly);
|
int x = ds.RasterXSize;
|
int y = ds.RasterYSize;
|
|
SpatialReference sr = ds.GetSpatialRef();
|
string srName = sr.GetName(); // GCS_Moon_2000
|
// string code = sr4326.GetAuthorityCode(null); // 4326
|
if (MOON200 != srName) return;
|
|
double[] lt = ImageToGeoSpace(ds, 0, 0); // -179.99998351102391, 89.999983511056882
|
double[] rb = ImageToGeoSpace(ds, x, y); // 179.999788852709, -89.9999026708093
|
|
string polygon = string.Format("", "");
|
}
|
|
// 从地理空间转换到像素空间
|
private int[] Geo2ImageSpace(Dataset ds, double x, double y)
|
{
|
double[] transform = new double[6];
|
ds.GetGeoTransform(transform); // 影像坐标变换参数
|
|
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 };
|
}
|
|
// 从地理空间转换到像素空间
|
private int[] Geo2ImageSpace(double[] transform, double x, double y)
|
{
|
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 };
|
}
|
|
// 从像素空间转换到地理空间
|
private double[] ImageToGeoSpace(Dataset ds, int row, int col)
|
{
|
double[] transform = new double[6];
|
ds.GetGeoTransform(transform);
|
|
double x = transform[0] + row * transform[1] + col * transform[2];
|
double y = transform[3] + row * transform[4] + col * transform[5];
|
|
return new double[] { x, y };
|
}
|
|
#region GDAL切图
|
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];
|
|
Geometry point = new Geometry(wkbGeometryType.wkbPoint);
|
point.AddPoint(xMin, yMin, 0);
|
if ("4326" == epsg)
|
{
|
point.AssignSpatialReference(sr4326);
|
return point;
|
}
|
|
point.AssignSpatialReference(ds.GetSpatialRef());
|
point.TransformTo(sr4326);
|
return point;
|
}
|
|
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] + (ds.RasterYSize * transform[1]);
|
|
Geometry point = new Geometry(wkbGeometryType.wkbPoint);
|
point.AddPoint(xMax, yMax, 0);
|
if ("4326" == epsg)
|
{
|
point.AssignSpatialReference(sr4326);
|
return point;
|
}
|
|
point.AssignSpatialReference(ds.GetSpatialRef());
|
point.TransformTo(sr4326);
|
return point;
|
}
|
|
public void GenerateTiles()
|
{
|
string outPath = "D:\\xyz\\ce\\xyz2";
|
string file = "D:\\xyz\\ce\\5_A1.tif";
|
if (!File.Exists(file)) return;
|
|
Stopwatch sw = new Stopwatch();
|
sw.Start(); // 程序开始时间
|
Dataset ds = null;
|
try
|
{
|
ds = Gdal.Open(file, Access.GA_Update);
|
if (null == ds) return;
|
|
double[] transform = new double[6];
|
ds.GetGeoTransform(transform);
|
int rasterCount = ds.RasterCount;
|
Console.WriteLine(string.Format("Origin = ({0}, {1})", transform[0], transform[3]));
|
Console.WriteLine(string.Format("Pixel Size = ({0}, {1})", transform[1], transform[5]));
|
|
// 4.1 首先获取原始影像的地理坐标范围
|
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);
|
Geometry imageBound = CreatePolygon(xmin, xmax, ymin, ymax, sr4326);
|
|
// 4.2 获取原始影像的像素分辨率
|
// 原始图像东西方向像素分辨率
|
double src_w_e_pixel_resolution = (xmax - xmin) / ds.RasterXSize;
|
// 原始图像南北方向像素分辨率
|
double src_n_s_pixel_resolution = (ymax - ymin) / ds.RasterYSize;
|
|
// 获取Band
|
Band in_band1 = ds.GetRasterBand(1);
|
Band in_band2 = ds.GetRasterBand(2);
|
Band in_band3 = ds.GetRasterBand(3);
|
//in_band1.Fill(0, 255); // GdalConst.GMF_NODATA
|
in_band1.SetNoDataValue(0); // -9999
|
in_band2.SetNoDataValue(0);
|
in_band3.SetNoDataValue(0);
|
|
int BufferSize = 256 * 256 * 5; // GdalConst.GDT_Int32
|
for (int zoom = 10; zoom <= 18; zoom++)
|
{
|
// 4.3 根据原始影像地理范围求解切片行列号 // 经纬度转瓦片编号
|
int tileRowMax = Lat2tile(ymin, zoom); // 纬度 -90 -> 90 lat
|
int tileRowMin = Lat2tile(ymax, zoom);
|
int tileColMin = Lon2tile(xmin, zoom); // 经度 -180 -> 180 lon
|
int tileColMax = Lon2tile(xmax, zoom);
|
|
Parallel.For(tileColMin, tileColMax + 1, (col) =>
|
{
|
Parallel.For(tileRowMin, tileRowMax + 1, (row) =>
|
{
|
// 4.4 求原始影像地理范围与指定缩放级别指定行列号的切片交集
|
double tempLatMin = Tile2lat(row + 1, zoom);
|
double tempLatMax = Tile2lat(row, zoom);
|
double tempLonMin = Tile2lon(col, zoom);
|
double tempLonMax = Tile2lon(col + 1, zoom);
|
|
Console.WriteLine(string.Format("{0}\\{1}\\{2}.png", zoom, col, row));
|
Geometry tileBound = CreatePolygon(tempLonMin, tempLonMax, tempLatMin, tempLatMax, sr4326);
|
Geometry intersect = tileBound.Intersection(imageBound);
|
if (null == intersect)
|
{
|
Console.WriteLine(string.Format("{0}\\{1}\\{2}.png，不存在", zoom, col, row));
|
return;
|
}
|
Envelope env = new Envelope();
|
intersect.GetEnvelope(env);
|
|
// 4.5 求解当前切片的像素分辨率(默认切片大小为256*256)
|
// 切片东西方向像素分辨率
|
double dst_w_e_pixel_resolution = (tempLonMax - tempLonMin) / 256;
|
// 切片南北方向像素分辨率
|
double dst_n_s_pixel_resolution = (tempLatMax - tempLatMin) / 256;
|
|
// 4.6 计算交集的像素信息
|
// 求切图范围和原始图像交集的起始点像素坐标
|
int offset_x = (int)((env.MinX - xmin) / src_w_e_pixel_resolution);
|
int offset_y = (int)Math.Abs((env.MaxY - ymax) / src_n_s_pixel_resolution);
|
|
// 求在切图地理范围内的原始图像的像素大小
|
int block_xsize = (int)((env.MaxX - env.MinX) / src_w_e_pixel_resolution);
|
int block_ysize = (int)((env.MaxY - env.MinY) / src_n_s_pixel_resolution);
|
|
// 求原始图像在切片内的像素大小
|
int image_Xbuf = (int)Math.Ceiling((env.MaxX - env.MinX) / dst_w_e_pixel_resolution);
|
int image_Ybuf = (int)Math.Ceiling(Math.Abs((env.MaxY - env.MinY) / dst_n_s_pixel_resolution));
|
|
// 求原始图像在切片中的偏移坐标
|
int imageOffsetX = (int)((env.MinX - tempLonMin) / dst_w_e_pixel_resolution);
|
int imageOffsetY = (int)Math.Abs((env.MaxY - tempLatMax) / dst_n_s_pixel_resolution);
|
imageOffsetX = imageOffsetX > 0 ? imageOffsetX : 0;
|
imageOffsetY = imageOffsetY > 0 ? imageOffsetY : 0;
|
|
// 4.7 使用GDAL的ReadRaster方法对影像指定范围进行读取与压缩。
|
// 推荐在切片前建立原始影像的金字塔文件，ReadRaster在内部实现中可直接读取相应级别的金字塔文件，提高效率。
|
int[] band1BuffData = new int[BufferSize]; // 256 * 256 * GdalConst.GDT_Int32
|
int[] band2BuffData = new int[BufferSize];
|
int[] band3BuffData = new int[BufferSize];
|
|
try
|
{
|
in_band1.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
in_band2.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
in_band3.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
|
// 4.8 将切片数据写入文件
|
// 使用gdal的MEM驱动在内存中创建一块区域存储图像数组
|
OSGeo.GDAL.Driver memDriver = Gdal.GetDriverByName("MEM");
|
Dataset msmDS = memDriver.Create("msmDS", 256, 256, 4, DataType.GDT_Int32, null);
|
Band dstBand1 = msmDS.GetRasterBand(1);
|
Band dstBand2 = msmDS.GetRasterBand(2);
|
Band dstBand3 = msmDS.GetRasterBand(3);
|
|
// 设置alpha波段数据,实现背景透明
|
Band alphaBand = msmDS.GetRasterBand(4);
|
int[] alphaData = new int[BufferSize];
|
for (int index = 0; index < alphaData.Length; index++)
|
{
|
if (band1BuffData[index] > 0) alphaData[index] = 255;
|
}
|
|
// 写各个波段数据
|
dstBand1.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
dstBand2.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
dstBand3.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
alphaBand.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, alphaData, image_Xbuf, image_Ybuf, 0, 0);
|
|
string path = Path.Combine(outPath, zoom.ToString(), col.ToString());
|
if (!Directory.Exists(path)) Directory.CreateDirectory(path);
|
|
// 使用PNG驱动将内存中的图像数组写入文件
|
string pngPath = path + "\\" + row + ".png";
|
OSGeo.GDAL.Driver pngDriver = Gdal.GetDriverByName("PNG");
|
Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS, 0, null, null, null);
|
|
msmDS.FlushCache();
|
pngDs.Dispose(); // 释放内存
|
}
|
catch (Exception ex)
|
{
|
Console.WriteLine(ex.Message);
|
}
|
});
|
});
|
}
|
}
|
catch (Exception ex)
|
{
|
Console.WriteLine(ex.Message);
|
}
|
finally
|
{
|
if (null != ds) ds.Dispose();
|
sw.Stop(); // 程序结束
|
Console.WriteLine("耗时" + Math.Round(sw.ElapsedMilliseconds / 1000.0, 2) + "s"); // 秒
|
}
|
}
|
|
public void GenerateTiles2()
|
{
|
string outPath = "D:\\xyz\\ce\\xyz";
|
string file = "D:\\xyz\\ce\\5_A1.tif";
|
//string file = "d:\\xyz\\dem\\dem\\33b.tif"; //string file = "D:\\xyz\\dq\\dq.vrt";
|
if (!File.Exists(file)) return;
|
|
Stopwatch sw = new Stopwatch();
|
sw.Start(); // 程序开始时间
|
Dataset ds = null;
|
try
|
{
|
ds = Gdal.Open(file, Access.GA_Update);
|
if (null == ds) return;
|
|
double[] transform = new double[6];
|
ds.GetGeoTransform(transform);
|
int rasterCount = ds.RasterCount;
|
Console.WriteLine(string.Format("Origin = ({0}, {1})", transform[0], transform[3]));
|
Console.WriteLine(string.Format("Pixel Size = ({0}, {1})", transform[1], transform[5]));
|
|
// 4.1 首先获取原始影像的地理坐标范围
|
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);
|
Geometry imageBound = CreatePolygon(xmin, xmax, ymin, ymax, sr4326);
|
|
// 4.2 获取原始影像的像素分辨率
|
// 原始图像东西方向像素分辨率
|
double src_w_e_pixel_resolution = (xmax - xmin) / ds.RasterXSize;
|
// 原始图像南北方向像素分辨率
|
double src_n_s_pixel_resolution = (ymax - ymin) / ds.RasterYSize;
|
|
// 获取Band
|
Band in_band1 = ds.GetRasterBand(1);
|
Band in_band2 = ds.GetRasterBand(2);
|
Band in_band3 = ds.GetRasterBand(3);
|
//in_band1.Fill(0, 255); // GdalConst.GMF_NODATA
|
in_band1.SetNoDataValue(0); // -9999
|
in_band2.SetNoDataValue(0);
|
in_band3.SetNoDataValue(0);
|
|
int BufferSize = 256 * 256 * 5; // GdalConst.GDT_Int32
|
for (int zoom = 10; zoom <= 18; zoom++)
|
{
|
// 4.3 根据原始影像地理范围求解切片行列号 // 经纬度转瓦片编号
|
int tileRowMax = Lat2tile(ymin, zoom); // 纬度 -90 -> 90 lat
|
int tileRowMin = Lat2tile(ymax, zoom);
|
int tileColMin = Lon2tile(xmin, zoom); // 经度 -180 -> 180 lon
|
int tileColMax = Lon2tile(xmax, zoom);
|
//int tileRowMax = lat2tile(latMin, zoom); // 纬度 -90 ——90 lat
|
//int tileRowMin = lat2tile(latMax, zoom); // 经度 -180 -- 180 lon
|
//int tileColMin = lon2tile(lonMin, zoom);
|
//int tileColMax = lon2tile(lonMax, zoom);
|
|
for (int col = tileColMin; col <= tileColMax; col++)
|
{
|
for (int row = tileRowMin; row <= tileRowMax; row++)
|
{
|
// 4.4 求原始影像地理范围与指定缩放级别指定行列号的切片交集
|
double tempLatMin = Tile2lat(row + 1, zoom);
|
double tempLatMax = Tile2lat(row, zoom);
|
double tempLonMin = Tile2lon(col, zoom);
|
double tempLonMax = Tile2lon(col + 1, zoom);
|
//double tempLatMin = tile2lat(row + 1, zoom);
|
//double tempLatMax = tile2lat(row, zoom);
|
//double tempLonMin = tile2lon(col, zoom);
|
//double tempLonMax = tile2lon(col + 1, zoom);
|
|
Console.WriteLine(string.Format("{0}\\{1}\\{2}.png", zoom, col, row));
|
Geometry tileBound = CreatePolygon(tempLonMin, tempLonMax, tempLatMin, tempLatMax, sr4326);
|
Geometry intersect = tileBound.Intersection(imageBound);
|
if (null == intersect)
|
{
|
Console.WriteLine(string.Format("{0}\\{1}\\{2}.png，不存在", zoom, col, row));
|
continue;
|
}
|
Envelope env = new Envelope();
|
intersect.GetEnvelope(env);
|
|
// 4.5 求解当前切片的像素分辨率(默认切片大小为256*256)
|
// 切片东西方向像素分辨率
|
double dst_w_e_pixel_resolution = (tempLonMax - tempLonMin) / 256;
|
// 切片南北方向像素分辨率
|
double dst_n_s_pixel_resolution = (tempLatMax - tempLatMin) / 256;
|
|
// 4.6 计算交集的像素信息
|
// 求切图范围和原始图像交集的起始点像素坐标
|
int offset_x = (int)((env.MinX - xmin) / src_w_e_pixel_resolution);
|
int offset_y = (int)Math.Abs((env.MaxY - ymax) / src_n_s_pixel_resolution);
|
|
// 求在切图地理范围内的原始图像的像素大小
|
int block_xsize = (int)((env.MaxX - env.MinX) / src_w_e_pixel_resolution);
|
int block_ysize = (int)((env.MaxY - env.MinY) / src_n_s_pixel_resolution);
|
|
// 求原始图像在切片内的像素大小
|
int image_Xbuf = (int)Math.Ceiling((env.MaxX - env.MinX) / dst_w_e_pixel_resolution);
|
int image_Ybuf = (int)Math.Ceiling(Math.Abs((env.MaxY - env.MinY) / dst_n_s_pixel_resolution));
|
|
// 求原始图像在切片中的偏移坐标
|
int imageOffsetX = (int)((env.MinX - tempLonMin) / dst_w_e_pixel_resolution);
|
int imageOffsetY = (int)Math.Abs((env.MaxY - tempLatMax) / dst_n_s_pixel_resolution);
|
imageOffsetX = imageOffsetX > 0 ? imageOffsetX : 0;
|
imageOffsetY = imageOffsetY > 0 ? imageOffsetY : 0;
|
|
// 4.7 使用GDAL的ReadRaster方法对影像指定范围进行读取与压缩。
|
// 推荐在切片前建立原始影像的金字塔文件，ReadRaster在内部实现中可直接读取相应级别的金字塔文件，提高效率。
|
int[] band1BuffData = new int[BufferSize]; // 256 * 256 * GdalConst.GDT_Int32
|
int[] band2BuffData = new int[BufferSize];
|
int[] band3BuffData = new int[BufferSize];
|
|
try
|
{
|
// ReadRaster(int xOff, int yOff, int xSize, int ySize, int[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace)
|
//in_band1.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf, GdalConst.GDT_Int32, band1BuffData, 0, 0);
|
//in_band2.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf, GdalConst.GDT_Int32, band2BuffData, 0, 0);
|
//in_band3.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf, GdalConst.GDT_Int32, band3BuffData, 0, 0);
|
in_band1.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
in_band2.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
in_band3.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
|
// 4.8 将切片数据写入文件
|
// 使用gdal的MEM驱动在内存中创建一块区域存储图像数组
|
OSGeo.GDAL.Driver memDriver = Gdal.GetDriverByName("MEM");
|
Dataset msmDS = memDriver.Create("msmDS", 256, 256, 4, DataType.GDT_Int32, null);
|
Band dstBand1 = msmDS.GetRasterBand(1);
|
Band dstBand2 = msmDS.GetRasterBand(2);
|
Band dstBand3 = msmDS.GetRasterBand(3);
|
|
// 设置alpha波段数据,实现背景透明
|
Band alphaBand = msmDS.GetRasterBand(4);
|
int[] alphaData = new int[BufferSize];
|
for (int index = 0; index < alphaData.Length; index++)
|
{
|
if (band1BuffData[index] > 0)
|
{
|
alphaData[index] = 255;
|
}
|
}
|
|
// 写各个波段数据
|
// WriteRaster(int xOff, int yOff, int xSize, int ySize, int[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace)
|
//dstBand1.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band1BuffData);
|
//dstBand2.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band2BuffData);
|
//dstBand3.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band3BuffData);
|
//alphaBand.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, alphaData);
|
dstBand1.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
dstBand2.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
dstBand3.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);
|
alphaBand.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, alphaData, image_Xbuf, image_Ybuf, 0, 0);
|
|
//String pngPath = "D:\\xyz\\temp" + "\\" + zoom + "c" + col + "r" + row + ".png";
|
//Console.WriteLine("pngPath=" + pngPath);
|
string path = Path.Combine(outPath, zoom.ToString(), col.ToString());
|
if (!Directory.Exists(path)) Directory.CreateDirectory(path);
|
|
// 使用PNG驱动将内存中的图像数组写入文件
|
string pngPath = path + "\\" + row + ".png";
|
OSGeo.GDAL.Driver pngDriver = Gdal.GetDriverByName("PNG");
|
Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS, 0, null, null, null);
|
|
// 释放内存
|
msmDS.FlushCache();
|
pngDs.Dispose();
|
}
|
catch (Exception ex)
|
{
|
Console.WriteLine(ex.Message);
|
}
|
}
|
}
|
}
|
}
|
catch (Exception ex)
|
{
|
Console.WriteLine(ex.Message);
|
}
|
finally
|
{
|
if (null != ds) ds.Dispose();
|
sw.Stop(); // 程序结束
|
Console.WriteLine("耗时" + Math.Round(sw.ElapsedMilliseconds / 1000.0, 2) + "s"); // 秒
|
}
|
}
|
|
/// <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>
|
/// 经度转瓦片编号
|
/// </summary>
|
public static int Lon2tile(double lon, int zoom)
|
{
|
return (int)(Math.Floor((lon + 180) / 360 * Math.Pow(2, zoom)));
|
}
|
|
/// <summary>
|
/// 纬度转瓦片编号
|
/// </summary>
|
public static int Lat2tile(double lat, int zoom)
|
{
|
return (int)(Math.Floor((1 - Math.Log(Math.Tan(lat * Math.PI / 180) + 1 / Math.Cos(lat * Math.PI / 180)) / Math.PI) / 2 * Math.Pow(2, zoom)));
|
}
|
|
/// <summary>
|
/// 瓦片编号转经度
|
/// </summary>
|
public static double Tile2lon(int col, int zoom)
|
{
|
return col / Math.Pow(2.0, zoom) * 360.0 - 180;
|
}
|
|
/// <summary>
|
/// 瓦片编号转纬度
|
/// </summary>
|
public static double Tile2lat(int row, int zoom)
|
{
|
double n = Math.PI - (2.0 * Math.PI * row) / Math.Pow(2.0, zoom);
|
return ToDegrees(Math.Atan(Math.Sinh(n)));
|
}
|
|
/// <summary>
|
/// 常数e
|
/// </summary>
|
public static readonly double E = 2.7182818284590452354;
|
|
/// <summary>
|
/// 常数Pi
|
/// </summary>
|
public static readonly double PI = 3.14159265358979323846;
|
|
/// <summary>
|
/// 度转弧度
|
/// </summary>
|
public static double ToRadians(double angdeg)
|
{
|
return angdeg / 180.0 * Math.PI;
|
}
|
|
/// <summary>
|
/// 弧度转度
|
/// </summary>
|
public static double ToDegrees(double angrad)
|
{
|
return angrad * 180.0 / Math.PI;
|
}
|
#endregion
|
}
|
}
|
