Commit 9176d225 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

GeoArray.projection is now always set to a WKT1 string (GDAL conform) no...


GeoArray.projection is now always set to a WKT1 string (GDAL conform) no matter if the input WKT has an EPSG code or not.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 3278114b
......@@ -40,7 +40,7 @@ from py_tools_ds.convenience.object_oriented import alias_property
from py_tools_ds.geo.coord_calc import get_corner_coordinates
from py_tools_ds.geo.coord_grid import snap_bounds_to_pixGrid
from py_tools_ds.geo.coord_trafo import mapXY2imXY, imXY2mapXY, transform_any_prj, reproject_shapelyGeometry
from py_tools_ds.geo.projection import prj_equal, WKT2EPSG, EPSG2WKT, isLocal
from py_tools_ds.geo.projection import prj_equal, WKT2EPSG, EPSG2WKT, isLocal, CRS
from py_tools_ds.geo.raster.conversion import raster2polygon
from py_tools_ds.geo.vector.topology \
import get_footprint_polygon, polyVertices_outside_poly, fill_holes_within_poly
......@@ -791,10 +791,9 @@ class GeoArray(object):
# for some reason GDAL reads arbitrary geotransforms as (0, 1, 0, 0, 0, 1) instead of (0, 1, 0, 0, 0, -1)
self._geotransform[5] = -abs(self._geotransform[5]) # => force ygsd to be negative
# temp conversion to EPSG needed because GDAL seems to modify WKT string when writing file to disk
# (e.g. using gdal_merge) -> conversion to EPSG and back undos that
# consequently use WKT1 strings here as GDAL always exports transformation results as WKT1
wkt = ds.GetProjection()
self._projection = EPSG2WKT(WKT2EPSG(wkt)) if not isLocal(wkt) else ''
self._projection = CRS(wkt).to_wkt(version="WKT1_GDAL") if not isLocal(wkt) else ''
if 'nodata' not in self._initParams or self._initParams['nodata'] is None:
band = ds.GetRasterBand(1)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment