package com.moon.server.service.data; import com.moon.server.entity.all.StaticData; import com.moon.server.entity.data.MetaFileEntity; import com.moon.server.helper.StringHelper; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; import org.gdal.gdal.ColorTable; import org.gdal.gdal.Dataset; import org.gdal.gdal.gdal; import org.gdal.gdalconst.gdalconst; import org.gdal.ogr.Geometry; import org.gdal.ogr.ogr; import org.gdal.osr.SpatialReference; import org.springframework.stereotype.Service; import java.io.File; /** * 读取栅格服务 * @author WWW * @date 2023-08-25 */ @Service public class ReadRasterService { private SpatialReference sr4326; private SpatialReference sr4490; private SpatialReference sr104903; private final static Log log = LogFactory.getLog(ReadRasterService.class); /** * 初始化空间引用 */ private void initSpatialReference() { if (null == sr4326) { sr4326 = new SpatialReference(); sr4326.ImportFromEPSG(StaticData.I4326); sr4490 = new SpatialReference(); sr4490.ImportFromEPSG(StaticData.I4490); sr104903 = new SpatialReference(StaticData.MOON_2000_WKT); } } /** * 读取栅格信息 */ public void readRasterInfo(MetaFileEntity mf, String file) { Dataset ds = null; try { File f = new File(file); if (!f.exists() || f.isDirectory()) { return; } ds = gdal.Open(file, gdalconst.GA_ReadOnly); if (null == ds || ds.getRasterCount() == 0 || null == ds.GetSpatialRef()) { return; } SpatialReference sr = ds.GetSpatialRef(); // 坐标系统 mf.setCoorSys(sr.GetName()); if (StaticData.MOON200.equals(mf.getCoorSys())) { // EPSG编码 mf.setEpsg("ESRI:" + StaticData.I104903); } else { // EPSG编码:PROJCS、GEOGCS、GEOGCS|UNIT 或 NULL String code = sr.GetAuthorityCode(null); mf.setEpsg(StringHelper.isEmpty(code) ? null : "EPSG:" + code); } // 行列数 mf.setGridsize(String.format("%d,%d", ds.getRasterXSize(), ds.getRasterYSize())); // 波段数 mf.setBands(String.valueOf(ds.getRasterCount())); // 数据类型 int dataType = ds.GetRasterBand(1).GetRasterDataType(); mf.setBandType(getDataType(dataType)); // 数据颜色表 ColorTable colorTable = ds.GetRasterBand(1).GetRasterColorTable(); mf.setCt(null == colorTable ? null : colorTable.toString()); // 高程基准 mf.sethDatum(null); // 分辨率 double[] tr = new double[6]; ds.GetGeoTransform(tr); mf.setResolution(String.format("%f,%f", tr[1], Math.abs(tr[5]))); if (tr[StaticData.I0] == 0.0 && tr[StaticData.I1] == 1.0 && tr[StaticData.I2] == 0.0 && tr[StaticData.I3] == 0.0 && tr[StaticData.I4] == 0.0 && tr[StaticData.I5] == 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 ((%f %f,%f %f,%f %f,%f %f,%f %f))')", xmin, ymax, xmax, ymax, xmax, ymin, xmin, ymin, xmin, ymax); mf.setGeom(geom); } catch (Exception ex) { log.error(ex.getMessage(), ex); } finally { if (null != ds) { ds.delete(); } } } /** * 获取数据类型 */ private String getDataType(int dataType) { if (dataType == gdalconst.GDT_Byte) { return "GDT_Byte"; } if (dataType == gdalconst.GDT_UInt16) { return "GDT_UInt16"; } if (dataType == gdalconst.GDT_Int16) { return "GDT_Int16"; } if (dataType == gdalconst.GDT_UInt32) { return "GDT_UInt32"; } if (dataType == gdalconst.GDT_Int32) { return "GDT_Int32"; } if (dataType == gdalconst.GDT_Float32) { return "GDT_Float32"; } if (dataType == gdalconst.GDT_Float64) { return "GDT_Float64"; } if (dataType == gdalconst.GDT_CInt16) { return "GDT_CInt16"; } if (dataType == gdalconst.GDT_CInt32) { return "GDT_CInt32"; } if (dataType == gdalconst.GDT_CFloat32) { return "GDT_CFloat32"; } if (dataType == gdalconst.GDT_CFloat64) { return "GDT_CFloat64"; } return null; } /** * 获取Dataset的最小点 */ private Geometry getMinPoint(Dataset ds) { double[] transform = new double[6]; ds.GetGeoTransform(transform); double xMin = transform[0]; double yMin = transform[3] - ds.getRasterYSize() * transform[1]; Geometry point = new Geometry(ogr.wkbPoint); point.AddPoint(xMin, yMin, 0); return transform(ds, point); } /** * 获取Dataset的最大点 */ 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.getRasterXSize() * transform[1]); double yMax = transform[3]; Geometry point = new Geometry(ogr.wkbPoint); point.AddPoint(xMax, yMax, 0); return transform(ds, point); } /** * 坐标转换 */ private Geometry transform(Dataset ds, Geometry point) { this.initSpatialReference(); point.AssignSpatialReference(ds.GetSpatialRef()); if (ds.GetSpatialRef().IsGeographic() > 0) { return point; } String srsName = ds.GetSpatialRef().GetName(); if (srsName.contains(StaticData.CGCS2000)) { point.TransformTo(sr4490); } else { point.TransformTo(sr4326); } point.SwapXY(); return point; } }