Commit c54349eb authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'enhancement/add_noepsg_compatibility' into 'master'

Add no-EPSG compatibility.

See merge request !24
parents 3278114b 553e360f
Pipeline #22138 passed with stages
in 19 minutes and 28 seconds
......@@ -2,6 +2,15 @@
History
=======
0.11.0 (2021-04-22)
-------------------
* GeoArray.projection is now always set to a WKT1 string (GDAL conform),
no matter if the input WKT has an EPSG code or not.
* Added compatibility of GeoArray.show_map() and GeoArray.show_footprint() with input WKTs that have no EPSG equivalent.
0.10.12 (2021-04-13)
--------------------
......
......@@ -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
......@@ -85,7 +85,6 @@ class GeoArray(object):
:param progress: show progress bars (default: True)
:param q: quiet mode (default: False)
"""
# TODO implement compatibility to GDAL VRTs
if not (isinstance(path_or_array, (str, np.ndarray, GeoArray)) or
issubclass(getattr(path_or_array, '__class__'), GeoArray)):
......@@ -226,7 +225,6 @@ class GeoArray(object):
@property
def rows(self):
"""Get the number of rows of the associated image array."""
return self.shape[0]
@property
......@@ -240,13 +238,11 @@ class GeoArray(object):
@property
def bands(self):
"""Get the number of bands of the associated image array."""
return self.shape[2] if len(self.shape) > 2 else 1
@property
def dtype(self):
"""Get the numpy data type of the associated image array."""
if self._dtype:
return self._dtype
elif self.is_inmem:
......@@ -258,7 +254,6 @@ class GeoArray(object):
@property
def geotransform(self):
"""Get the GDAL GeoTransform of the associated image, e.g., (283500.0, 5.0, 0.0, 4464500.0, 0.0, -5.0)"""
if self._geotransform:
return self._geotransform
elif not self.is_inmem:
......@@ -283,13 +278,11 @@ class GeoArray(object):
@property
def xgsd(self):
"""Get the X resolution in units of the given or detected projection."""
return self.geotransform[1]
@property
def ygsd(self):
"""Get the Y resolution in units of the given or detected projection."""
return abs(self.geotransform[5])
@property
......@@ -326,8 +319,8 @@ class GeoArray(object):
@property
def epsg(self):
# type: (None) -> int
"""Get the EPSG code of the projection of the GeoArray."""
return WKT2EPSG(self.projection)
@epsg.setter
......@@ -791,10 +784,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)
......@@ -1343,10 +1335,14 @@ class GeoArray(object):
assert self.projection, "A projection is needed for a map visualization. Got '%s'." % self.projection
# get image to plot
# (reproject to LonLat as workaround in case self.epsg is None because cartopy relies on an existing EPSG code)
nodataVal = nodataVal if nodataVal is not None else self.nodata
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, boundsMap=boundsMap,
boundsMapPrj=boundsMapPrj, res_factor=res_factor,
nodataVal=nodataVal)
gA2plot = GeoArray(*self._get_plottable_image(xlim, ylim, band, boundsMap=boundsMap,
boundsMapPrj=boundsMapPrj, res_factor=res_factor,
nodataVal=nodataVal,
out_prj=None if self.epsg is not None else 4326),
nodata=nodataVal)
image2plot = gA2plot[:]
# create map
def get_cartopy_crs_from_epsg(epsg_code):
......@@ -1358,19 +1354,21 @@ class GeoArray(object):
else:
raise NotImplementedError('The show_map() method currently does not support the given projection.')
crs_in = get_cartopy_crs_from_epsg(self.epsg)
crs_out = get_cartopy_crs_from_epsg(out_epsg if out_epsg is not None else self.epsg)
crs_in = get_cartopy_crs_from_epsg(gA2plot.epsg)
crs_out = get_cartopy_crs_from_epsg(out_epsg if out_epsg is not None else gA2plot.epsg)
fig = plt.figure(figsize=figsize)
ax = plt.axes(projection=crs_out)
ax.set_extent(self.box.boundsMap, crs=crs_in)
ax.set_extent(gA2plot.box.boundsMap, crs=crs_in)
palette, vmin, vmax = self._get_cmap_vmin_vmax(cmap, vmin, vmax, pmin, pmax, image2plot, nodataVal)
palette, vmin, vmax = gA2plot._get_cmap_vmin_vmax(cmap, vmin, vmax, pmin, pmax, image2plot, nodataVal)
if nodataVal is not None and np.std(image2plot) != 0: # do not show nodata
image2plot = np.ma.masked_equal(image2plot, nodataVal)
ax.imshow(image2plot, cmap=palette, interpolation=interpolation, vmin=vmin, vmax=vmax,
origin='upper', transform=crs_in, extent=list(self.box.boundsMap))
origin='upper', transform=crs_in,
extent=list(gA2plot.box.boundsMap)
)
# show grid lines
try:
......@@ -1403,7 +1401,7 @@ class GeoArray(object):
import folium
import geojson
lonlatPoly = reproject_shapelyGeometry(self.footprint_poly, self.epsg, 4326)
lonlatPoly = reproject_shapelyGeometry(self.footprint_poly, self.prj, 4326)
m = folium.Map(location=tuple(np.array(lonlatPoly.centroid.coords.xy).flatten())[::-1])
gjs = geojson.Feature(geometry=lonlatPoly, properties={})
......
......@@ -454,6 +454,45 @@ class Test_GeoArray(TestCase):
def test_show_map(self):
self.gA.show_map()
def test_show_map_noepsg(self):
wkt_noepsg = \
"""
PROJCRS["BU MEaSUREs Lambert Azimuthal Equal Area - SA - V01",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["unnamed",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",-15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-60,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1]]]
"""
wkt_noepsg = ' '.join(wkt_noepsg.split())
gA = GeoArray(self.gA[:], self.gA.gt, wkt_noepsg)
gA.show_map()
def test_show_histogram(self):
self.gA.show_histogram()
......
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