Commit 2159a21e authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

SensorMapGeometryTransformer.compute_areadefinition_sensor2map() now directly...


SensorMapGeometryTransformer.compute_areadefinition_sensor2map() now directly uses pyresample instead of GDAL if the target resolution is given.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 84ab0604
Pipeline #8573 failed with stage
in 53 seconds
......@@ -26,6 +26,7 @@
from typing import Union, List, Tuple, Optional
import os
import warnings
from tempfile import TemporaryDirectory
import numpy as np
......@@ -38,7 +39,7 @@ from py_tools_ds.processing.shell import subcall_with_output
# NOTE: In case of ImportError: dlopen: cannot load any more object with static TLS,
# one could add 'from pykdtree.kdtree import KDTree' here (before pyresample import)
from pyresample.geometry import AreaDefinition, SwathDefinition
from pyresample.geometry import AreaDefinition, SwathDefinition, DynamicAreaDefinition, create_area_def
from pyresample.bilinear import resample_bilinear
from pyresample.kd_tree import resample_nearest, resample_gauss, resample_custom
......@@ -140,102 +141,120 @@ class SensorMapGeometryTransformer(object):
:param tgt_res: target X/Y resolution (e.g., (30, 30))
:return:
"""
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())
return gdal_ds
with TemporaryDirectory() as td:
path_xycoords = os.path.join(td, 'xy_coords.bsq')
path_xycoords_vrt = os.path.join(td, 'xy_coords.vrt')
path_data = os.path.join(td, 'data.bsq')
path_datavrt = os.path.join(td, 'data.vrt')
path_data_out = os.path.join(td, 'data_out.bsq')
# write X/Y coordinate array
if tgt_epsg == 4326:
xy_coords = np.dstack([self.swath_definition.lons,
self.swath_definition.lats])
# xy_coords = np.dstack([self.swath_definition.lons[::10, ::10],
# self.swath_definition.lats[::10, ::10]])
else:
xy_coords = np.dstack(list(transform_coordArray(EPSG2WKT(4326), EPSG2WKT(tgt_epsg),
self.swath_definition.lons,
self.swath_definition.lats)))
write_numpy_to_image(xy_coords, path_xycoords, 'ENVI')
# create VRT for X/Y coordinate array
ds_xy_coords = gdal.Open(path_xycoords)
drv_vrt = gdal.GetDriverByName("VRT")
# noinspection PyUnusedLocal
vrt = raiseErr_if_empty(drv_vrt.CreateCopy(path_xycoords_vrt, ds_xy_coords))
del ds_xy_coords, vrt
# create VRT for one data band
mask_band = np.ones((data.shape[:2]), np.int32)
write_numpy_to_image(mask_band, path_data, 'ENVI')
ds_data = gdal.Open(path_data)
vrt = raiseErr_if_empty(drv_vrt.CreateCopy(path_datavrt, ds_data))
vrt.SetMetadata({"X_DATASET": path_xycoords_vrt,
"Y_DATASET": path_xycoords_vrt,
"X_BAND": "1",
"Y_BAND": "2",
"PIXEL_OFFSET": "0",
"LINE_OFFSET": "0",
"PIXEL_STEP": "1",
"LINE_STEP": "1",
"SRS": EPSG2WKT(tgt_epsg),
}, "GEOLOCATION")
vrt.FlushCache()
del ds_data, vrt
subcall_with_output("gdalwarp '%s' '%s' "
'-geoloc '
'-t_srs EPSG:%d '
'-srcnodata 0 '
'-r near '
'-of ENVI '
'-dstnodata none '
'-et 0 '
'-overwrite '
'-te %s '
'%s' % (path_datavrt, path_data_out, tgt_epsg,
' '.join([str(i) for i in tgt_extent]),
' -tr %s %s' % tgt_res if tgt_res else '',),
v=True)
# get output X/Y size
ds_out = raiseErr_if_empty(gdal.Open(path_data_out))
x_size = ds_out.RasterXSize
y_size = ds_out.RasterYSize
out_gt = ds_out.GetGeoTransform()
# noinspection PyUnusedLocal
ds_out = None
# add 1 px buffer around out_extent to avoid cutting the output image
x_size += 2
y_size += 2
out_gt = list(out_gt)
out_gt[0] -= out_gt[1]
out_gt[3] += abs(out_gt[5])
out_gt = tuple(out_gt)
xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=out_gt, cols=x_size, rows=y_size))
out_extent = xmin, ymin, xmax, ymax
if tgt_res:
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 '
'%s. Provide a target extent where the coordinates are multiples of the pixel sizes to '
'avoid this.' % (str(out_res), str(tgt_res)), UserWarning)
# get area_definition
area_definition = AreaDefinition(area_id='',
description='',
proj_id='',
projection=get_proj4info(proj=tgt_prj),
width=x_size,
height=y_size,
area_extent=list(out_extent),
)
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())
return gdal_ds
with TemporaryDirectory() as td:
path_xycoords = os.path.join(td, 'xy_coords.bsq')
path_xycoords_vrt = os.path.join(td, 'xy_coords.vrt')
path_data = os.path.join(td, 'data.bsq')
path_datavrt = os.path.join(td, 'data.vrt')
path_data_out = os.path.join(td, 'data_out.bsq')
# write X/Y coordinate array
if tgt_epsg == 4326:
xy_coords = np.dstack([self.swath_definition.lons,
self.swath_definition.lats])
# xy_coords = np.dstack([self.swath_definition.lons[::10, ::10],
# self.swath_definition.lats[::10, ::10]])
else:
xy_coords = np.dstack(list(transform_coordArray(EPSG2WKT(4326), EPSG2WKT(tgt_epsg),
self.swath_definition.lons,
self.swath_definition.lats)))
write_numpy_to_image(xy_coords, path_xycoords, 'ENVI')
# create VRT for X/Y coordinate array
ds_xy_coords = gdal.Open(path_xycoords)
drv_vrt = gdal.GetDriverByName("VRT")
# noinspection PyUnusedLocal
vrt = raiseErr_if_empty(drv_vrt.CreateCopy(path_xycoords_vrt, ds_xy_coords))
del ds_xy_coords, vrt
# create VRT for one data band
mask_band = np.ones((data.shape[:2]), np.int32)
write_numpy_to_image(mask_band, path_data, 'ENVI')
ds_data = gdal.Open(path_data)
vrt = raiseErr_if_empty(drv_vrt.CreateCopy(path_datavrt, ds_data))
vrt.SetMetadata({"X_DATASET": path_xycoords_vrt,
"Y_DATASET": path_xycoords_vrt,
"X_BAND": "1",
"Y_BAND": "2",
"PIXEL_OFFSET": "0",
"LINE_OFFSET": "0",
"PIXEL_STEP": "1",
"LINE_STEP": "1",
"SRS": EPSG2WKT(tgt_epsg),
}, "GEOLOCATION")
vrt.FlushCache()
del ds_data, vrt
subcall_with_output("gdalwarp '%s' '%s' "
'-geoloc '
'-t_srs EPSG:%d '
'-srcnodata 0 '
'-r near '
'-of ENVI '
'-dstnodata none '
'-et 0 '
'-overwrite '
'-te %s '
'%s' % (path_datavrt, path_data_out, tgt_epsg,
' '.join([str(i) for i in tgt_extent]),
' -tr %s %s' % tgt_res if tgt_res else '',),
v=True)
# get output X/Y size
ds_out = raiseErr_if_empty(gdal.Open(path_data_out))
x_size = ds_out.RasterXSize
y_size = ds_out.RasterYSize
out_gt = ds_out.GetGeoTransform()
# noinspection PyUnusedLocal
ds_out = None
# add 1 px buffer around out_extent to avoid cutting the output image
x_size += 2
y_size += 2
out_gt = list(out_gt)
out_gt[0] -= out_gt[1]
out_gt[3] += abs(out_gt[5])
out_gt = tuple(out_gt)
xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=out_gt, cols=x_size, rows=y_size))
out_extent = xmin, ymin, xmax, ymax
# get area_definition
area_definition = AreaDefinition(area_id='',
description='',
proj_id='',
projection=get_proj4info(proj=tgt_prj),
width=x_size,
height=y_size,
area_extent=list(out_extent),
)
return area_definition
......
......@@ -26,8 +26,6 @@
from typing import Union, Tuple
import multiprocessing
import sys
import warnings
import numpy as np
from py_tools_ds.geo.projection import WKT2EPSG, proj4_to_WKT
......@@ -97,13 +95,6 @@ class SensorMapGeometryTransformer3D(object):
self.opts['nprocs'] = opts.get('nprocs', multiprocessing.cpu_count())
self.mp_alg = ('bands' if self.lons.shape[2] >= opts['nprocs'] else 'tiles') if mp_alg == 'auto' else mp_alg
# override self.mp_alg if SensorMapGeometryTransformer3D is called by nosetest or unittest
is_called_by_nose_cmd = 'nosetest' in sys.argv[0]
if self.opts['nprocs'] > 1 and self.mp_alg == 'bands' and is_called_by_nose_cmd:
warnings.warn("mp_alg='bands' causes deadlocks if SensorMapGeometryTransformer3D is called within a "
"nosetest console call. Using mp_alg='tiles'.")
self.mp_alg = 'tiles'
@staticmethod
def _to_map_geometry_2D(kwargs_dict: dict) -> Tuple[np.ndarray, tuple, str, int]:
assert [var is not None for var in (_global_shared_lons, _global_shared_lats, _global_shared_data)]
......
......@@ -47,20 +47,20 @@ class Test_SensorMapGeometryTransformer(TestCase):
self.dem_sensor_geo = LoadFile(os.path.join(tests_path, 'data', 'dem_sensor_geo.bsq'))
self.lons = LoadFile(os.path.join(tests_path, 'data', 'lons_full_vnir.bsq'))
self.lats = LoadFile(os.path.join(tests_path, 'data', 'lats_full_vnir.bsq'))
self.dem_area_extent_coarse_subset_utm = [622613.864409047, # LL_x
self.dem_area_extent_coarse_subset_utm = (622613.864409047, # LL_x
5254111.40255343, # LL_x
660473.864409047, # LL_x
5269351.40255343] # UR_y
5269351.40255343) # UR_y
self.expected_dem_area_extent_lonlat = [10.685733901515151, # LL_x
self.expected_dem_area_extent_lonlat = (10.685733901515151, # LL_x
47.44113415492957, # LL_y
11.073066098484848, # UR_x
47.54576584507042] # UR_y
47.54576584507042) # UR_y
self.expected_dem_area_extent_utm = [626938.928052, # LL_x
self.expected_dem_area_extent_utm = (626938.928052, # LL_x
5256253.56579, # LL_y
656188.928052, # UR_x
5267203.56579] # UR_y
5267203.56579) # UR_y
def test_to_sensor_geometry(self):
for rsp_alg in rsp_algs:
......@@ -147,20 +147,20 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
self.lons_3D = np.dstack([lons, lons]) # TODO use different lons per band here
self.lats_3D = np.dstack([lats, lats]) # TODO use different lats per band here
self.dem_area_extent_coarse_subset_utm = [622613.864409047, # LL_x
self.dem_area_extent_coarse_subset_utm = (622613.864409047, # LL_x
5254111.40255343, # LL_x
660473.864409047, # LL_x
5269351.40255343] # UR_y
5269351.40255343) # UR_y
self.expected_dem_area_extent_lonlat = [10.685733901515151, # LL_x
self.expected_dem_area_extent_lonlat = (10.685733901515151, # LL_x
47.44113415492957, # LL_y
11.073066098484848, # UR_x
47.54576584507042] # UR_y
47.54576584507042) # UR_y
self.expected_dem_area_extent_utm = [626938.928052, # LL_x
self.expected_dem_area_extent_utm = (626938.928052, # LL_x
5256253.56579, # LL_y
656188.928052, # UR_x
5267203.56579] # UR_y
5267203.56579) # UR_y
def test_to_map_geometry_lonlat_3D_geolayer(self):
for rsp_alg in rsp_algs:
......
Supports Markdown
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