Commit 8160ea9f authored by Daniel Scheffler's avatar Daniel Scheffler

GeoArray.set_gdalDataset_meta(): Bugfix for returning gt with positive ygsd in...

GeoArray.set_gdalDataset_meta(): Bugfix for returning gt with positive ygsd in case of arbitrary coordinates.
parent 7cf18eb8
......@@ -685,6 +685,11 @@ class GeoArray(object):
self._shape = tuple([ds.RasterYSize, ds.RasterXSize] + ([ds.RasterCount] if ds.RasterCount > 1 else []))
self._dtype = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType)
self._geotransform = ds.GetGeoTransform()
# for some reason GDAL reads arbitrary geotransforms as (0, 1, 0, 0, 0, 1) instead of (0, 1, 0, 0, 0, -1)
# => force ygsd to be negative
self._geotransform = tuple(list(self._geotransform)[:5] + [-abs(self._geotransform[5])])
# 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
self._projection = EPSG2WKT(WKT2EPSG(ds.GetProjection()))
......@@ -709,9 +714,15 @@ class GeoArray(object):
# TODO add the remaining global metadata
# avoid double-call of set_gdalDataset_meta by setting self._metadata to default value
self._metadata = \
self._metadata if self._metadata is not None else GeoDataFrame(columns=range(self.bands))
# fill metadata
self.metadata[b] = meta_gs
del band
del ds, band
del ds
self._gdalDataset_meta_already_set = True
......@@ -1086,7 +1097,7 @@ class GeoArray(object):
assert self.geotransform and tuple(self.geotransform) != (0, 1, 0, 0, 0, -1), \
'A valid geotransform is needed for a map visualization. Got %s.' % list(self.geotransform)
assert self.projection, 'A projection is needed for a map visualization. Got %s.' % self.projection
assert self.projection, "A projection is needed for a map visualization. Got '%s'." % self.projection
# get image to plot
nodataVal = nodataVal if nodataVal is not None else self.nodata
......
......@@ -71,7 +71,7 @@ class Test_GeoarrayAppliedOnPathArray(unittest.TestCase):
(R_exp, C_exp, B_exp) = expected_shape = (10, 11, 2)
expected_result = (3, R_exp, C_exp, B_exp) # dimensions, rows, columns, bands
expected_dtype = np.dtype('float32')
given_geotransform = (365985.0, 30.0, 0.0, 5916615.0, 0.0, -30.0)
given_geotransform = [365985.0, 30.0, 0.0, 5916615.0, 0.0, -30.0]
expected_resolution = (30, 30)
expected_grid = [[365985.0, 366015.0], [5916615.0, 5916585.0]]
given_pszProj4_string = '+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs'
......
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