Python gdal投影信息

Question

我正在使用 gdal_polygonize python 脚本对光栅图像进行多边形化并将多边形存储在 postgres 数据库中。到目前为止一切正常。

gdal_polygonize.py abia.tif -f PostgreSQL PG:"dbname='abiaDB' host='127.0.0.1' port='5434' user='postgres' password='****'" mylayer

将数据存储在数据库中后，我从那里将其导出为 GeoJson 文件，我想用 Leaflat 显示该文件。因此 Leaflat 仅适用于投影 EPSG:4326。因此我使用了 proj4.js 插件，它从 Gauss-Kruger zone 4 EPSG:31468 投影进行转换。所以我必须像这样在代码中手动定义原始投影：

proj4.defs('EPSG:31468', '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs');

L.Proj.geoJson(data, {
  'pointToLayer': function(feature, latlng) {
    return L.marker(latlng).bindPopup(feature.properties.name);
  }
 }).addTo(map);

有没有一种方法可以让我说 python 脚本也应该将投影信息存储在数据库中。我的目标是可视化更加自动化，因此当有其他图像时，使用其他投影应该从数据库中获取投影信息。有没有办法可以说 gdal_polygonize 函数应该将信息存储在额外的行或类似的东西中？

多边形化 Java

gdal.AllRegister();
ogr.RegisterAll();
args = gdal.GeneralCmdLineProcessor(args);

//Open source file
Dataset hDataset = gdal.Open(args[0], gdalconstConstants.GA_ReadOnly);        
Band rasterBand = hDataset.GetRasterBand(1);
Band maskBand = rasterBand.GetMaskBand();

Driver driver = ogr.GetDriverByName("Memory");
DataSource dataSource =  driver.CreateDataSource("mem");

SpatialReference srs = null;
if(!hDataset.GetProjectionRef().isEmpty())
    srs = new SpatialReference(hDataset.GetProjectionRef());

Layer outputLayer = dataSource.CreateLayer("mylayer", srs);
FieldDefn field_def = new FieldDefn("DN",ogr.OFTInteger);
outputLayer.CreateField(field_def);
gdal.Polygonize(rasterBand, maskBand, outputLayer, 0, new Vector<>(), new TermProgressCallback());   

//Transformation       
DataSource dataDest = driver.CreateDataSource("mem2");
//Create destination projection
SpatialReference dst = new SpatialReference();
dst.ImportFromEPSG(4326);

 CoordinateTransformation ct = CoordinateTransformation.CreateCoordinateTransformation(srs, dst);

//Write data to database
DataSource dataSourceDb = ogr.Open( "PG:dbname='abiaDB' host='127.0.0.1' port='5434' user='postgres' password='****'", 1 );
dataSourceDb.CopyLayer(outputLayer, "mylayer");   

Answer 1

一个想法是您可以修改gdal_polygonize.py 以重新投影结果。使用 MEM 驱动程序存储来自dst_layer 的临时多边形化结果，然后使用 osr 模块重新投影。

另一个更简单的想法是使用postgis创建一个数据库VIEW，它可以投影SRID=4326的表的几何列，即

CREATE VIEW mytable_latlong AS
  SELECT gid, ST_Transform(geom, 4326) AS geom, ...
  FROM mytable