Commit 7d171c08 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added first working version of SensorMapGeometryTransformer.

parent 9bf1431b
......@@ -3,7 +3,10 @@
import numpy as np
import warnings
import multiprocessing
import os
from tempfile import TemporaryDirectory
# custom
try:
from osgeo import gdal
from osgeo import gdalnumeric
......@@ -11,17 +14,20 @@ except ImportError:
import gdal
import gdalnumeric
# custom
import rasterio
from rasterio.warp import reproject as rio_reproject
from rasterio.warp import calculate_default_transform as rio_calc_transform
from rasterio.warp import Resampling
from pyresample.geometry import AreaDefinition, SwathDefinition
from pyresample.kd_tree import resample_nearest, resample_gauss
from pyresample.utils import get_area_def
from ...dtypes.conversion import dTypeDic_NumPy2GDAL
from ..projection import WKT2EPSG, isProjectedOrGeographic, prj_equal
from ..coord_trafo import pixelToLatLon
from ..coord_trafo import pixelToLatLon, get_proj4info, proj4_to_dict
from ..coord_calc import corner_coord_to_minmax, get_corner_coordinates
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...io.raster.writer import write_numpy_to_image
from ...processing.progress_mon import ProgressBar
from ...compatibility.gdal import get_gdal_func
......@@ -485,3 +491,102 @@ def warp_ndarray(ndarray, in_gt, in_prj=None, out_prj=None, out_dtype=None,
warnings.warn('The output bounds of warp_ndarray do not match the requested bounds!')
return res_arr, res_gt, res_prj
class SensorMapGeometryTransformer(object):
def __init__(self, data, lons, lats, resamp_alg='nearest', **opts):
# type: (np.ndarray, np.ndarray, np.ndarray, str, dict) -> None
self.data = data
self.resamp_alg = resamp_alg
self.opts = dict(radius_of_influence=30,
sigmas=15)
self.opts.update(opts)
self.area_definition = None
self.swath_definition = SwathDefinition(lons=lons, lats=lats)
self.area_extent = [np.min(lons), np.min(lats), np.max(lons), np.max(lats)]
def compute_output_shape(self):
with TemporaryDirectory() as td:
path_lons_lats = os.path.join(td, 'lons_lats.bsq')
path_lons_lats_vrt = os.path.join(td, 'lons_lats.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 lons and lats
# lons_lats = np.dstack([self.swath_definition.lons[::10, ::10], self.swath_definition.lats[::10, ::10]])
lons_lats = np.dstack([self.swath_definition.lons, self.swath_definition.lats])
write_numpy_to_image(lons_lats, path_lons_lats, 'ENVI')
# create VRT for lons and lats
ds_lons_lats = gdal.Open(path_lons_lats)
drv_vrt = gdal.GetDriverByName("VRT")
vrt = drv_vrt.CreateCopy(path_lons_lats_vrt, ds_lons_lats)
del ds_lons_lats, vrt
# create VRT for one data band
mask_band = np.ones((self.data.shape[:2]), np.int32)
write_numpy_to_image(mask_band, path_data, 'ENVI')
ds_data = gdal.Open(path_data)
vrt = drv_vrt.CreateCopy(path_datavrt, ds_data)
vrt.SetMetadata({"X_DATASET": path_lons_lats_vrt,
"Y_DATASET": path_lons_lats_vrt,
"X_BAND": "1",
"Y_BAND": "2",
"PIXEL_OFFSET": "0",
"LINE_OFFSET": "0",
"PIXEL_STEP": "1",
"LINE_STEP": "1",
}, "GEOLOCATION")
vrt.FlushCache()
del ds_data, vrt
os.system('gdalwarp -geoloc -t_srs EPSG:4326 %s %s -srcnodata 0 -r near -of ENVI '
'-dstnodata none -et 0 -overwrite' % (path_datavrt, path_data_out))
# get output X/Y size
ds_out = gdal.Open(path_data_out)
x_size = ds_out.RasterXSize
y_size = ds_out.RasterYSize
del ds_out
return x_size, y_size
def get_area_definition(self, proj4_args, cols_out, rows_out):
"""Get output area definition."""
area_def_out = get_area_def(area_id='',
area_name='',
proj_id='',
proj4_args=proj4_args,
x_size=cols_out,
y_size=rows_out,
area_extent=self.area_extent # xmin, ymin, xmax, ymax
)
return area_def_out
def _resample(self, source_geo_def: object, target_geo_def: object):
if self.resamp_alg == 'nearest':
opts = {k: v for k, v in self.opts.items() if k not in ['sigmas']}
result = resample_nearest(source_geo_def, self.data, target_geo_def, **opts)
else:
opts = {k: v for k, v in self.opts.items()}
result = resample_gauss(source_geo_def, self.data, target_geo_def, **opts)
return result
def to_map_geometry(self, tgt_prj):
cols_out, rows_out = self.compute_output_shape()
self.area_definition = self.get_area_definition(proj4_args=get_proj4info(proj=tgt_prj),
cols_out=cols_out, rows_out=rows_out,
)
return self._resample(self.swath_definition, self.area_definition)
def to_sensor_geometry(self, src_prj, src_extent):
proj4_args = proj4_to_dict(get_proj4info(proj=src_prj))
self.area_definition = AreaDefinition('', '', '', proj4_args, self.data.shape[1], self.data.shape[0],
src_extent)
return self._resample(self.area_definition, self.swath_definition)
......@@ -16,7 +16,7 @@ with open("py_tools_ds/version.py") as version_file:
exec(version_file.read(), version)
requirements = ['gdal', 'numpy', 'shapely', 'six', 'rasterio', 'pandas', 'geopandas',
'scikit-image', 'pyproj', 'spectral']
'scikit-image', 'pyproj', 'spectral', 'pyresample']
setup_requirements = [] # TODO(danschef): put setup requirements (distutils extensions, etc.) here
test_requirements = requirements + ["coverage", "nose", "nose2", "nose-htmloutput", "rednose"]
......
......@@ -24,6 +24,7 @@ dependencies:
- sphinx-argparse
- six
- spectral
- pyresample
- flake8
- pycodestyle<2.4.0 # fixes ImportError: module 'pycodestyle' has no attribute 'break_around_binary_operator'
- pylint
......
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