From e49d001acc57c9cda43a2b522234ad92c8bd4fcd Mon Sep 17 00:00:00 2001 From: Peter Urban Date: Mon, 29 Jul 2024 10:18:54 +0200 Subject: [PATCH] > 0.23.4 - overview.nav_plot.create_figure function: reproject background image to latlon by default --- meson.build | 2 +- python/tests/overview/test_navplot.py | 39 ++++++++++++ .../pingprocessing/overview/nav_plot.py | 60 ++++++++++++++++-- .../background_maps/test_latlon.tiff | Bin 0 -> 1773 bytes unittest_data/background_maps/test_utm.tiff | Bin 0 -> 1761 bytes 5 files changed, 96 insertions(+), 5 deletions(-) create mode 100644 python/tests/overview/test_navplot.py create mode 100644 unittest_data/background_maps/test_latlon.tiff create mode 100644 unittest_data/background_maps/test_utm.tiff diff --git a/meson.build b/meson.build index 0269855..41d417a 100644 --- a/meson.build +++ b/meson.build @@ -9,7 +9,7 @@ project( 'cpp', license: 'MPL-2.0', - version: '0.8.0', + version: '0.8.1', default_options: ['warning_level=2', 'buildtype=release', 'cpp_std=c++20'], meson_version: '>=1.3.2' #first version with clang-cl openmp support ) diff --git a/python/tests/overview/test_navplot.py b/python/tests/overview/test_navplot.py new file mode 100644 index 0000000..f1efbcc --- /dev/null +++ b/python/tests/overview/test_navplot.py @@ -0,0 +1,39 @@ +import themachinethatgoesping as Ping +import rasterio as rio +import os + +class TestNavPlot: + """ + This class contains unit tests for navigation plot functions. + """ + + def test_create_figure_background_image_projection(self): + src_utm = '../../../unittest_data/background_maps/test_utm.tiff' # map in UTM projection EPSG:32631 + src_latlon = '../../../unittest_data/background_maps/test_latlon.tiff' # map in latlon projection EPSG:4326 + + src_utm = os.path.join(os.path.dirname(__file__), src_utm) + src_latlon = os.path.join(os.path.dirname(__file__), src_latlon) + + # default projection should be EPSG:4326 + _,_,crs = Ping.pingprocessing.overview.nav_plot.create_figure('nav', background_image_path=src_latlon, return_crs=True) + assert crs == rio.crs.CRS.from_epsg(4326) + + # utm should be converted to default projection EPSG:4326 + _,_,crs = Ping.pingprocessing.overview.nav_plot.create_figure('nav', background_image_path=src_utm, return_crs=True) + assert crs == rio.crs.CRS.from_epsg(4326) + + # if dst_crs is set to None, default projecttion should be preserved + _,_,crs = Ping.pingprocessing.overview.nav_plot.create_figure('nav', background_image_path=src_latlon, return_crs=True, dst_crs = None) + assert crs == rio.crs.CRS.from_epsg(4326) + + # if dst_crs is set to None, default projecttion should be preserved + _,_,crs = Ping.pingprocessing.overview.nav_plot.create_figure('nav', background_image_path=src_utm, return_crs=True, dst_crs = None) + assert crs == rio.crs.CRS.from_epsg(32631) + + # if dst_crs is set to another projection, this projection should be used + _,_,crs = Ping.pingprocessing.overview.nav_plot.create_figure('nav', background_image_path=src_latlon, return_crs=True, dst_crs = 'EPSG:32631') + assert crs == rio.crs.CRS.from_epsg(32631) + + # if dst_crs is set to another projection, this projection should be used + _,_,crs = Ping.pingprocessing.overview.nav_plot.create_figure('nav', background_image_path=src_utm, return_crs=True, dst_crs = 'EPSG:32631') + assert crs == rio.crs.CRS.from_epsg(32631) \ No newline at end of file diff --git a/python/themachinethatgoesping/pingprocessing/overview/nav_plot.py b/python/themachinethatgoesping/pingprocessing/overview/nav_plot.py index 320316c..51f301f 100644 --- a/python/themachinethatgoesping/pingprocessing/overview/nav_plot.py +++ b/python/themachinethatgoesping/pingprocessing/overview/nav_plot.py @@ -6,13 +6,50 @@ from matplotlib import pyplot as plt import rasterio.plot as rioplt +import rasterio.warp as riowarp import rasterio as rio +from rasterio.io import MemoryFile +from contextlib import contextmanager from themachinethatgoesping.pingprocessing.core.progress import get_progress_iterator +#Adapted from: https://gis.stackexchange.com/questions/443822/how-do-you-reproject-a-raster-using-rasterio-in-memory +@contextmanager +def reproject_raster(src, dst_crs): + + src_crs = src.crs + transform, width, height = riowarp.calculate_default_transform(src_crs, dst_crs, src.width, src.height, *src.bounds) + kwargs = src.meta.copy() + + kwargs.update({ + 'crs': dst_crs, + 'transform': transform, + 'width': width, + 'height': height}) + + with MemoryFile() as memfile: + with memfile.open(**kwargs) as dst: + for i in range(1, src.count + 1): + riowarp.reproject( + source=rio.band(src, i), + destination=rio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=transform, + dst_crs=dst_crs, + resampling=riowarp.Resampling.nearest) + with memfile.open() as dataset: # Reopen as DatasetReader + yield dataset # Note yield not return as we're a contextmanager + def create_figure( - name: str, aspect: str = "equal", close_plots: bool = True, background_image_path: str = None, **kwargs + name: str, + aspect: str = "equal", + close_plots: bool = True, + background_image_path: str = None, + dst_crs = 'EPSG:4326', + return_crs = False, + **kwargs ) -> Tuple[plt.Figure, plt.Axes]: """ Create a figure with a given name and aspect ratio. @@ -38,8 +75,6 @@ def create_figure( # initialize axis ax.grid(True, linestyle="--", color="gray", alpha=0.5) - ax.set_xlabel("longitude") - ax.set_ylabel("latitude") ax.set_title(name) ax.set_aspect(aspect) @@ -47,8 +82,25 @@ def create_figure( background_map = rio.open(background_image_path) _kwargs = {"cmap": "Greys_r"} _kwargs.update(kwargs) - rioplt.show(background_map, ax=ax, **_kwargs) + + if dst_crs is None or dst_crs == background_map.crs: + rioplt.show(background_map, ax=ax, **_kwargs) + dst_crs = background_map.crs + else: + with reproject_raster(background_map, dst_crs) as reprojected_map: + rioplt.show(reprojected_map, ax=ax, **_kwargs) + dst_crs = reprojected_map.crs + + if dst_crs.is_geographic: + ax.set_xlabel("longitude") + ax.set_ylabel("latitude") + elif dst_crs.is_projected: + ax.set_xlabel("easting") + ax.set_ylabel("northing") + + if return_crs: + return fig, ax, dst_crs return fig, ax diff --git a/unittest_data/background_maps/test_latlon.tiff b/unittest_data/background_maps/test_latlon.tiff new file mode 100644 index 0000000000000000000000000000000000000000..d6d28706d5e87c5c6dd4804b2659c079e3814175 GIT binary patch literal 1773 zcmd5+T}YE*6h7a!=G>Odm}Xc_(hDtuVnIQK>uU3I{prP^uD&ixSs`6mP+DUK)rz1} z=nt0^!P5Q^6Zy`&n+h!>VR`4$erKtbQiRp>Ex-8sZ_!0BobNpEd(QKm^S&RPii(9) zGlR%MBpFGfG$Kid4b!*?G2=8=^UOM&E1A+<98;d|JbxNPoo4HD=IV0RPGci-a`xIu z^><(m`|M0Mj2V%aWwXUOZ*n=k1S`h6xi=3OLIetcr)UumrA zunMC2BGF~|`w7o6-Dmr? ziyTVb99&-l{sL(CDZXTvUb^JVqEFWO BVQl~a literal 0 HcmV?d00001 diff --git a/unittest_data/background_maps/test_utm.tiff b/unittest_data/background_maps/test_utm.tiff new file mode 100644 index 0000000000000000000000000000000000000000..e8fca39524614a8859c6eada62c2caca19e25f93 GIT binary patch literal 1761 zcmd5+T}YEr7(V+pb*-72bN0ua8Mq}{vMwqZ*oEnYqP3)g(N#^aDsKvZAT`IXQb7>3 zZfa5xHd$dwMSky5L2zP_G3p|JN)Wb-kX$B4J>S~o%gofP&N#f!ch37h?{m)g?XlVR zv&sa7Q9x-KeoiBY;9HRk&cU|4dgWj#Y z&me4w!yQURotC>?Tlncd6=EizbF5ECmOGPA)htk!I{|fit_kkqVS6X#`crwkDs4?2 zL2oubiaC;c*6J2X=VLs}t>&YFMXJ>e(q!Wtwe5XJuU@KG<%yUwxW%#^=bj;FOER<% zQ=dT{cKTDnw?`}N8-d^d#>b=zB#F;v%`iv?pb%8oF`&x+E5)FQ09f3_18+p vI{X#B^)G)dfU(B?^*C_cU#|hj{q-quLzlpk;eF%4@$vtAZGPdw<-vXe8yiMv literal 0 HcmV?d00001