From 06325bf8dc86d4aafd18362463294971d4844ab4 Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 12:46:35 +0100 Subject: [PATCH 1/8] adding test case for two UTM zones --- tests/conftest.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/conftest.py b/tests/conftest.py index 5b5edb3..fb919ce 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -254,6 +254,21 @@ def _get_polygons(): return _get_polygons +@pytest.fixture() +def get_polygons_utm_zones(get_project_root_path): + """ + Returns path to an area of interest at the boundary of two UTM zones + """ + def _get_polygons(): + + testdata_dir = get_project_root_path.joinpath('data') + testdata_polys = testdata_dir.joinpath( + Path('sample_polygons').joinpath('utm_zones_western-ch.gpkg') + ) + return testdata_polys + return _get_polygons + + @pytest.fixture() def get_scene_collection(get_bandstack): """fixture returing a SceneCollection with three scenes""" From 3b60d56cd5d9b0f7e09209c75eed361f6a7cc746 Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 12:46:48 +0100 Subject: [PATCH 2/8] adding test case for two UTM zones --- data/sample_polygons/utm_zones_western-ch.gpkg | Bin 0 -> 106496 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 data/sample_polygons/utm_zones_western-ch.gpkg diff --git a/data/sample_polygons/utm_zones_western-ch.gpkg b/data/sample_polygons/utm_zones_western-ch.gpkg new file mode 100644 index 0000000000000000000000000000000000000000..f7a113d5636c93e0356ea27633db4bd8a9379cef GIT binary patch literal 106496 zcmeI5Uu+vkdcc?ZN6D5%=MtGwd_IqvQz+IeiKPBgmYmYmN}?l?3c0kc_>8>`Iixls zcj?`gB^kx_BqeE^qP?QmhqSmBXc0eLacCZ{&808NV_yzwfdWN{)lNcMJpV<*pY?oh9@vH87zh9X zAOHk_01yBIKmZ5;0U!VbfB+Eqk_eo7mSpVl^pMZC6CX0Ytvz>7eAv78C6xx*0s$ZZ z1b_e#00KY&2mpbHLEw#7|Ddxo!1T5JZa))S3`K&Y7tjb=x*8ehR?-Wdohzy2m00?Q z&l%@q7r5ktE!@Y>MfuEfRnEtzR~ETcGCp_1cNuFME~3!XSjfjtj80CDMka?tfMl<|kG8N}nZ}@1zKGruG437D3-fR*x5gHHDn3Ys~nImN`s2IzeQg+S9 z2AhaaMS>wJo=L>{8@?Pa5OE?o91Km4jgN*S;qg!;GB!HeBxxiR4ma0ftg1uu3K!$! z$;1ucqM}PO7LuFYdK z%E&pKCq~Ml0#*vDSYDH|D4j|tk_m3ZimM*jg|dpG<|J7p&w5py$10Yy7-`#@j!QRw zWo_<@OrPuAdaab0sOnggMA<6jSJlj`2dE}=Kp}h2{CU3X0&&#WrJb|a7Sd>HjH;Kz zx|F3JbL!tMmoy`*9!XAxP-JS%?vj=W@y4mp=p^yPsp0VWWH>xJKGDn%CWE7f$2EHC z#AuiTew4mw*nX+vH$K*~y%}tIhmWnal!ks~%+$DONsd@9h-&mTB*`d;vrKiK7%(HMXJUSE#4TU0nC^$VjF&!Qs9v=(85{gmH#*m3}piWAOHk_ z01yBIKmZ5;0U!VbfB+EqG6;0F_IH^x2QdEsGWbHCKmZ5;0U!VbfB+Bx0zd!=00AHX z1daoN+W6o6{{Lr;=df00fRJf%@nFKW98YKd!ohUO)f{00AHX1b_e#00KY& z2mk>f00e-*VF*0ua`*S0Y3b}VfAUXWyHHw5F0C#k6W02FzhpeWJWNg?0|bBo5C8%| z00;m9AOHk_01yBIKmZ6l5d?bM`ui$t|Ivv6SpV;duo@@{1b_e#00KY&2mk>f00e*l z5C8%|;0Oeq-~VS$9U%)y0tA2n5C8%|00;m9AOHk_01yBIKmZ6lF$C!Ie;EHiG3EkA zfdCKy0zd!=00AHX1b_e#00KY&2pow3eg7Yf|Br+Y5&;1q00e*l5C8%|00;m9AOHk_ z01$X$2*C6IC&pBuC=dVwKmZ5;0U!VbfB+Bx0zd!=0D&VD=;{5G@w9xx^!{)6zxMv9 z_cyv--G9=xcJjTG>+WxLooj!;eY0(<_M=>4?i&ki95B0%6V5ZJxqai4j4 zsPpynx+)b4SS^&-3xX~cu}~68(3WsVQA$D)-^N8@qns0UTrUwzaZ#Q_d@8=Mz@?DC zKE5#rKf27#CsQ2Ctjy7*N*>ipLCHj89)3jR5f_auqEzx4Nxl$IpftC{#dws8#?#zj z^l~!Av&e7ih!h#=;;TiB{0;i}{phx&6iH1qlvQL+S4BzIr|l|aMOjvK^eRS0tZ7JJ z6J-?i2hg=ej%YoY;I0kZ6rMrBz(vF*=5}X#O_=S6fGNouR`q&V_XE?+p$apeNOLLv z0r8o6J|288_H=mDLoP;=b9hU;Rg`orh#R_MoD21&giw7*?_R9i?Hw9&y%{lVUmst8 z-ebXzUS8r*eJmO@G7)64gb~kO=TRa_{xVBTtj+fE1kce{UyZAH&8C#o*U)~kYgj6* z>5b8b-|dssh}$n=S=u80xg}`Q9gOCZ8Cs#LN@nDC^T=2-k>-gDC3sX`7j(UZR#Nfh zXlfO`#I4#Hur>z?u(o7*lB`q`Y_zCv#}^XyIRsFOnu6q2g9E zSrXbP>qJ@5OWH$YS!X6AZ&E$lnbr(s*0<3(T3OU3u_&lGFKF8u8mtHg7EX3}Z;T(> zPYL-(U#0fwew}f5%0YJ ztjudjeFBtjW|spxo;()y4Hc78HsspH zVUC-RW|qt=PaxXS;mwR3h|6AX#qaDv?!B+JxV?VA>y6D?-5m9@PFfr5WJTE<+5r_J zcZ;GxUCAlKN;%T$OL+;aW`x;#*&Nogs#Mk`a`R?faILEQ{cKUxbRkRb!wTg3M(!E; z_;T$Ehb+NCax<6Lr4k+q7G&+My-GN>?j?%kUJJ(2Chj-b3poJ|fR2;!NZI84bKB?L8t)y*awwYr>%v+Zn@XL8dG%;=FSSpPYVQ+SMJ#jW$jVK zaTGO2E>%r0-#MI@WXZ^)t~u7XyFgFp-!omz$4t)$-39mS9fNKE()#U|G3H}WdvCG( z)9&KQPnatw?siSzd%Q+6`u%%LJt~RaRj>O@h79&D?l(TM-zpw;xYs0)>KLUm;7iUU zyVL`YS&lyBb0|^T(2OBpK_&e#>6+?BeYD9LBR@x?OVzPo{ItpT%%RvGar9>|TR2|h zKf7w^0a>p09jcw0F`V?4&UFI^0h0++$j_5Q-{L(U>6P(>+!l1e>fO-ww`r+ z`}~`> zFY?vXYwf00e+Qa{@b;T)oWh zX?H9A*VfX~varImc9QQG{Le?Tu5dAYr^f&Kb`n^je6DzTm`ETF!oNgnT)u)g4 z-OPONAIAEVI2Y~vk$JeyMEl+{50$X%ZzmXj>gn&D{9S2N9sY)*9OXV>0pD2C;(+*Mu0nEZ;5tRHn>1sLfJ&4(1ZxJ}ppZ}t2U*2>e<*nY?g0<`C%0v`ywNT2VWhd~k|7(zTpxCzA)2 zXZ_9Q4$SZWf6jP5e?*moVt@b;00KY&2mk>f00e*l5C8%|00;nqM@gWStU=IajQ>v# zG2}mNKmZ5;0U!VbfB+Bx0zd!=*a`g3imQ$J?$=2WUAc>Ht}9>tyZ^j$_W3vCfB3^# z{JZJ*&0zQh0zd!=9Df3I{Lh>^{(6FjKmZ5;0U!VbfB+Bx0zd!=00AHX1c1QN3DEKX zQ_MFR&tLWab?-*cXFW>ye{_HI^0A z#;X+_KTg=^M-8>$=NvCz+$Sv>to}Aa(6UMyBc8j?Bl-^eL}qD;HOduav1CTDS*5s9 zl5GL1vMES8;kH=Zz}5>LSzRh&O&3dLG#BMLKEBLVb3}98d^EGfqd`qq^HgxqKd?G9 zP#PM@@dJy~1IyC``B(g`Usg8#fdE^|33;hV)@sOExmUBIab=OLNuk?f*;q1>=2Ou) zd3$75EvZ5&FA4c|G@ne7_c$idOWZ0NtYsP?Z;+hlQd}a&rICZtpp*+VW>#*RSyhXJ zvkUB=9CUmAe%D^6(VmXN`}emUB8pW}*eY8R%InDSCe3urfvSmltH$utB|bE)`zcdm zIXNe)IblPU>Y77|*h|p2>zhd}ANu0Yk*!jZWK>x=#F;})(jVB_^1IKZ`#bkMmG`e! zc0vvpFXRbHkN)$XcBnG7VozF`r5e`El zOU&(df8BlNMbds>+P5YfW#U+P?66g-h+B23G)a3~HR@7T1D&5m3GSMy6N)D0(4cWj zDU#RgrfC3z4hMvBcB2?%=HR%0mV5jx0F&x@CyN?EBFe2p|RWO44?xlkv zBb_@#&6Q{!RyUpLRg2i$t~?UE?mladcQ!q&8x#jY)DiFUgWl)O3y!*d`+t`$()n zh9FstGKu)h8P2W*i)_PF7C8$K1XlVxXvJ5JdRM?{_|^tp)UfXKD~@{CMGn&8FiTPf zdOTMdSNB^?_gmxxHXr~5fB+Bx0zd!=00AHX1c1OJA@D}L#m$_VJ^NOy@+7CSta(Yz zW{0y%Y2;Qx(uQ+-&U!#IbAB!v<5yQWL}SoO=JHZJhWz^_^CN$HHO+I&{!1rjj1;%A zN=KNpzR}_E)Nm*~G%t}yad)uq(g}oSBr;UirM!gIOEV+36B=huDcAyMMjC>sTv3z@ z8)P|ol0)^7hSTRyB+S?u=aODD0*pNk5nm+_7saejq^gHTj*`5tjTj2%G1f@5v7=FB z5UAuNxj@3IhcrB^YE+&aX=PED#3FfKs%>jXCMUP?yrPyweb#qvT{r4&2yU*^g6XG` zkR&au%4=gJ%-ma{8UdKQ6r33`4{E_yDA`+qWPra$PDy>t11mg-wX7Jvk))XmM6K)8yWXV&8f-08Rq-?$@7B1z*qK0Qi>Ou`+ zQbjvF!vG@{vNPrkAhF7f(MT`RMfLykkvPo*7t-VTN5YJ5C8%|00;m9 zAOHk_01yBIKmZ5;fyYOHZRzy(Gt8`uF&_Mz(+%NZcyuTf8VW`DP;h#5VmdrNJU%}D MN(<5R@QeBVKSNV%+W-In literal 0 HcmV?d00001 From bcf75a197d908d4c4f47c26da3a6ba54e272d55f Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 12:47:01 +0100 Subject: [PATCH 3/8] testing grid alignment at two UTM zones --- tests/mapper/test_grid_alignment.py | 113 ++++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) create mode 100644 tests/mapper/test_grid_alignment.py diff --git a/tests/mapper/test_grid_alignment.py b/tests/mapper/test_grid_alignment.py new file mode 100644 index 0000000..f7fffa3 --- /dev/null +++ b/tests/mapper/test_grid_alignment.py @@ -0,0 +1,113 @@ +''' +Ensure that all scenes in a SceneCollection have the same spatial +extent, pixel size and number of rows and columns. We test this +behavior for Sentinel-2 scenes at the boundary of two UTM zones. +''' + +import geopandas as gpd +import pytest + +from datetime import datetime +from eodal.config import get_settings +from eodal.core.raster import RasterCollection +from eodal.core.sensors import Sentinel2 +from eodal.mapper.feature import Feature +from eodal.mapper.filter import Filter +from eodal.mapper.mapper import Mapper, MapperConfigs + +settings = get_settings() +settings.USE_STAC = True + + +def test_grid_alignment(get_polygons_utm_zones): + time_start = datetime(2018, 3, 2) + time_end = datetime(2018, 4, 6) + metadata_filters = [ + Filter('cloudy_pixel_percentage', '<', 70), + Filter('processing_level', '==', 'Level-2A') + ] + feature = Feature.from_geoseries( + gds=gpd.read_file(get_polygons_utm_zones()).geometry) + mapper_configs = MapperConfigs( + collection='sentinel2-msi', + time_start=time_start, + time_end=time_end, + feature=feature, + metadata_filters=metadata_filters + ) + + mapper = Mapper(mapper_configs) + mapper.query_scenes() + + def resample(ds: RasterCollection, **kwargs): + return ds.resample(inplace=False, **kwargs) + + scene_kwargs = { + 'scene_constructor': Sentinel2.from_safe, + 'scene_constructor_kwargs': { + 'band_selection': ['red', 'red_edge_1'], + 'apply_scaling': False + }, + 'scene_modifier': resample, + 'scene_modifier_kwargs': {'target_resolution': 10} + } + mapper.load_scenes(scene_kwargs=scene_kwargs) + + # all scenes should have the same extent and, since, we resampled + # all bands to a common spatial extent, the same pixel size and number + # of rows and columns + scoll = mapper.data + + pixres_x, pixres_y, ulx, uly, nrows, ncols = [], [], [], [], [], [] + for _, scene in scoll: + for _, band in scene: + geo_info = band.geo_info + pixres_x.append(geo_info.pixres_x) + pixres_y.append(geo_info.pixres_y) + ulx.append(geo_info.ulx) + uly.append(geo_info.uly) + nrows.append(band.nrows) + ncols.append(band.ncols) + + assert len(set(pixres_x)) == 1 + assert len(set(pixres_y)) == 1 + assert len(set(ulx)) == 1 + assert len(set(uly)) == 1 + assert len(set(nrows)) == 1 + assert len(set(ncols)) == 1 + + # this should also work when the bands are not bandstacked + # i.e., have a different pixel size + scene_kwargs = { + 'scene_constructor': Sentinel2.from_safe, + 'scene_constructor_kwargs': { + 'band_selection': ['red', 'red_edge_1'], + 'apply_scaling': False + } + } + mapper.load_scenes(scene_kwargs=scene_kwargs) + + scoll = mapper.data + + band_names = scoll[scoll.timestamps[0]].band_names + # the upper left corner must be always the same + ulx, uly = [], [] + for band_name in band_names: + pixres_x, pixres_y, nrows, ncols = [], [], [], [] + for _, scene in scoll: + band = scene[band_name] + geo_info = band.geo_info + pixres_x.append(geo_info.pixres_x) + pixres_y.append(geo_info.pixres_y) + ulx.append(geo_info.ulx) + uly.append(geo_info.uly) + nrows.append(band.nrows) + ncols.append(band.ncols) + + assert len(set(pixres_x)) == 1 + assert len(set(pixres_y)) == 1 + assert len(set(nrows)) == 1 + assert len(set(ncols)) == 1 + + assert len(set(ulx)) == 1 + assert len(set(uly)) == 1 From 56dae5b0ab54a48adf0f0d1336fa1bdb26450064 Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 12:47:20 +0100 Subject: [PATCH 4/8] implemented heuristic for reference scene (#90) --- eodal/mapper/mapper.py | 70 +++++++++++++++++++++++++++--------------- 1 file changed, 45 insertions(+), 25 deletions(-) diff --git a/eodal/mapper/mapper.py b/eodal/mapper/mapper.py index f3cd1bc..79eba5c 100644 --- a/eodal/mapper/mapper.py +++ b/eodal/mapper/mapper.py @@ -49,6 +49,7 @@ from eodal.metadata.database.querying import find_raw_data_by_bbox from eodal.metadata.utils import reconstruct_path from eodal.utils.exceptions import STACError +from blinker._utilities import reference settings = get_settings() logger = settings.logger @@ -705,12 +706,34 @@ def _load_scenes_collection( )) bounds_gdf = pd.concat(bounds_list) total_bounds = bounds_gdf.total_bounds + + # ..versionadd 0.2.3 (https://github.com/EOA-team/eodal/issues/90) + # we need to ensure that all scenes have not only the same extent but + # also the same pixel size and, hence, number of rows and columns + # The heuristic is: + # 1) is there a since that was not reprojected (_duplicated is False) + # 2) If no: use the first scene as reference scene + # 3) If yes: use the first scene that was not reprojected + + # is there at a least a single scene that was not reprojected + if not self.metadata._duplicated.all(): + reference_time = self.metadata[ + ~self.metadata._duplicated][self.time_column].values[0] + timestamps = scoll.timestamps + reference_timestamp_idx = np.argmin( + [pd.to_datetime(x) - reference_time for x in timestamps]) + reference_scene = scoll[ + timestamps[reference_timestamp_idx]] + else: + reference_scene = scoll[scoll.timestamps[0]] # loop over scenes and project them onto the total bounds for _, scene in scoll: - if scene.is_bandstack(): - # get a band of the SceneCollection - reference_band = scene[scene.band_names[0]] + # loop over the bands and reproject them + # onto the target grid + for band_name, band in scene: + # get the reference band + reference_band = reference_scene[band_name] geo_info = reference_band.geo_info # get the transform of the destination extent minx, maxy = total_bounds[0], total_bounds[-1] @@ -721,28 +744,25 @@ def _load_scenes_collection( d=0, e=geo_info.pixres_y, f=maxy) - - # loop over the bands and reproject them - # onto the target grid - for _, band in scene: - dst_shape = (max(nrows), max(ncols)) - destination = np.zeros( - dst_shape, - dtype=band.values.dtype - ) - # determine nodata - if not np.isnan(band.nodata): - dst_nodata = band.nodata - else: - # if band nodata is NaN we set - # no-data to None (rasterio default) - dst_nodata = None - band.reproject( - inplace=True, - target_crs=reference_band.crs, - dst_transform=dst_transform, - destination=destination, - dst_nodata=dst_nodata) + + dst_shape = (max(nrows), max(ncols)) + destination = np.zeros( + dst_shape, + dtype=band.values.dtype + ) + # determine nodata + if not np.isnan(band.nodata): + dst_nodata = band.nodata + else: + # if band nodata is NaN we set + # no-data to None (rasterio default) + dst_nodata = None + band.reproject( + inplace=True, + target_crs=reference_band.crs, + dst_transform=dst_transform, + destination=destination, + dst_nodata=dst_nodata) self.data = scoll From 984ffd7b3e9fc2d666d940192f09d9306c6d44c8 Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 15:44:40 +0100 Subject: [PATCH 5/8] removed invalid import statement --- eodal/mapper/mapper.py | 1 - 1 file changed, 1 deletion(-) diff --git a/eodal/mapper/mapper.py b/eodal/mapper/mapper.py index 79eba5c..30448f9 100644 --- a/eodal/mapper/mapper.py +++ b/eodal/mapper/mapper.py @@ -49,7 +49,6 @@ from eodal.metadata.database.querying import find_raw_data_by_bbox from eodal.metadata.utils import reconstruct_path from eodal.utils.exceptions import STACError -from blinker._utilities import reference settings = get_settings() logger = settings.logger From 11f8f753392ea039245466e9f4c86e2996c3ab2e Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 16:26:43 +0100 Subject: [PATCH 6/8] updating copyright year 2022 -> 2023 --- eodal/config/settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/eodal/config/settings.py b/eodal/config/settings.py index 444fc25..402485b 100644 --- a/eodal/config/settings.py +++ b/eodal/config/settings.py @@ -7,7 +7,7 @@ The ``Settings`` class uses ``pydantic``. This means all attributes of the class can be **overwritten** using environmental variables or a `.env` file. -Copyright (C) 2022 Lukas Valentin Graf +Copyright (C) 2022, 2023 Lukas Valentin Graf This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by From 04d944554c309dfa115b55a813e4066c50dfc94e Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 16:27:11 +0100 Subject: [PATCH 7/8] moved `USE_STAC` statement inside test function --- tests/mapper/test_grid_alignment.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/tests/mapper/test_grid_alignment.py b/tests/mapper/test_grid_alignment.py index f7fffa3..a76fcfa 100644 --- a/tests/mapper/test_grid_alignment.py +++ b/tests/mapper/test_grid_alignment.py @@ -16,10 +16,13 @@ from eodal.mapper.mapper import Mapper, MapperConfigs settings = get_settings() -settings.USE_STAC = True def test_grid_alignment(get_polygons_utm_zones): + """test grid alignment for Sentinel-2""" + + settings.USE_STAC = True + time_start = datetime(2018, 3, 2) time_end = datetime(2018, 4, 6) metadata_filters = [ @@ -70,7 +73,9 @@ def resample(ds: RasterCollection, **kwargs): ncols.append(band.ncols) assert len(set(pixres_x)) == 1 + assert pixres_x[0] == 10 assert len(set(pixres_y)) == 1 + assert pixres_y[0] == -10 assert len(set(ulx)) == 1 assert len(set(uly)) == 1 assert len(set(nrows)) == 1 @@ -108,6 +113,14 @@ def resample(ds: RasterCollection, **kwargs): assert len(set(pixres_y)) == 1 assert len(set(nrows)) == 1 assert len(set(ncols)) == 1 + # only the red band should have a spatial resolution of 10 m + # as we did not resample the 20 m bands + if band_name == 'B04': + assert pixres_x[0] == 10. + assert pixres_y[0] == -10. + else: + assert pixres_x[0] == 20. + assert pixres_y[0] == -20. assert len(set(ulx)) == 1 assert len(set(uly)) == 1 From 25ff686d443ecd595e68b37bdb70031216f1a660 Mon Sep 17 00:00:00 2001 From: lukas Date: Sun, 19 Nov 2023 16:27:33 +0100 Subject: [PATCH 8/8] FIX: only call resampling function if it's really required --- eodal/mapper/mapper.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/eodal/mapper/mapper.py b/eodal/mapper/mapper.py index 30448f9..413739b 100644 --- a/eodal/mapper/mapper.py +++ b/eodal/mapper/mapper.py @@ -453,11 +453,18 @@ def _process_scene( scene = scene_modifier.__call__(scene, **scene_modifier_kwargs) # reproject scene if necessary - scene.reproject( - target_crs=item.target_epsg, - interpolation_method=reprojection_method, - inplace=True - ) + # ..versionadd 0.2.3 check for the EPSG to make sure no unnecessary + # operations are undertaken to save runtime and avoid floating + # point inaccuracies + epsg_scene = [b.geo_info.epsg for _, b in scene] + intersection = set(epsg_scene).intersection(set([item.target_epsg])) + # we need to reproject only if the intersection returns an empty set. + if len(intersection) == 0: + scene.reproject( + target_crs=item.target_epsg, + interpolation_method=reprojection_method, + inplace=True + ) return scene