From ac5d4c7a9c7ed6c8d41f1d6e6d18b069b7842dd5 Mon Sep 17 00:00:00 2001 From: George Spill Date: Fri, 26 Jul 2024 11:31:03 +0100 Subject: [PATCH 1/5] Update point to handle inverted latitude --- rio_tiler/reader.py | 39 ++++++++++++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 9 deletions(-) diff --git a/rio_tiler/reader.py b/rio_tiler/reader.py index 5ac6db53..739692fb 100644 --- a/rio_tiler/reader.py +++ b/rio_tiler/reader.py @@ -96,7 +96,9 @@ def read( resampling_method: RIOResampling = "nearest", reproject_method: WarpResampling = "nearest", unscale: bool = False, - post_process: Optional[Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray]] = None, + post_process: Optional[ + Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray] + ] = None, ) -> ImageData: """Low level read function. @@ -199,7 +201,9 @@ def read( values = dataset.read( indexes=indexes, window=window, - out_shape=(len(indexes), height, width) if height and width else None, + out_shape=( + (len(indexes), height, width) if height and width else None + ), resampling=io_resampling, boundless=boundless, ) @@ -247,7 +251,9 @@ def read( stats = [] for ix in indexes: tags = dataset.tags(ix) - if all(stat in tags for stat in ["STATISTICS_MINIMUM", "STATISTICS_MAXIMUM"]): + if all( + stat in tags for stat in ["STATISTICS_MINIMUM", "STATISTICS_MAXIMUM"] + ): stat_min = float(tags.get("STATISTICS_MINIMUM")) stat_max = float(tags.get("STATISTICS_MAXIMUM")) stats.append((stat_min, stat_max)) @@ -310,7 +316,9 @@ def part( resampling_method: RIOResampling = "nearest", reproject_method: WarpResampling = "nearest", unscale: bool = False, - post_process: Optional[Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray]] = None, + post_process: Optional[ + Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray] + ] = None, ) -> ImageData: """Read part of a dataset. @@ -358,8 +366,12 @@ def part( src_bounds = transform_bounds( src_dst.crs, dst_crs, *src_dst.bounds, densify_pts=21 ) - x_overlap = max(0, min(src_bounds[2], bounds[2]) - max(src_bounds[0], bounds[0])) - y_overlap = max(0, min(src_bounds[3], bounds[3]) - max(src_bounds[1], bounds[1])) + x_overlap = max( + 0, min(src_bounds[2], bounds[2]) - max(src_bounds[0], bounds[0]) + ) + y_overlap = max( + 0, min(src_bounds[3], bounds[3]) - max(src_bounds[1], bounds[1]) + ) cover_ratio = (x_overlap * y_overlap) / ( (bounds[2] - bounds[0]) * (bounds[3] - bounds[1]) ) @@ -507,7 +519,9 @@ def point( resampling_method: RIOResampling = "nearest", reproject_method: WarpResampling = "nearest", unscale: bool = False, - post_process: Optional[Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray]] = None, + post_process: Optional[ + Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray] + ] = None, ) -> PointData: """Read a pixel value for a point. @@ -559,9 +573,16 @@ def point( xs, ys = transform_coords(coord_crs, dataset.crs, [lon], [lat]) lon, lat = xs[0], ys[0] + dataset_min_lon, dataset_min_lat, dataset_max_lon, dataset_max_lat = ( + dataset.bounds + ) + # check if latitude is inverted + dataset_min_lat, dataset_max_lat = min(dataset_min_lat, dataset_max_lat), max( + dataset_min_lat, dataset_max_lat + ) if not ( - (dataset.bounds[0] < lon < dataset.bounds[2]) - and (dataset.bounds[1] < lat < dataset.bounds[3]) + (dataset_min_lon < lon < dataset_max_lon) + and (dataset_min_lat < lat < dataset_max_lat) ): raise PointOutsideBounds("Point is outside dataset bounds") From cc93e4f0dc77ebcc261787336add9cc9bc13c0b1 Mon Sep 17 00:00:00 2001 From: George Spill Date: Fri, 26 Jul 2024 11:44:41 +0100 Subject: [PATCH 2/5] add formatting --- rio_tiler/reader.py | 29 +++++++++-------------------- 1 file changed, 9 insertions(+), 20 deletions(-) diff --git a/rio_tiler/reader.py b/rio_tiler/reader.py index 739692fb..18f410b1 100644 --- a/rio_tiler/reader.py +++ b/rio_tiler/reader.py @@ -96,9 +96,7 @@ def read( resampling_method: RIOResampling = "nearest", reproject_method: WarpResampling = "nearest", unscale: bool = False, - post_process: Optional[ - Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray] - ] = None, + post_process: Optional[Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray]] = None, ) -> ImageData: """Low level read function. @@ -251,9 +249,7 @@ def read( stats = [] for ix in indexes: tags = dataset.tags(ix) - if all( - stat in tags for stat in ["STATISTICS_MINIMUM", "STATISTICS_MAXIMUM"] - ): + if all(stat in tags for stat in ["STATISTICS_MINIMUM", "STATISTICS_MAXIMUM"]): stat_min = float(tags.get("STATISTICS_MINIMUM")) stat_max = float(tags.get("STATISTICS_MAXIMUM")) stats.append((stat_min, stat_max)) @@ -316,9 +312,7 @@ def part( resampling_method: RIOResampling = "nearest", reproject_method: WarpResampling = "nearest", unscale: bool = False, - post_process: Optional[ - Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray] - ] = None, + post_process: Optional[Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray]] = None, ) -> ImageData: """Read part of a dataset. @@ -366,12 +360,8 @@ def part( src_bounds = transform_bounds( src_dst.crs, dst_crs, *src_dst.bounds, densify_pts=21 ) - x_overlap = max( - 0, min(src_bounds[2], bounds[2]) - max(src_bounds[0], bounds[0]) - ) - y_overlap = max( - 0, min(src_bounds[3], bounds[3]) - max(src_bounds[1], bounds[1]) - ) + x_overlap = max(0, min(src_bounds[2], bounds[2]) - max(src_bounds[0], bounds[0])) + y_overlap = max(0, min(src_bounds[3], bounds[3]) - max(src_bounds[1], bounds[1])) cover_ratio = (x_overlap * y_overlap) / ( (bounds[2] - bounds[0]) * (bounds[3] - bounds[1]) ) @@ -519,9 +509,7 @@ def point( resampling_method: RIOResampling = "nearest", reproject_method: WarpResampling = "nearest", unscale: bool = False, - post_process: Optional[ - Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray] - ] = None, + post_process: Optional[Callable[[numpy.ma.MaskedArray], numpy.ma.MaskedArray]] = None, ) -> PointData: """Read a pixel value for a point. @@ -577,8 +565,9 @@ def point( dataset.bounds ) # check if latitude is inverted - dataset_min_lat, dataset_max_lat = min(dataset_min_lat, dataset_max_lat), max( - dataset_min_lat, dataset_max_lat + dataset_min_lat, dataset_max_lat = ( + min(dataset_min_lat, dataset_max_lat), + max(dataset_min_lat, dataset_max_lat), ) if not ( (dataset_min_lon < lon < dataset_max_lon) From 2dd04165f553470532ab635727b1399533cde6d2 Mon Sep 17 00:00:00 2001 From: George Spill Date: Fri, 26 Jul 2024 14:06:15 +0100 Subject: [PATCH 3/5] add warning and test --- rio_tiler/reader.py | 6 ++++++ tests/test_reader.py | 8 ++++++++ 2 files changed, 14 insertions(+) diff --git a/rio_tiler/reader.py b/rio_tiler/reader.py index 18f410b1..50e4ad46 100644 --- a/rio_tiler/reader.py +++ b/rio_tiler/reader.py @@ -565,6 +565,12 @@ def point( dataset.bounds ) # check if latitude is inverted + if dataset_min_lat > dataset_max_lat: + warnings.warn( + "BoundingBox of the dataset is inverted (minLat > maxLat).", + UserWarning, + ) + dataset_min_lat, dataset_max_lat = ( min(dataset_min_lat, dataset_max_lat), max(dataset_min_lat, dataset_max_lat), diff --git a/tests/test_reader.py b/tests/test_reader.py index 9f996d1c..1bb3539e 100644 --- a/tests/test_reader.py +++ b/tests/test_reader.py @@ -35,6 +35,7 @@ COG_NODATA_FLOAT_NAN = os.path.join( os.path.dirname(__file__), "fixtures", "cog_nodata_float_nan.tif" ) +COG_INVERTED = os.path.join(PREFIX, "inverted_lat.tif") @pytest.fixture(autouse=True) @@ -850,3 +851,10 @@ def test_tile_read_nodata_float(): prev = reader.read(src_dst, max_size=100) assert prev.mask[0, 0] == 0 assert not numpy.all(prev.mask) + + +def test_inverted_latitude_point(): + """Make sure we can read a point from a file with inverted latitude.""" + with rasterio.open(COG_INVERTED) as src_dst: + pt = reader.point(src_dst, [-104.77519499, 38.95367054]) + assert pt.data[0] == -9999.0 From 4c1e4aa0653812833d54ab7f2b87a7556054d5a4 Mon Sep 17 00:00:00 2001 From: George Spill Date: Fri, 26 Jul 2024 14:22:02 +0100 Subject: [PATCH 4/5] fix test --- tests/test_reader.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_reader.py b/tests/test_reader.py index 1bb3539e..ece1186a 100644 --- a/tests/test_reader.py +++ b/tests/test_reader.py @@ -35,7 +35,7 @@ COG_NODATA_FLOAT_NAN = os.path.join( os.path.dirname(__file__), "fixtures", "cog_nodata_float_nan.tif" ) -COG_INVERTED = os.path.join(PREFIX, "inverted_lat.tif") +COG_INVERTED = os.path.join(os.path.dirname(__file__), "fixtures", "inverted_lat.tif") @pytest.fixture(autouse=True) From c85946b3992096d0602dfce6477f29b59e7faf2f Mon Sep 17 00:00:00 2001 From: vincentsarago Date: Wed, 31 Jul 2024 16:14:10 +0200 Subject: [PATCH 5/5] update changelog --- CHANGES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGES.md b/CHANGES.md index b27c12e2..868db7cd 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,6 +1,7 @@ # unreleased * better error message for `TileOutsideBounds` errors (author @abarciauskas-bgse, https://github.com/cogeotiff/rio-tiler/pull/712) +* handle of inverted latitude in `reader.point` (author @georgespill, https://github.com/cogeotiff/rio-tiler/pull/716) # 6.6.1 (2024-05-17)