java20 通过gdal3.6.4，geotools实现遥感影像分类，坐标转换，输出_java后台4326和32649坐标转换

前台是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 。