Commit 9217744d authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added compatibility of GeoArray.show_map() and GeoArray.show_footprint() with...


Added compatibility of GeoArray.show_map() and GeoArray.show_footprint() with input WKTs that have no EPSG equivalent.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 9176d225
......@@ -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
......@@ -1342,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):
......@@ -1357,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:
......@@ -1402,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