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;
|
import java.util.ArrayList;
|
import java.util.List;
|
|
/**
|
* 读取栅格服务
|
* @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);
|
// 设置最值
|
setMinAndMax(ds, mf);
|
|
// 分辨率
|
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;
|
}
|
|
/**
|
* 设置最值
|
* GDALRasterBand::GetHistogram 统计直方图
|
*/
|
private void setMinAndMax(Dataset ds, MetaFileEntity mf) {
|
List<Double> minList = new ArrayList<>();
|
List<Double> maxList = new ArrayList<>();
|
|
for (int i = 1; i <= ds.getRasterCount(); i++) {
|
Double[] min = new Double[1];
|
Double[] max = new Double[1];
|
ds.GetRasterBand(i).GetMinimum(min);
|
ds.GetRasterBand(i).GetMaximum(max);
|
|
minList.add(min[0]);
|
maxList.add(max[0]);
|
}
|
|
mf.setMin(StringHelper.join(minList, ","));
|
mf.setMax(StringHelper.join(maxList, ","));
|
}
|
}
|
