Coverage for enpt/processors/dem_preprocessing/dem_preprocessing.py: 96%
54 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-03-07 11:39 +0000
« prev ^ index » next coverage.py v7.4.1, created at 2024-03-07 11:39 +0000
1# -*- coding: utf-8 -*-
3# EnPT, EnMAP Processing Tool - A Python package for pre-processing of EnMAP Level-1B data
4#
5# Copyright (C) 2018-2024 Karl Segl (GFZ Potsdam, segl@gfz-potsdam.de), Daniel Scheffler
6# (GFZ Potsdam, danschef@gfz-potsdam.de), Niklas Bohn (GFZ Potsdam, nbohn@gfz-potsdam.de),
7# Stéphane Guillaso (GFZ Potsdam, stephane.guillaso@gfz-potsdam.de)
8#
9# This software was developed within the context of the EnMAP project supported
10# by the DLR Space Administration with funds of the German Federal Ministry of
11# Economic Affairs and Energy (on the basis of a decision by the German Bundestag:
12# 50 EE 1529) and contributions from DLR, GFZ and OHB System AG.
13#
14# This program is free software: you can redistribute it and/or modify it under
15# the terms of the GNU General Public License as published by the Free Software
16# Foundation, either version 3 of the License, or (at your option) any later
17# version. Please note the following exception: `EnPT` depends on tqdm, which
18# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files
19# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore".
20# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE.
21#
22# This program is distributed in the hope that it will be useful, but WITHOUT
23# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
24# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
25# details.
26#
27# You should have received a copy of the GNU Lesser General Public License along
28# with this program. If not, see <https://www.gnu.org/licenses/>.
30"""EnPT pre-processing module for digital elevation models."""
32from typing import Union, Tuple # noqa: F401
33from multiprocessing import cpu_count
34import numpy as np
35from pyproj import CRS
37from geoarray import GeoArray
38from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry, transform_any_prj
39from py_tools_ds.geo.vector.topology import get_footprint_polygon, get_overlap_polygon
41from ..spatial_transform import Geometry_Transformer, get_UTMEPSG_from_LonLat, get_center_coord
43__author__ = 'Daniel Scheffler'
46class DEM_Processor(object):
47 def __init__(self, dem_path_geoarray: Union[str, GeoArray],
48 enmapIm_cornerCoords: Tuple[Tuple[float, float]],
49 CPUs: int = None):
50 self.dem = GeoArray(dem_path_geoarray)
51 self.enmapIm_cornerCoords = enmapIm_cornerCoords
52 self.CPUs = CPUs or cpu_count()
54 self._set_nodata_if_not_provided()
55 self._validate_input()
57 def _validate_input(self):
58 # check geocoding of DEM
59 if not self.dem.is_map_geo:
60 raise ValueError((self.dem.gt, self.dem.prj),
61 'The provided digital elevation model has no valid geo-coding or projection.')
63 # check if provided projection is WGS-84
64 ell = CRS(self.dem.prj).datum.name
65 if not ell.startswith('World Geodetic System 1984'):
66 raise ValueError(ell, "The digital elevation model must be provided with 'WGS84' as geographic datum.")
68 # check overlap
69 dem_ll_mapPoly = reproject_shapelyGeometry(self.dem.footprint_poly, prj_src=self.dem.epsg, prj_tgt=4326)
70 enmapIm_ll_mapPoly = get_footprint_polygon(self.enmapIm_cornerCoords, fix_invalid=True)
71 overlap_dict = get_overlap_polygon(dem_ll_mapPoly, enmapIm_ll_mapPoly)
72 overlap_perc = round(overlap_dict['overlap percentage'], 4)
74 if overlap_perc < 100:
75 # compute minimal extent in user provided projection
76 cornersXY = np.array([transform_any_prj(4326, self.dem.epsg, x, y) for x, y in self.enmapIm_cornerCoords])
77 xmin, xmax = cornersXY[:, 0].min(), cornersXY[:, 0].max()
78 ymin, ymax = cornersXY[:, 1].min(), cornersXY[:, 1].max()
80 raise ValueError('The provided digital elevation model does not cover the EnMAP image completely '
81 '(only around %.1f percent). The minimal needed extent in the provided projection is: \n'
82 'xmin: %s, xmax: %s, ymin: %s, ymax: %s.' % (overlap_perc, xmin, xmax, ymin, ymax))
84 def _set_nodata_if_not_provided(self):
85 # noinspection PyProtectedMember
86 if self.dem._nodata is None:
87 self.dem.nodata = -9999
89 @classmethod
90 def get_flat_dem_from_average_elevation(cls, corner_coords_lonlat, average_elevation, xres=30, yres=30):
91 """Return a 'flat DEM' instance of DEM_Processor.
93 (a GeoArray fully covering the given coorner coordinates with the average elevation as pixel values)
95 :param corner_coords_lonlat: corner coordinates to be covered by the output DEM
96 :param average_elevation: average elevation in meters
97 :param xres: x-resolution in meters
98 :param yres: y-resolution in meters
99 :return:
100 """
101 # compute the dimensions of the flat output DEM
102 tgt_utm_epsg = get_UTMEPSG_from_LonLat(*get_center_coord(corner_coords_lonlat))
103 corner_coords_utm = [transform_any_prj(prj_src=4326, prj_tgt=tgt_utm_epsg, x=x, y=y)
104 for x, y in corner_coords_lonlat]
105 x_all, y_all = list(zip(*corner_coords_utm))
106 cols = int(np.ceil((max(x_all) - min(x_all)) / xres)) + 2
107 rows = int(np.ceil((max(y_all) - min(y_all)) / yres)) + 2
109 # get a GeoArray instance
110 dem_gA = GeoArray(np.full((rows, cols), fill_value=average_elevation),
111 geotransform=(np.floor(min(x_all)) - xres, xres, 0, np.ceil(max(y_all)) + yres, 0, -yres),
112 projection=CRS(tgt_utm_epsg).to_wkt())
114 return cls(dem_gA, corner_coords_lonlat)
116 def fill_gaps(self):
117 pass
119 def compute_slopes(self):
120 # compute on map geometry (as provided)
121 pass
123 def compute_aspect(self):
124 # compute on map geometry (as provided)
125 pass
127 def to_sensor_geometry(self,
128 lons: np.ndarray,
129 lats: np.ndarray):
130 GT = Geometry_Transformer(lons=lons, lats=lats, backend='gdal', resamp_alg='bilinear', nprocs=self.CPUs)
131 data_sensorgeo = GT.to_sensor_geometry(self.dem)
133 return GeoArray(data_sensorgeo)