前台是openlayer vue 项目,后台是spring boot 。需求是前台页面通过设置空间范围,然后对范围内的影像,进行分类输出。中间踩过一些坑,将新路历程记录下来,一来是回顾记录下来,二来,也希望给其他人一些参考。
一、前台传入geojson图形和影像,以及分类等相关参数
//获取geojson数据
new GeoJSON().writeFeatures(layer.getSource().getFeatures())
{"geoJSON":{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[122.75304112963302,42.499201220100524],[122.8043874485347,42.499201220100524],[122.8043874485347,42.52820855610342],[122.75304112963302,42.52820855610342],[122.75304112963302,42.499201220100524]]]},"properties":null}]}}
二、利用gdal,用geojson数据裁剪影像。
public void imageCutByShp() throws Exception {
gdal.AllRegister()//注册驱动,否则下边执行报错
Dataset source = gdal.Open("E:\\a.tif"); //tif文件路径
Vector vector = new Vector();
vector.add("-cutline");
vector.add("E:\\cut.shp");//shp文件路径
//或者是geojson字符串都可以
//String jsonObject = JSONObject.toJSONString(geoJSON);
// vector.add(jsonObject); //geojson的字符串
vector.add("-crop_to_cutline");
WarpOptions warpAppOptions = new WarpOptions(vector);
Dataset[] datasets = new Dataset[]{source};
Dataset output = gdal.Warp("E:\\output.tif", datasets, warpAppOptions); //第一个参数是生成结果路径和文件名
}
三、调用分类算法,将原始影像tif生成分类后的tif。
(小编是集成用其他小伙伴利用python深度学习完成的工具,对应的exe,实现分类)
四、通过gdal将栅格矢量化输出
try {
String output = "E:\\a.shp";
Dataset dataset = gdal.Open(output);
Band band = dataset.GetRasterBand(1);
SpatialReference prj = new SpatialReference();
if (!dataset.GetProjectionRef().isEmpty()) {
prj.ImportFromWkt(dataset.GetProjectionRef());
}
Driver driver = gdal.GetDriverByName("ESRI Shapefile");
Dataset dsOut = driver.Create(output, dataset.getRasterXSize(),
dataset.getRasterYSize(), 1, gdalconst.GDT_Byte);
Layer layer = dsOut.CreateLayer("layer", prj);
FieldDefn field = new FieldDefn("value", ogr.OFTReal);
//第二个参数是字段类型,ogr.OFTInteger和ogr.OFTReal都是属性字段的数据类型,但它们的区别在于存储的数据类型不同。ogr.OFTInteger表示整数类型,而ogr.OFTReal表示浮点数类型。如果您需要存储整数类型的属性值,可以使用ogr.OFTInteger,如果需要存储浮点数类型的属性值,可以使用ogr.OFTReal。在创建属性字段时,需要根据实际情况选择合适的数据类型。
layer.CreateField(field);
//矢量化
gdal.Polygonize(band, null, layer, 0);
Feature feature = layer.GetNextFeature();
while (null != feature) {
int value = feature.GetFieldAsInteger("value");
if (value == 0) {
layer.DeleteFeature(feature.GetFID());
}
feature = layer.GetNextFeature();
}
layer.SyncToDisk();
File[] files = new File[5];
files[0] = new File(output);
files[1] = new File(output.replace(".shp", ".dbf"));
files[2] = new File(output.replace(".shp", ".prj"));
files[3] = new File(output.replace(".shp", ".shx"));
dataset.delete();
dsOut.delete();
} catch (Exception exception) {
log.error("影像转shp失败:" + exception.toString());
}
以上是tif转shp,小编也试过tif直接转geojson,将驱动改成GeoJSON的,Driver driver = gdal.GetDriverByName("GeoJSON");File[] files = new File[1];files[0] = new File(output);也能转成功,但是layer.GetNextFeature()一直是null,得到的数据里也没有将栅格像素值转到geojson属性里。
五、坐标转换
由于前台openlayer显示的坐标系统是4326的,分析的影像坐标系统是cgcs2000的,前台页面通过openlayer转坐标系统,图形偏的好远,不成功,所以在后台利用geotools中的JTS进行坐标转换
5.1从shp中读取SimpleFeatureCollection
public SimpleFeatureCollection readShp() throws Exception {
SimpleFeatureCollection simpleFeatureCollection = null;
ShapefileDataStore shpDataStore = null;
try {
String input = "E:\\a.shp'";
File file = new File(input);
shpDataStore = new ShapefileDataStore(file.toURL());
//设置字符编码
Charset charset = "GBK"
shpDataStore.setCharset(charset);
String typeName = shpDataStore.getTypeNames()[0];
SimpleFeatureSource featureSource = null;
featureSource = shpDataStore.getFeatureSource(typeName);
simpleFeatureCollection = featureSource.getFeatures();
} catch (Exception exception) {
log.error("读取shp失败:" + exception.toString());
throw new Exception ("读取shp失败:" + exception.toString());
}
return simpleFeatureCollection;
}
5.2坐标转换
/**
* @param simpleFeatureCollection 要素集
* @param currentCRS 当前的坐标系统
* @param targetCRS 转换的坐标系统
*/
public SimpleFeatureCollection project(SimpleFeatureCollection simpleFeatureCollection, CoordinateReferenceSystem currentCRS, CoordinateReferenceSystem targetCRS) throws MismatchedDimensionException, TransformException, FactoryException, TemplateException {
// 当前要素集空间参考
if (null == simpleFeatureCollection || simpleFeatureCollection.size() == 0) {
return simpleFeatureCollection;
}
if (null == targetCRS) {
throw new Exception("坐标转换错误:targetCRS 空了");
}
if(currentCRS==null){
currentCRS = CRS.decode(CRS.lookupIdentifier(simpleFeatureCollection.getSchema().getCoordinateReferenceSystem(), false), true);
}
MathTransform transform = null;
if (!currentCRS.equals(targetCRS)) {
// 判断是否可直接转换
if (CRS.isCompatible(currentCRS, targetCRS, true)) {
// 如果可直接转换 转换完成再叠加
transform = CRS.findMathTransform(currentCRS, targetCRS, true);
} else {
// 否则抛出异常 坐标系不能直接转换。
throw new Exception("坐标转换错误");
}
}
// transform 不为空说明需要转换坐标
SimpleFeatureType resultFeatureType = null;
if (transform != null) {
resultFeatureType = SimpleFeatureTypeBuilder.retype(simpleFeatureCollection.getSchema(), targetCRS);
SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(resultFeatureType);
try (SimpleFeatureIterator iterator1 = simpleFeatureCollection.features()) {
List<SimpleFeature> resultFeatures = new ArrayList<>();
while (iterator1.hasNext()) {
SimpleFeature feature = iterator1.next();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
Geometry newGeometry = JTS.transform(geometry, transform);
featureBuilder.add(newGeometry);
if (feature.getAttributes().size() > 1) {
featureBuilder.addAll(feature.getAttributes().subList(1, feature.getAttributes().size()));
}
SimpleFeature newFeature = featureBuilder.buildFeature(null);
resultFeatures.add(newFeature);
}
simpleFeatureCollection = new ListFeatureCollection(resultFeatureType, resultFeatures);
}
}
return simpleFeatureCollection;
}
使用featureBuilder.add()方式添加,一定要注意一定按照SimpleFeatureType给的字段顺序进行赋值!!!
六、将simpleFeatureCollection 转换成geojson返回
public String featureCollection2GeoJsonStr(SimpleFeatureCollection featureCollection) throws IOException {
FeatureJSON featureJSON = new FeatureJSON(new GeometryJSON(15));// 避免精度丢失
FeatureCollection features = (FeatureCollection) featureCollection;
return featureJSON.toString(features);
}
//显示geojson数据
resultLayer.getSource().addFeatures(new GeoJSON().readFeatures(geojson))
这里需要注意的是由于SimpleFeatureCollection 是 org.geotools.data.simple命名空间下的,featureJSON.toString()方法需要传入的参数FeatureCollection是org.geotools.feature命名空间下的,所以将org.geotools.data.simple.SimpleFeatureCollection 转成org.geotools.feature.FeatureCollection 。