From 09df3b428eb17cdd5102374a8e8bfb42c780b350 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lvaro?= Date: Fri, 15 Jul 2022 23:13:32 +0200 Subject: [PATCH 1/6] z3nth10n/feature/water_layers: now we can interpolate between water and land rasters --- ddd/geo/elevation.py | 12 ++++++------ ddd/geo/georaster.py | 31 ++++++++++++++++++++++++------- examples/geoterrain.py | 3 ++- examples/p.log | 25 +++++++++++++++++++++++++ examples/periodictable.py | 4 ++-- 5 files changed, 59 insertions(+), 16 deletions(-) create mode 100644 examples/p.log diff --git a/ddd/geo/elevation.py b/ddd/geo/elevation.py index 3a5fc9e..99ed2a6 100644 --- a/ddd/geo/elevation.py +++ b/ddd/geo/elevation.py @@ -43,13 +43,13 @@ def elevation(self, point): logger.warn("Non finite elevation value found at %s: %s", point, value) value = 0 # - 0.01 - if value < -1000.0 or value > 10000.0: - #raise DDDException("Suspicious value for elevation: %s. Aborting." % value) - # (Sea values in EUDEM11 are found to be -3.573423841207179e+38) - value = 0 + # if value < -1000.0 or value > 10000.0: + # #raise DDDException("Suspicious value for elevation: %s. Aborting." % value) + # # (Sea values in EUDEM11 are found to be -3.573423841207179e+38) + # value = 0 - if value is None: - value = 0 + # if value is None: + # value = 0 return value diff --git a/ddd/geo/georaster.py b/ddd/geo/georaster.py index 71619c7..e6951bc 100644 --- a/ddd/geo/georaster.py +++ b/ddd/geo/georaster.py @@ -15,6 +15,8 @@ from ddd.core import settings from ddd.core.exception import DDDException import numpy as np +import sys +import traceback # Get instance of logger for this module @@ -122,12 +124,12 @@ def value_interpolated(self, point, k=3): return self.value_simple(point) # Correct EUDEM11 <10000 values - try: - height_matrix = np.maximum(height_matrix, 0) - except TypeError as e: - # Fix None values (empty / missing) - # TODO: Except on sea, this should possibly get the average from surrounding values - return 0 + # try: + # height_matrix = np.maximum(height_matrix, 0) + # except TypeError as e: + # # Fix None values (empty / missing) + # # TODO: Except on sea, this should possibly get the average from surrounding values + # return 0 if k == 1: interpolated = interp2d([0, 1], [0, 1], height_matrix, 'linear') # , copy=False @@ -150,7 +152,10 @@ def __init__(self, tiles_config): self.tiles_config = tiles_config self._last_tile = None + self._last_landmask_tile = None self._last_tile_config = None + self._last_landmask_config = None + self._last_waterlayer_tile = None # Cache of transformers for different georaster tile CRS self._transformers = {} @@ -191,6 +196,14 @@ def tile_from_point(self, point): # point[1] >= cc['bounds_wgs84_xy'][1] and point[1] < cc['bounds_wgs84_xy'][3]): if (projected_point[0] >= cc['bounds'][0] and projected_point[0] < cc['bounds'][2] and projected_point[1] >= cc['bounds'][1] and projected_point[1] < cc['bounds'][3]): + # TODO: the config must be ordered + if 'is_water_mask' in cc and cc['is_water_mask'] == True: + self._last_landmask_config = cc + self._last_landmask_tile = GeoRasterTile.load(cc['path'], cc['crs']) + continue + if 'is_water_layer' in cc and cc['is_water_layer'] == True: + self._last_waterlayer_tile = GeoRasterTile.load(cc['path'], cc['crs']) + continue self._last_tile_config = cc self._last_tile = GeoRasterTile.load(cc['path'], cc['crs']) return self._last_tile @@ -201,7 +214,11 @@ def value(self, point, interpolate=True): tile = self.tile_from_point(point) if tile is None: raise DDDException("No raster tile found for point: %s" % (point, )) - return tile.value(point, interpolate) + landmask_value=self._last_landmask_tile.value(point, interpolate) + if self._last_landmask_config['water_value'] == landmask_value: + return self._last_waterlayer_tile.value(point, interpolate) + else: + return tile.value(point, interpolate) def area(self, bounds): """ diff --git a/examples/geoterrain.py b/examples/geoterrain.py index cd0cf73..aa5911d 100644 --- a/examples/geoterrain.py +++ b/examples/geoterrain.py @@ -33,7 +33,8 @@ def pipeline_start(pipeline, root): item = item.material(ddd.mats.terrain) item = item.smooth(angle=math.pi/12) item = ddd.uv.map_cubic(item) - item.show() + # item.show() + item.save('/tmp/geo.glb') root.append(item) diff --git a/examples/p.log b/examples/p.log new file mode 100644 index 0000000..cbe966e --- /dev/null +++ b/examples/p.log @@ -0,0 +1,25 @@ +Loading config from: /root/.ddd.conf +2022-04-06 17:35:26,457 DDD logging initialized. +2022-04-06 17:35:26,457 Running ddd.core.commands.run.RunCommand +2022-04-06 17:35:26,457 Properties: {} +/root/earthplusplus/ddd/venv/lib/python3.9/site-packages/pkg_resources/__init__.py:123: PkgResourcesDeprecationWarning: 0.6.5-devel is an invalid version and will not be supported in a future release + warnings.warn( +2022-04-06 17:35:27,693 Loaded 66 materials. +2022-04-06 17:35:27,874 Running DDD123 script. +2022-04-06 17:35:27,874 Appending to path: +2022-04-06 17:35:27,875 Running pipeline: Pipeline(name='DDD Server Build Pipeline') (14 configured tasks) +2022-04-06 17:35:27,875 Running DDDTask(14-cache) (obj=None) +2022-04-06 17:35:27,875 Continuing from cached state: /tmp/ddd.cache +2022-04-06 17:35:27,875 Loading DDD Pipeline cache data from: /tmp/ddd.cache +2022-04-06 17:35:27,935 Running DDDTask(15-type_metal) (obj=7) +2022-04-06 17:35:27,938 Running DDDTask(16-type_noble_gas) (obj=7) +2022-04-06 17:35:27,940 Running DDDTask(17-type_halogen) (obj=4) +2022-04-06 17:35:27,942 Running DDDTask(18-create_base) (obj=90) +2022-04-06 17:35:30,162 Running DDDTask(19-fonts) (obj=None) +2022-04-06 17:35:30,163 Running DDDTask(20-font) (obj=90) +2022-04-06 17:35:30,722 Running DDDTask(21-texts_front) (obj=90) +2022-04-06 17:35:32,229 Running DDDTask(22-position) (obj=90) +2022-04-06 17:35:32,266 Running DDDTask(23-show) (obj=None) +2022-04-06 17:35:32,266 Saving to: /tmp/periodictable.glb (DDDObject3(Elements3_1163, faces=0, children=90)) +2022-04-06 17:35:36,392 This pipeline uses stage caching! use --cache-clear to force a full rebuild. +2022-04-06 17:35:36,392 Pipeline processing time: 0:08.5 m diff --git a/examples/periodictable.py b/examples/periodictable.py index b5627ba..2846359 100644 --- a/examples/periodictable.py +++ b/examples/periodictable.py @@ -251,8 +251,8 @@ def position(root, obj): def show(root, logger): """ """ - #root.find("/Elements3").save("/tmp/periodictable.glb") - root.find("/Elements3").show() + root.find("/Elements3").save("/tmp/periodictable.glb") + #root.find("/Elements3").show() logger.warn("This pipeline uses stage caching! use --cache-clear to force a full rebuild.") From 331b46aa64272e3f34ea48d8d356218405afd24d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lvaro?= Date: Fri, 15 Jul 2022 23:14:11 +0200 Subject: [PATCH 2/6] remove p.log --- examples/p.log | 25 ------------------------- 1 file changed, 25 deletions(-) delete mode 100644 examples/p.log diff --git a/examples/p.log b/examples/p.log deleted file mode 100644 index cbe966e..0000000 --- a/examples/p.log +++ /dev/null @@ -1,25 +0,0 @@ -Loading config from: /root/.ddd.conf -2022-04-06 17:35:26,457 DDD logging initialized. -2022-04-06 17:35:26,457 Running ddd.core.commands.run.RunCommand -2022-04-06 17:35:26,457 Properties: {} -/root/earthplusplus/ddd/venv/lib/python3.9/site-packages/pkg_resources/__init__.py:123: PkgResourcesDeprecationWarning: 0.6.5-devel is an invalid version and will not be supported in a future release - warnings.warn( -2022-04-06 17:35:27,693 Loaded 66 materials. -2022-04-06 17:35:27,874 Running DDD123 script. -2022-04-06 17:35:27,874 Appending to path: -2022-04-06 17:35:27,875 Running pipeline: Pipeline(name='DDD Server Build Pipeline') (14 configured tasks) -2022-04-06 17:35:27,875 Running DDDTask(14-cache) (obj=None) -2022-04-06 17:35:27,875 Continuing from cached state: /tmp/ddd.cache -2022-04-06 17:35:27,875 Loading DDD Pipeline cache data from: /tmp/ddd.cache -2022-04-06 17:35:27,935 Running DDDTask(15-type_metal) (obj=7) -2022-04-06 17:35:27,938 Running DDDTask(16-type_noble_gas) (obj=7) -2022-04-06 17:35:27,940 Running DDDTask(17-type_halogen) (obj=4) -2022-04-06 17:35:27,942 Running DDDTask(18-create_base) (obj=90) -2022-04-06 17:35:30,162 Running DDDTask(19-fonts) (obj=None) -2022-04-06 17:35:30,163 Running DDDTask(20-font) (obj=90) -2022-04-06 17:35:30,722 Running DDDTask(21-texts_front) (obj=90) -2022-04-06 17:35:32,229 Running DDDTask(22-position) (obj=90) -2022-04-06 17:35:32,266 Running DDDTask(23-show) (obj=None) -2022-04-06 17:35:32,266 Saving to: /tmp/periodictable.glb (DDDObject3(Elements3_1163, faces=0, children=90)) -2022-04-06 17:35:36,392 This pipeline uses stage caching! use --cache-clear to force a full rebuild. -2022-04-06 17:35:36,392 Pipeline processing time: 0:08.5 m From 828f7b7dc6c7a6828d1c701e15c671746fe4cc7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lvaro?= Date: Sat, 16 Jul 2022 00:02:17 +0200 Subject: [PATCH 3/6] z3nth10n/feature/water_layers: bug fixing some conditions --- ddd/geo/georaster.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/ddd/geo/georaster.py b/ddd/geo/georaster.py index e6951bc..f6b9aba 100644 --- a/ddd/geo/georaster.py +++ b/ddd/geo/georaster.py @@ -191,34 +191,42 @@ def tile_from_point(self, point): for cc in self.tiles_config: projected_point = self._transform_wgs84_to(point, cc['crs']) + # print(cc) #if (point[0] >= cc['bounds_wgs84_xy'][0] and point[0] < cc['bounds_wgs84_xy'][2] and # point[1] >= cc['bounds_wgs84_xy'][1] and point[1] < cc['bounds_wgs84_xy'][3]): if (projected_point[0] >= cc['bounds'][0] and projected_point[0] < cc['bounds'][2] and projected_point[1] >= cc['bounds'][1] and projected_point[1] < cc['bounds'][3]): # TODO: the config must be ordered - if 'is_water_mask' in cc and cc['is_water_mask'] == True: + if self._last_landmask_tile is None and 'is_water_mask' in cc and cc['is_water_mask'] == True: self._last_landmask_config = cc self._last_landmask_tile = GeoRasterTile.load(cc['path'], cc['crs']) + # print(self._last_landmask_tile) + # print(f"is_water_mask: {cc}") continue - if 'is_water_layer' in cc and cc['is_water_layer'] == True: + if self._last_waterlayer_tile is None and 'is_water_layer' in cc and cc['is_water_layer'] == True: self._last_waterlayer_tile = GeoRasterTile.load(cc['path'], cc['crs']) + # print(f"is_water_layer: {cc}") + # print(self._last_waterlayer_tile) continue self._last_tile_config = cc self._last_tile = GeoRasterTile.load(cc['path'], cc['crs']) return self._last_tile - return None + return None # self._last_waterlayer_tile # TODO def value(self, point, interpolate=True): tile = self.tile_from_point(point) + + if self._last_landmask_tile is not None: + landmask_value=self._last_landmask_tile.value(point, interpolate) + + if self._last_landmask_config is not None and self._last_landmask_config['water_value'] == landmask_value: + return self._last_waterlayer_tile.value(point, interpolate) + if tile is None: raise DDDException("No raster tile found for point: %s" % (point, )) - landmask_value=self._last_landmask_tile.value(point, interpolate) - if self._last_landmask_config['water_value'] == landmask_value: - return self._last_waterlayer_tile.value(point, interpolate) - else: - return tile.value(point, interpolate) + return tile.value(point, interpolate) def area(self, bounds): """ From b236c1543c3b858431357963988c6fe7107277f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lvaro?= Date: Wed, 20 Jul 2022 23:03:19 +0200 Subject: [PATCH 4/6] z3nth10n/feature/water_layers: more powerful implementation you can assign water/landmasks to two variables without problems, the only dependency is that the mask must cover the entire map --- ddd/geo/georaster.py | 78 ++++++++++++++++++++++++++++---------------- 1 file changed, 50 insertions(+), 28 deletions(-) diff --git a/ddd/geo/georaster.py b/ddd/geo/georaster.py index f6b9aba..6b3fc79 100644 --- a/ddd/geo/georaster.py +++ b/ddd/geo/georaster.py @@ -148,14 +148,18 @@ class GeoRasterLayer: Groups several GeoRasterTile configurations. """ + landmaskMapCache=None + def __init__(self, tiles_config): self.tiles_config = tiles_config self._last_tile = None - self._last_landmask_tile = None self._last_tile_config = None - self._last_landmask_config = None + self._last_landmask_tile = None + self._last_landmask_tile_config = None self._last_waterlayer_tile = None + self._last_waterlayer_tile_config = None + GeoRasterLayer.landmaskMapCache = {} # Cache of transformers for different georaster tile CRS self._transformers = {} @@ -175,57 +179,75 @@ def _transform_wgs84_to(self, point, crs): projected_point = tile_crs_transformer.transform(point[0], point[1]) return projected_point - def tile_from_point(self, point): + def tile_from_point(self, point, interpolate=True): """ Returns the GeoRasterTile for the given point. This method caches the last accessed tile, as it is more likely to be hit next. """ - if self._last_tile: - projected_point = self._transform_wgs84_to(point, self._last_tile.crs) - if (projected_point[0] >= self._last_tile_config['bounds'][0] and projected_point[0] < self._last_tile_config['bounds'][2] and - projected_point[1] >= self._last_tile_config['bounds'][1] and projected_point[1] < self._last_tile_config['bounds'][3]): - return self._last_tile + key = self._transform_wgs84_to(point, 'epsg:3857') + key = (round(key[0]), round(key[1])) - for cc in self.tiles_config: + if self._last_tile or self._last_waterlayer_tile and self._last_landmask_tile: + tile = self._last_tile if self._last_tile else self._last_waterlayer_tile + projected_point = self._transform_wgs84_to(point, tile.crs) + is_water = self._last_waterlayer_tile is not None and key in GeoRasterLayer.landmaskMapCache + if is_water: + landmask_value = GeoRasterLayer.landmaskMapCache[key] + if self._last_landmask_tile_config['water_value'] == landmask_value: + return (self._last_waterlayer_tile, self._last_waterlayer_tile_config) + + if not 'is_water_mask' in self._last_tile_config and (projected_point[0] >= self._last_tile_config['bounds'][0] and projected_point[0] < self._last_tile_config['bounds'][2] and + projected_point[1] >= self._last_tile_config['bounds'][1] and projected_point[1] < self._last_tile_config['bounds'][3]): + return (self._last_tile, self._last_tile_config) + + for cc in self.tiles_config: projected_point = self._transform_wgs84_to(point, cc['crs']) - # print(cc) #if (point[0] >= cc['bounds_wgs84_xy'][0] and point[0] < cc['bounds_wgs84_xy'][2] and # point[1] >= cc['bounds_wgs84_xy'][1] and point[1] < cc['bounds_wgs84_xy'][3]): if (projected_point[0] >= cc['bounds'][0] and projected_point[0] < cc['bounds'][2] and projected_point[1] >= cc['bounds'][1] and projected_point[1] < cc['bounds'][3]): # TODO: the config must be ordered - if self._last_landmask_tile is None and 'is_water_mask' in cc and cc['is_water_mask'] == True: - self._last_landmask_config = cc + is_mask_related=False + if self._last_landmask_tile is None and 'is_water_mask' in cc: + self._last_landmask_tile_config = cc self._last_landmask_tile = GeoRasterTile.load(cc['path'], cc['crs']) - # print(self._last_landmask_tile) - # print(f"is_water_mask: {cc}") - continue - if self._last_waterlayer_tile is None and 'is_water_layer' in cc and cc['is_water_layer'] == True: + is_mask_related=True + + if self._last_waterlayer_tile is None and 'is_water_layer' in cc: + self._last_waterlayer_tile_config = cc self._last_waterlayer_tile = GeoRasterTile.load(cc['path'], cc['crs']) - # print(f"is_water_layer: {cc}") - # print(self._last_waterlayer_tile) + is_mask_related=True + + if self._last_landmask_tile is not None and self._last_waterlayer_tile is not None and 'is_water_mask' in cc and cc['is_water_mask'] == True: + landmask_value = self._last_landmask_tile.value(point, interpolate) + GeoRasterLayer.landmaskMapCache[key] = landmask_value + + if self._last_landmask_tile_config['water_value'] == landmask_value: + return (self._last_waterlayer_tile, self._last_waterlayer_tile_config) + + # if configs for watermask/landmask is already set, the last condition from above was already not met and we are still on those configs we must skip + if self._last_landmask_tile is not None and 'is_water_mask' in cc or self._last_waterlayer_tile is not None and 'is_water_layer' in cc: continue + + if is_mask_related: + continue + self._last_tile_config = cc self._last_tile = GeoRasterTile.load(cc['path'], cc['crs']) - return self._last_tile + return (self._last_tile, self._last_tile_config) - return None # self._last_waterlayer_tile # TODO + return (None, None) # self._last_waterlayer_tile # TODO def value(self, point, interpolate=True): - tile = self.tile_from_point(point) - - if self._last_landmask_tile is not None: - landmask_value=self._last_landmask_tile.value(point, interpolate) - - if self._last_landmask_config is not None and self._last_landmask_config['water_value'] == landmask_value: - return self._last_waterlayer_tile.value(point, interpolate) + tile, config = self.tile_from_point(point, interpolate) + path=config['path'] if config else 'None' if tile is None: - raise DDDException("No raster tile found for point: %s" % (point, )) + raise DDDException(f"[{path}] No raster tile found for point: {point}") return tile.value(point, interpolate) def area(self, bounds): From bdfef7dca85002596b75408eb1d71a792c511359 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lvaro?= Date: Wed, 20 Jul 2022 23:23:11 +0200 Subject: [PATCH 5/6] z3nth10n/feature/water_layers: solve last minor issues --- ddd/geo/georaster.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/ddd/geo/georaster.py b/ddd/geo/georaster.py index 6b3fc79..0b29149 100644 --- a/ddd/geo/georaster.py +++ b/ddd/geo/georaster.py @@ -197,11 +197,12 @@ def tile_from_point(self, point, interpolate=True): if is_water: landmask_value = GeoRasterLayer.landmaskMapCache[key] if self._last_landmask_tile_config['water_value'] == landmask_value: - return (self._last_waterlayer_tile, self._last_waterlayer_tile_config) + return self._last_waterlayer_tile - if not 'is_water_mask' in self._last_tile_config and (projected_point[0] >= self._last_tile_config['bounds'][0] and projected_point[0] < self._last_tile_config['bounds'][2] and + if self._last_tile_config and not 'is_water_mask' in self._last_tile_config and (projected_point[0] >= self._last_tile_config['bounds'][0] and projected_point[0] < self._last_tile_config['bounds'][2] and projected_point[1] >= self._last_tile_config['bounds'][1] and projected_point[1] < self._last_tile_config['bounds'][3]): - return (self._last_tile, self._last_tile_config) + return self._last_tile + # return (self._last_tile, self._last_tile_config) for cc in self.tiles_config: projected_point = self._transform_wgs84_to(point, cc['crs']) @@ -227,7 +228,8 @@ def tile_from_point(self, point, interpolate=True): GeoRasterLayer.landmaskMapCache[key] = landmask_value if self._last_landmask_tile_config['water_value'] == landmask_value: - return (self._last_waterlayer_tile, self._last_waterlayer_tile_config) + return self._last_waterlayer_tile + # return (self._last_waterlayer_tile, self._last_waterlayer_tile_config) # if configs for watermask/landmask is already set, the last condition from above was already not met and we are still on those configs we must skip if self._last_landmask_tile is not None and 'is_water_mask' in cc or self._last_waterlayer_tile is not None and 'is_water_layer' in cc: @@ -238,16 +240,17 @@ def tile_from_point(self, point, interpolate=True): self._last_tile_config = cc self._last_tile = GeoRasterTile.load(cc['path'], cc['crs']) - return (self._last_tile, self._last_tile_config) + return self._last_tile + # return (self._last_tile, self._last_tile_config) - return (None, None) # self._last_waterlayer_tile # TODO + return None def value(self, point, interpolate=True): - tile, config = self.tile_from_point(point, interpolate) - path=config['path'] if config else 'None' + tile = self.tile_from_point(point, interpolate) if tile is None: - raise DDDException(f"[{path}] No raster tile found for point: {point}") + # {tile.layer.GetDescription()} + raise DDDException(f"No raster tile found for point: {point}") return tile.value(point, interpolate) def area(self, bounds): From 9e1f261c01cec5e6bd6ffb25d681b6f196186e07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lvaro?= Date: Tue, 26 Jul 2022 23:15:25 +0200 Subject: [PATCH 6/6] z3nth10n/feature/water_layers: interpolating between layer as smooth as possible --- ddd/geo/georaster.py | 45 ++++++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/ddd/geo/georaster.py b/ddd/geo/georaster.py index 0b29149..0854f22 100644 --- a/ddd/geo/georaster.py +++ b/ddd/geo/georaster.py @@ -159,7 +159,6 @@ def __init__(self, tiles_config): self._last_landmask_tile_config = None self._last_waterlayer_tile = None self._last_waterlayer_tile_config = None - GeoRasterLayer.landmaskMapCache = {} # Cache of transformers for different georaster tile CRS self._transformers = {} @@ -179,7 +178,7 @@ def _transform_wgs84_to(self, point, crs): projected_point = tile_crs_transformer.transform(point[0], point[1]) return projected_point - def tile_from_point(self, point, interpolate=True): + def tile_from_point(self, point): """ Returns the GeoRasterTile for the given point. @@ -193,16 +192,14 @@ def tile_from_point(self, point, interpolate=True): tile = self._last_tile if self._last_tile else self._last_waterlayer_tile projected_point = self._transform_wgs84_to(point, tile.crs) - is_water = self._last_waterlayer_tile is not None and key in GeoRasterLayer.landmaskMapCache - if is_water: - landmask_value = GeoRasterLayer.landmaskMapCache[key] - if self._last_landmask_tile_config['water_value'] == landmask_value: + if self._last_waterlayer_tile is not None: + landmask_value = self._last_landmask_tile.value(point, False) + if self._last_landmask_tile_config['water_value'] >= landmask_value: return self._last_waterlayer_tile if self._last_tile_config and not 'is_water_mask' in self._last_tile_config and (projected_point[0] >= self._last_tile_config['bounds'][0] and projected_point[0] < self._last_tile_config['bounds'][2] and projected_point[1] >= self._last_tile_config['bounds'][1] and projected_point[1] < self._last_tile_config['bounds'][3]): return self._last_tile - # return (self._last_tile, self._last_tile_config) for cc in self.tiles_config: projected_point = self._transform_wgs84_to(point, cc['crs']) @@ -223,11 +220,10 @@ def tile_from_point(self, point, interpolate=True): self._last_waterlayer_tile = GeoRasterTile.load(cc['path'], cc['crs']) is_mask_related=True - if self._last_landmask_tile is not None and self._last_waterlayer_tile is not None and 'is_water_mask' in cc and cc['is_water_mask'] == True: - landmask_value = self._last_landmask_tile.value(point, interpolate) - GeoRasterLayer.landmaskMapCache[key] = landmask_value - - if self._last_landmask_tile_config['water_value'] == landmask_value: + if self._last_landmask_tile is not None and self._last_waterlayer_tile is not None: # and 'is_water_mask' in cc and cc['is_water_mask'] == True: + landmask_value = self._last_landmask_tile.value(point, False) + + if self._last_landmask_tile_config['water_value'] >= landmask_value: return self._last_waterlayer_tile # return (self._last_waterlayer_tile, self._last_waterlayer_tile_config) @@ -246,12 +242,33 @@ def tile_from_point(self, point, interpolate=True): return None def value(self, point, interpolate=True): - tile = self.tile_from_point(point, interpolate) + tile = self.tile_from_point(point) if tile is None: # {tile.layer.GetDescription()} raise DDDException(f"No raster tile found for point: {point}") - return tile.value(point, interpolate) + # else: + # # key = self._transform_wgs84_to(point, 'epsg:3857') + # # key = (round(key[0]), round(key[1])) + # maskvalue=self._last_landmask_tile.value(point, False) + # print(f"[{maskvalue}] {tile.layer.GetDescription()} -> {point} || {tile.value(point, interpolate)} {tile.value(point, False)}") # , {self._last_waterlayer_tile.layer.GetDescription()}") + # # print(f"[{maskvalue}, {self._last_landmask_tile.layer.GetDescription()}] {tile.layer.GetDescription()} -> {point}") + + non_interpolated=tile.value(point, False) + interpolated=tile.value(point, interpolate) + + if abs(non_interpolated - 32768) <= 10 or abs(non_interpolated + 32768) <= 10: + if self._last_waterlayer_tile is not None: + return self._last_waterlayer_tile.value(point, interpolate) + return 0 + + if non_interpolated == 0 and abs(interpolated) >= 100: + return 0 + + if abs(non_interpolated) > 0 and abs(interpolated / non_interpolated) >= 3: + return non_interpolated + + return interpolated def area(self, bounds): """