diff --git a/python/idsse_common/idsse/common/grid_proj.py b/python/idsse_common/idsse/common/grid_proj.py index 6880c9a2..282d98d8 100644 --- a/python/idsse_common/idsse/common/grid_proj.py +++ b/python/idsse_common/idsse/common/grid_proj.py @@ -1,14 +1,17 @@ """Module that wraps pyproj (cartographic) library and transforms objects into other data forms""" -# ------------------------------------------------------------------------------- +# ---------------------------------------------------------------------------------- # Created on Mon Jul 31 2023 # # Copyright (c) 2023 Colorado State University. All rights reserved. (1) +# Copyright (c) 2023 Regents of the University of Colorado. All rights reserved. (2) # # Contributors: # Mackenzie Grimes (1) +# Geary Layne (2) # -# ------------------------------------------------------------------------------- +# ---------------------------------------------------------------------------------- # pylint: disable=invalid-name +# cspell:word fliplr, flipud from typing import Self, Tuple, Union, Optional from enum import Enum @@ -24,10 +27,15 @@ class RoundingMethod(Enum): - """Transformations to apply to calculated pixel values when casting to ints""" + """Transformations indicators to be applied to pixel values when casting to ints""" ROUND = 'ROUND' FLOOR = 'FLOOR' +class Flip(Enum): + """Flip axis indicators to be applied flipping pixel orientation""" + BOTH = 'BOTH' + HORIZONTAL = 'HORIZONTAL' + VERTICAL = 'VERTICAL' class GridProj: """ @@ -36,8 +44,8 @@ class GridProj: """ def __init__(self, crs: CRS, - latitude: float, - longitude: float, + lower_left_lat: float, + lower_left_lon: float, width: float, height: float, dx: float, @@ -46,8 +54,8 @@ def __init__(self, self._trans = Transformer.from_crs(crs.geodetic_crs, crs) self._x_offset = 0 self._y_offset = 0 - if latitude is not None: - self._x_offset, self._y_offset = self._trans.transform(longitude, latitude) + if lower_left_lat is not None: + self._x_offset, self._y_offset = self._trans.transform(lower_left_lon, lower_left_lat) self._w = width self._h = height self._dx = dx @@ -79,18 +87,58 @@ def from_proj_grid_spec(cls, proj_string: str, grid_string: str) -> Self: int(grid_args['w']), int(grid_args['h']), grid_args['dx'], grid_args['dy']) - def map_proj_to_pixel( + @property + def width(self): + """Provides access grid space width""" + return self._w + + @property + def height(self): + """Provides access grid space height""" + return self._h + + @property + def shape(self): + """Provides access grid space shape: (width, height)""" + return self._w, self._h + + def flip(self, flip: Flip = Flip.BOTH): + """Reverse the order of pixels for a given orientation + + Args: + flip (Flip): The default, flip=BOTH, will flip over both axes. + """ + if flip is Flip.BOTH: + self._x_offset, self._y_offset = self.map_pixel_to_crs(self.width, self.height) + self._dx *= -1 + self._dy *= -1 + elif flip is Flip.HORIZONTAL: + self._x_offset, _ = self.map_pixel_to_crs(self.width, 0) + self._dx *= -1 + elif flip is Flip.VERTICAL: + _, self._y_offset = self.map_pixel_to_crs(0, self.height) + self._dy *= -1 + + def fliplr(self): + """Reverse the order of the pixels left to right""" + self.flip(Flip.HORIZONTAL) + + def flipud(self): + """Reverse the order of the pixels up to down""" + self.flip(Flip.VERTICAL) + + def map_geo_to_pixel( self, - x: float, - y: float, + lon: float, + lat: float, rounding: Optional[Union[str, RoundingMethod]] = None, precision: int = 0 ) -> Pixel: """Map projection to a pixel. Args: - x (float): x geographic coordinate - y (float): y geographic coordinate + lon (float): x geographic coordinate + lat (float): y geographic coordinate rounding (Optional[Union[str, RoundingMethod]]): ROUND to apply round_half_away() to pixel values, FLOOR to apply math.floor(). @@ -103,43 +151,75 @@ def map_proj_to_pixel( Pixel: x, y values of pixel, rounded to ints if rounding is passed, otherwise floats """ # pylint: disable=not-an-iterable - return self.map_geo_to_pixel( - *self._trans.transform(x, y), + return self.map_crs_to_pixel( + *self._trans.transform(lon, lat), rounding=rounding, precision=precision ) - def map_pixel_to_proj(self, x: float, y: float) -> Tuple[float, float]: - """Map pixel to a projection""" + def map_pixel_to_geo(self, x: float, y: float) -> Tuple[float, float]: + """Map pixel x,y to a projection + + Args: + x (float): x coordinate in pixel space + y (float): y coordinate in pixel space + + Returns: + Tuple[float, float]: Geographic coordinate as lon,lat + """ return self._trans.transform( - *self.map_pixel_to_geo(x, y), + *self.map_pixel_to_crs(x, y), direction=TransformDirection.INVERSE ) - def map_proj_to_geo(self, x, y) -> Tuple[float, float]: - """Map projection to geographic coordinates""" - return self._trans.transform(x, y) + def map_geo_to_crs(self, lon: float, lat: float) -> Tuple[float, float]: + """Map geographic coordinate (lon, lat) to CRS - def map_pixel_to_geo(self, x: float, y: float) -> Tuple[float, float]: - """Map pixel to geographic coordinates""" + Args: + lon (float): x geographic coordinate + lat (float): y geographic coordinate + + Returns: + Tuple[float, float]: Coordinate Reference System + """ + return self._trans.transform(lon, lat) + + def map_pixel_to_crs(self, x: float, y: float) -> Tuple[float, float]: + """Map pixel space (x,y) to Coordinate Reference System + + Args: + x (float): x coordinate in pixel space + y (float): y coordinate in pixel space + + Returns: + Tuple[float, float]: Coordinate Reference System coordinate + """ return x * self._dx + self._x_offset, y * self._dy + self._y_offset - def map_geo_to_proj(self, x: float, y: float) -> Tuple[float, float]: - """Map geographic coordinates to a projection""" + def map_crs_to_geo(self, x: float, y: float) -> Tuple[float, float]: + """Map Coordinate Reference System (x,y) to Geographical space (lon,lat) + + Args: + x (float): x coordinate in CRS space + y (float): y coordinate in CRS space + + Returns: + Tuple[float, float]: Geographic coordinate as lon,lat + """ return self._trans.transform(x, y, direction=TransformDirection.INVERSE) - def map_geo_to_pixel( + def map_crs_to_pixel( self, x: float, y: float, rounding: Optional[Union[str, RoundingMethod]] = None, precision: int = 0 ) -> Pixel: - """Map geographic coordinates to pixel x and y + """Map Coordinate Reference System (x,y) coordinates to pixel x and y Args: - x (float): x geographic coordinate - y (float): y geographic coordinate + x (float): x CRS coordinate + y (float): y CRS coordinate rounding (Optional[Union[str, RoundingMethod]]): ROUND to apply round_half_away() to pixel values, FLOOR to apply math.floor(). diff --git a/python/idsse_common/test/test_grid_proj.py b/python/idsse_common/test/test_grid_proj.py index 78adae7d..7fb0216a 100644 --- a/python/idsse_common/test/test_grid_proj.py +++ b/python/idsse_common/test/test_grid_proj.py @@ -1,13 +1,15 @@ """Test suite for grid_proj.py""" -# -------------------------------------------------------------------------------- +# ---------------------------------------------------------------------------------- # Created on Wed Aug 2 2023 # # Copyright (c) 2023 Colorado State University. All rights reserved. (1) +# Copyright (c) 2023 Regents of the University of Colorado. All rights reserved. (2) # # Contributors: # Mackenzie Grimes (1) +# Geary Layne (2) # -# -------------------------------------------------------------------------------- +# ---------------------------------------------------------------------------------- # pylint: disable=missing-function-docstring,redefined-outer-name,invalid-name,protected-access import math @@ -31,17 +33,20 @@ ) GRID_SPEC_WITHOUT_LOWER_LEFT = '+dx=2539.703 +dy=2539.703 +w=2345 +h=1597' +WIDTH = 2345 +HEIGHT = 1597 + EXAMPLE_PIXELS = [ (0, 0), (0, 1), (2000, 1500) ] -EXAMPLE_PROJECTIONS = [ +EXAMPLE_LON_LAT = [ (-126.2766, 19.229), (-126.28209649530142, 19.251224896946418), (-71.12860760432866, 54.09108720984939) ] -EXAMPLE_GEOS = [ +EXAMPLE_CRS = [ (-3275807.350733357, -260554.63043505285), (-3275807.3507333556, -258014.92743505232), (1803598.6492666428, 3548999.8695649463) @@ -63,7 +68,7 @@ def grid_proj() -> GridProj: def test_from_proj_grid_spec(grid_proj: GridProj): assert isinstance(grid_proj, GridProj) - assert (grid_proj._x_offset, grid_proj._y_offset) == approx_tuple(EXAMPLE_GEOS[0]) + assert (grid_proj._x_offset, grid_proj._y_offset) == approx_tuple(EXAMPLE_CRS[0]) assert grid_proj._h == 1597 assert grid_proj._w == 2345 assert grid_proj._dx == approx(2539.703) @@ -78,24 +83,56 @@ def test_from_proj_grid_spec_with_offset(): proj_with_offset = GridProj.from_proj_grid_spec(PROJ_SPEC_WITH_OFFSET, GRID_SPEC_WITHOUT_LOWER_LEFT) - proj_xy = proj_with_offset.map_pixel_to_proj(*EXAMPLE_PIXELS[0]) - assert proj_xy == approx_tuple(EXAMPLE_PROJECTIONS[0]) + proj_xy = proj_with_offset.map_pixel_to_geo(*EXAMPLE_PIXELS[0]) + assert proj_xy == approx_tuple(EXAMPLE_LON_LAT[0]) + + +# test properties +def test_get_width(grid_proj: GridProj): + assert grid_proj.width == WIDTH + + +def test_get_height(grid_proj: GridProj): + assert grid_proj.height == HEIGHT + + +def test_get_shape(grid_proj: GridProj): + assert grid_proj.shape == (WIDTH, HEIGHT) + + +# test flips +def test_fliplr(grid_proj: GridProj): + lower_left_lon_lat = grid_proj.map_pixel_to_geo(0, 0) + grid_proj.fliplr() + assert grid_proj.map_pixel_to_geo(grid_proj.width, 0) == lower_left_lon_lat + + +def test_flipud(grid_proj: GridProj): + lower_left_lon_lat = grid_proj.map_pixel_to_geo(0, 0) + grid_proj.flipud() + assert grid_proj.map_pixel_to_geo(0, grid_proj.height) == lower_left_lon_lat + + +def test_flip_both_lr_ud(grid_proj: GridProj): + lower_left_lon_lat = grid_proj.map_pixel_to_geo(0, 0) + grid_proj.flip() + assert grid_proj.map_pixel_to_geo(grid_proj.width, grid_proj.height) == lower_left_lon_lat # transformation methods testing -def test_map_proj_to_pixel_round_half_away(grid_proj: GridProj): - for index, proj in enumerate(EXAMPLE_PROJECTIONS): - pixel_xy = grid_proj.map_proj_to_pixel( - *proj, +def test_map_crs_to_pixel_round_half_away(grid_proj: GridProj): + for index, crs_xy in enumerate(EXAMPLE_CRS): + pixel_xy = grid_proj.map_crs_to_pixel( + *crs_xy, rounding=RoundingMethod.ROUND ) assert pixel_xy == EXAMPLE_PIXELS[index] -def test_map_proj_to_pixel_round_floor(grid_proj: GridProj): - for index, proj in enumerate(EXAMPLE_PROJECTIONS): - i, j = grid_proj.map_proj_to_pixel( - *proj, +def test_map_crs_to_pixel_round_floor(grid_proj: GridProj): + for index, crs_xy in enumerate(EXAMPLE_CRS): + i, j = grid_proj.map_crs_to_pixel( + *crs_xy, rounding=RoundingMethod.FLOOR ) # due to math imprecision internal to pyproj.transform(), some test results are a bit @@ -103,77 +140,77 @@ def test_map_proj_to_pixel_round_floor(grid_proj: GridProj): assert (approx(i, abs=1), approx(j, abs=1)) == EXAMPLE_PIXELS[index] -def test_map_proj_to_geo(grid_proj: GridProj): - for index, proj in enumerate(EXAMPLE_PROJECTIONS): - geo_xy = grid_proj.map_proj_to_geo(*proj) - assert geo_xy == approx_tuple(EXAMPLE_GEOS[index]) +def test_map_geo_to_crs(grid_proj: GridProj): + for index, lon_lat in enumerate(EXAMPLE_LON_LAT): + geo_xy = grid_proj.map_geo_to_crs(*lon_lat) + assert geo_xy == approx_tuple(EXAMPLE_CRS[index]) -def test_map_pixel_to_geo(grid_proj: GridProj): +def test_map_pixel_to_crs(grid_proj: GridProj): for index, pixel in enumerate(EXAMPLE_PIXELS): - geo_x, geo_y = grid_proj.map_pixel_to_geo(*pixel) - assert (geo_x, geo_y) == approx_tuple(EXAMPLE_GEOS[index]) + geo_x, geo_y = grid_proj.map_pixel_to_crs(*pixel) + assert (geo_x, geo_y) == approx_tuple(EXAMPLE_CRS[index]) -def test_map_pixel_to_proj(grid_proj: GridProj): +def test_map_pixel_to_geo(grid_proj: GridProj): for index, pixel in enumerate(EXAMPLE_PIXELS): - proj_x, proj_y = grid_proj.map_pixel_to_proj(*pixel) - assert (proj_x, proj_y) == approx_tuple(EXAMPLE_PROJECTIONS[index]) + proj_x, proj_y = grid_proj.map_pixel_to_geo(*pixel) + assert (proj_x, proj_y) == approx_tuple(EXAMPLE_LON_LAT[index]) -def test_map_geo_to_proj(grid_proj: GridProj): - for index, geo in enumerate(EXAMPLE_GEOS): - proj_x, proj_y = grid_proj.map_geo_to_proj(*geo) - assert (proj_x, proj_y) == approx_tuple(EXAMPLE_PROJECTIONS[index]) +def test_map_crs_to_geo(grid_proj: GridProj): + for index, geo in enumerate(EXAMPLE_CRS): + proj_x, proj_y = grid_proj.map_crs_to_geo(*geo) + assert (proj_x, proj_y) == approx_tuple(EXAMPLE_LON_LAT[index]) -def test_geo_to_pixel_no_rounding(grid_proj: GridProj): - for index, geo in enumerate(EXAMPLE_GEOS): - i, j = grid_proj.map_geo_to_pixel(*geo) +def test_crs_to_pixel_no_rounding(grid_proj: GridProj): + for index, geo in enumerate(EXAMPLE_CRS): + i, j = grid_proj.map_crs_to_pixel(*geo) # round result, which will not be precisely the integer that was passed assert (round_half_away(i, 6), round_half_away(j, 6)) == EXAMPLE_PIXELS[index] -def test_geo_to_pixel_floor(grid_proj: GridProj): - for index, geo in enumerate(EXAMPLE_GEOS): - floor_ij = grid_proj.map_geo_to_pixel(*geo, RoundingMethod.FLOOR) - i, j = grid_proj.map_geo_to_pixel(*geo) +def test_crs_to_pixel_floor(grid_proj: GridProj): + for index, geo in enumerate(EXAMPLE_CRS): + floor_ij = grid_proj.map_crs_to_pixel(*geo, RoundingMethod.FLOOR) + i, j = grid_proj.map_crs_to_pixel(*geo) assert (math.floor(i), math.floor(j)) == floor_ij assert (round_half_away(i), round_half_away(j)) == EXAMPLE_PIXELS[index] -def test_geo_to_pixel_round(grid_proj: GridProj): - for index, geo in enumerate(EXAMPLE_GEOS): - i, j = grid_proj.map_geo_to_pixel(*geo, RoundingMethod.ROUND) +def test_crs_to_pixel_round(grid_proj: GridProj): + for index, geo in enumerate(EXAMPLE_CRS): + i, j = grid_proj.map_crs_to_pixel(*geo, RoundingMethod.ROUND) assert (i, j) == EXAMPLE_PIXELS[index] -def test_geo_to_pixel_round_str(grid_proj: GridProj): - i, j = grid_proj.map_geo_to_pixel(*EXAMPLE_GEOS[0], 'round') +def test_crs_to_pixel_round_str(grid_proj: GridProj): + i, j = grid_proj.map_crs_to_pixel(*EXAMPLE_CRS[0], 'round') assert (i, j) == EXAMPLE_PIXELS[0] - i, j = grid_proj.map_geo_to_pixel(*EXAMPLE_GEOS[1], 'ROUND') + i, j = grid_proj.map_crs_to_pixel(*EXAMPLE_CRS[1], 'ROUND') assert (i, j) == EXAMPLE_PIXELS[1] def test_compound_transformations_stay_consistent(grid_proj: GridProj): # start with pixel, convert to projection initial_pixel = (EXAMPLE_PIXELS[2][0], EXAMPLE_PIXELS[2][1]) - proj_x, proj_y = grid_proj.map_pixel_to_proj(*initial_pixel) - initial_proj = (proj_x, proj_y) + proj_x, proj_y = grid_proj.map_pixel_to_geo(*initial_pixel) + initial_geo = (proj_x, proj_y) # convert projection to geographic coordinates - geo_x, geo_y = grid_proj.map_proj_to_geo(proj_x, proj_y) - initial_geo = geo_x, geo_y + geo_x, geo_y = grid_proj.map_geo_to_crs(proj_x, proj_y) + initial_crs = geo_x, geo_y # convert geographic coordinates back to pixel, full circle, and data should be unchanged - pixel_x, pixel_y = grid_proj.map_geo_to_pixel(geo_x, geo_y) + pixel_x, pixel_y = grid_proj.map_crs_to_pixel(geo_x, geo_y) assert (round_half_away(pixel_x, 6), round_half_away(pixel_y, 6)) == initial_pixel # convert pixel back to geographic coordinates - geo_x, geo_y = grid_proj.map_pixel_to_geo(pixel_x, pixel_y) - assert (geo_x, geo_y) == approx_tuple(initial_geo) + geo_x, geo_y = grid_proj.map_pixel_to_crs(pixel_x, pixel_y) + assert (geo_x, geo_y) == approx_tuple(initial_crs) # convert geographic coordinates back to projection - proj_x, proj_y = grid_proj.map_geo_to_proj(geo_x, geo_y) - assert (proj_x, proj_y) == approx_tuple(initial_proj) + proj_x, proj_y = grid_proj.map_crs_to_geo(geo_x, geo_y) + assert (proj_x, proj_y) == approx_tuple(initial_geo)