Commit 2865b5c9 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added paramter tgt_coordgrid to to_map_geometry methods to automatically move...


Added paramter tgt_coordgrid to to_map_geometry methods to automatically move the output extent to a given coordinate grid.
compute_areadefinition_sensor2map() now also adds 1 pixel around the output extent in the pyresample version just like in the GDAL version.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 2251e5d4
Pipeline #8615 failed with stage
in 1 minute and 6 seconds
......@@ -34,6 +34,7 @@ import gdal
from py_tools_ds.geo.projection import EPSG2WKT, WKT2EPSG, proj4_to_WKT
from py_tools_ds.geo.coord_trafo import get_proj4info, proj4_to_dict, transform_coordArray, transform_any_prj
from py_tools_ds.geo.coord_calc import corner_coord_to_minmax, get_corner_coordinates
from py_tools_ds.geo.coord_grid import find_nearest
from py_tools_ds.io.raster.writer import write_numpy_to_image
from py_tools_ds.processing.shell import subcall_with_output
......@@ -141,17 +142,26 @@ class SensorMapGeometryTransformer(object):
:param tgt_res: target X/Y resolution (e.g., (30, 30))
:return:
"""
# get the target extent if not given
# (this also makes sure that create_area_def does not return a DynamicAreaDefinition)
tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
tgt_extent = tgt_extent or self._get_target_extent(tgt_epsg)
if tgt_res:
# add 1 px buffer around out_extent to avoid cutting the output image
xmin, ymin, xmax, ymax = tgt_extent
tgt_extent = (xmin - tgt_res[0], ymin - tgt_res[1], xmax + tgt_res[0], ymax + tgt_res[1])
# get the area definition
# NOTE: This would return a DynamicAreaDefinition if the extent is not provided
# -> could be transformed to an AreaDefinition by using its .freeze() method and
# passing LonLats and resolution
area_definition = \
create_area_def(area_id='',
projection=get_proj4info(proj=tgt_prj),
area_extent=tgt_extent,
resolution=tgt_res)
if isinstance(area_definition, DynamicAreaDefinition):
area_definition = area_definition.freeze(lonslats=(self.lons, self.lats),
resolution=tgt_res)
out_res = (area_definition.pixel_size_x, area_definition.pixel_size_y)
if tgt_res and tgt_res != out_res:
warnings.warn('With respect to the Lon/Lat arrays the pixel size was set to %s instead of the desired '
......@@ -159,9 +169,6 @@ class SensorMapGeometryTransformer(object):
'avoid this.' % (str(out_res), str(tgt_res)), UserWarning)
else:
tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
tgt_extent = tgt_extent or self._get_target_extent(tgt_epsg)
def raiseErr_if_empty(gdal_ds):
if not gdal_ds:
raise Exception(gdal.GetLastErrorMsg())
......@@ -314,6 +321,7 @@ class SensorMapGeometryTransformer(object):
tgt_prj: Union[str, int] = None,
tgt_extent: Tuple[float, float, float, float] = None,
tgt_res: Tuple = None,
tgt_coordgrid: Tuple[Tuple, Tuple] = None,
area_definition: AreaDefinition = None) -> Tuple[np.ndarray, tuple, str]:
"""Transform the input sensor geometry array into map geometry.
......@@ -321,6 +329,8 @@ class SensorMapGeometryTransformer(object):
:param tgt_prj: target projection (WKT or 'epsg:1234' or <EPSG_int>)
:param tgt_extent: extent coordinates of output map geometry array (LL_x, LL_y, UR_x, UR_y) in the tgt_prj
:param tgt_res: target X/Y resolution (e.g., (30, 30))
:param tgt_coordgrid: target coordinate grid ((x, x), (y, y)):
if given, the output extent is moved to this coordinate grid
:param area_definition: an instance of pyresample.geometry.AreaDefinition;
OVERRIDES tgt_prj, tgt_extent and tgt_res; saves computation time
"""
......@@ -332,6 +342,9 @@ class SensorMapGeometryTransformer(object):
raise ValueError(data.shape, 'Expected a sensor geometry data array with %d rows and %d columns.'
% self.lons.shape[:2])
if tgt_coordgrid:
tgt_res = _get_validated_tgt_res(tgt_coordgrid, tgt_res)
# get area_definition
if area_definition:
self.area_definition = area_definition
......@@ -339,6 +352,12 @@ class SensorMapGeometryTransformer(object):
if not tgt_prj:
raise ValueError(tgt_prj, 'Target projection must be given if area_definition is not given.')
tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
tgt_extent = tgt_extent or self._get_target_extent(tgt_epsg)
if tgt_coordgrid:
tgt_extent = _move_extent_to_coordgrid(tgt_extent, *tgt_coordgrid)
self.area_definition = self.compute_areadefinition_sensor2map(
data, tgt_prj=tgt_prj, tgt_extent=tgt_extent, tgt_res=tgt_res)
......@@ -398,3 +417,25 @@ def _corner_coords_lonlat_to_extent(corner_coords_ll: List, tgt_epsg: int):
tgt_extent = [np.min(x_coords), np.min(y_coords), np.max(x_coords), np.max(y_coords)]
return tgt_extent
def _move_extent_to_coordgrid(extent: Tuple[float, float, float, float],
tgt_xgrid: Tuple[float, float],
tgt_ygrid: Tuple[float, float]):
tgt_xgrid, tgt_ygrid = np.array(tgt_xgrid), np.array(tgt_ygrid)
xmin, ymin, xmax, ymax = extent
tgt_xmin = find_nearest(tgt_xgrid, xmin, roundAlg='off', extrapolate=True)
tgt_xmax = find_nearest(tgt_xgrid, xmax, roundAlg='on', extrapolate=True)
tgt_ymin = find_nearest(tgt_ygrid, ymin, roundAlg='off', extrapolate=True)
tgt_ymax = find_nearest(tgt_ygrid, ymax, roundAlg='on', extrapolate=True)
return tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax
def _get_validated_tgt_res(tgt_coordgrid, tgt_res):
exp_tgt_res = np.ptp(tgt_coordgrid[0]), np.ptp(tgt_coordgrid[1])
if tgt_res and exp_tgt_res != tgt_res:
raise ValueError('The target resolution must be compliant to the target coordinate grid if given.')
return exp_tgt_res
......@@ -31,7 +31,9 @@ import numpy as np
from py_tools_ds.geo.projection import WKT2EPSG, EPSG2WKT, proj4_to_WKT
from py_tools_ds.geo.coord_trafo import get_proj4info, transform_any_prj
from .transformer_2d import SensorMapGeometryTransformer, _corner_coords_lonlat_to_extent
from .transformer_2d import \
SensorMapGeometryTransformer, _corner_coords_lonlat_to_extent, \
_move_extent_to_coordgrid, _get_validated_tgt_res
from pyresample import AreaDefinition
......@@ -110,7 +112,9 @@ class SensorMapGeometryTransformer3D(object):
return data_mapgeo, out_gt, out_prj, kwargs_dict['band_idx']
def _get_common_target_extent(self, tgt_epsg):
def _get_common_target_extent(self,
tgt_epsg: int,
tgt_coordgrid: Tuple[Tuple, Tuple] = None):
if tgt_epsg == 4326:
corner_coords_ll = [[self.lons[0, 0, :].min(), self.lats[0, 0, :].max()], # common UL_xy
[self.lons[0, -1, :].max(), self.lats[0, -1, :].max()], # common UR_xy
......@@ -136,16 +140,20 @@ class SensorMapGeometryTransformer3D(object):
# get the extent
common_tgt_extent = (min(X_prj), min(Y_prj), max(X_prj), max(Y_prj))
if tgt_coordgrid:
common_tgt_extent = _move_extent_to_coordgrid(common_tgt_extent, *tgt_coordgrid)
return common_tgt_extent
def _get_common_area_definition(self,
data: np.ndarray,
tgt_prj: Union[str, int],
tgt_extent: Tuple[float, float, float, float] = None,
tgt_res: Tuple = None) -> AreaDefinition:
tgt_res: Tuple = None,
tgt_coordgrid: Tuple[Tuple, Tuple] = None) -> AreaDefinition:
# get common target extent
tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
tgt_extent = tgt_extent or self._get_common_target_extent(tgt_epsg)
tgt_extent = tgt_extent or self._get_common_target_extent(tgt_epsg, tgt_coordgrid=tgt_coordgrid)
SMGT2D = SensorMapGeometryTransformer(lons=self.lons[:, :, 0], # does NOT affect the computed area definition
lats=self.lats[:, :, 0], # -> only needed for __init__
......@@ -164,6 +172,7 @@ class SensorMapGeometryTransformer3D(object):
tgt_prj: Union[str, int],
tgt_extent: Tuple[float, float, float, float] = None,
tgt_res: Tuple = None,
tgt_coordgrid: Tuple[Tuple, Tuple] = None,
area_definition: AreaDefinition = None) -> Tuple[np.ndarray, tuple, str]:
"""Transform the input sensor geometry array into map geometry.
......@@ -171,8 +180,10 @@ class SensorMapGeometryTransformer3D(object):
:param tgt_prj: target projection (WKT or 'epsg:1234' or <EPSG_int>)
:param tgt_extent: extent coordinates of output map geometry array (LL_x, LL_y, UR_x, UR_y) in the tgt_prj
:param tgt_res: target X/Y resolution (e.g., (30, 30))
:param tgt_coordgrid: target coordinate grid ((x, x), (y, y)):
if given, the output extent is moved to this coordinate grid
:param area_definition: an instance of pyresample.geometry.AreaDefinition;
OVERRIDES tgt_prj, tgt_extent and tgt_res; saves computation time
OVERRIDES tgt_prj, tgt_extent, tgt_res and tgt_coordgrid; saves computation time
"""
if data.ndim != 3:
raise ValueError(data.ndim, "'data' must have 3 dimensions.")
......@@ -181,13 +192,16 @@ class SensorMapGeometryTransformer3D(object):
raise ValueError(data.shape, 'Expected a sensor geometry data array with %d rows, %d columns and %d bands.'
% self.lons.shape)
if tgt_coordgrid:
tgt_res = _get_validated_tgt_res(tgt_coordgrid, tgt_res)
init_opts = self.opts.copy()
if self.mp_alg == 'bands':
del init_opts['nprocs'] # avoid sub-multiprocessing
# get common area_definition
if not area_definition:
area_definition = self._get_common_area_definition(data, tgt_prj, tgt_extent, tgt_res)
area_definition = self._get_common_area_definition(data, tgt_prj, tgt_extent, tgt_res, tgt_coordgrid)
args = [dict(
resamp_alg=self.resamp_alg,
......
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