-
Notifications
You must be signed in to change notification settings - Fork 0
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
feat: IDSSE-256: GridProj #5
Changes from 3 commits
b2c059b
6fdc41c
58bbd45
5d216eb
c60d13b
89d7611
9a4c6eb
e6776d9
6b7350f
e7c7be4
9578861
8aa1685
1dd8a63
252f9c9
8c574f0
d60c7ce
9de285e
63618b0
ce4c78d
b287da4
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,86 @@ | ||
"""Module that wraps pyproj (cartographic/geodetic) library and transforms objects into other data forms""" | ||
# ------------------------------------------------------------------------------- | ||
# Created on Mon Jul 31 2023 | ||
# | ||
# Copyright (c) 2023 Colorado State University. All rights reserved. (1) | ||
# | ||
# Contributors: | ||
# Mackenzie Grimes (1) | ||
# | ||
# ------------------------------------------------------------------------------- | ||
# pylint: disable=invalid-name | ||
|
||
from typing import Self, Tuple, Any, Optional | ||
|
||
from pyproj import CRS, Transformer | ||
from pyproj.enums import TransformDirection | ||
|
||
|
||
class GridProj: | ||
"""Wrapper for pyproj instance with methods to convert to and from geographic coordinates, pixels, etc.""" | ||
def __init__(self, | ||
crs: CRS, | ||
latitude: float, | ||
longitude: float, | ||
width: float, | ||
height: float, | ||
dx: float, | ||
dy: Optional[float] = None): | ||
# pylint: disable=too-many-arguments,unpacking-non-sequence | ||
self._trans = Transformer.from_crs(crs.geodetic_crs, crs) | ||
self._x_offset, self._y_offset = self._trans.transform(longitude, latitude) | ||
self._w = width | ||
self._h = height | ||
self._dx = dx | ||
self._dy = dy if dy else dx | ||
|
||
@classmethod | ||
def from_proj_grid_spec(cls, proj_string: str, grid_string: str) -> Self: | ||
""" Create GridProj instance from projection grid specs | ||
|
||
Args: | ||
proj_string (str): pyproj projection string | ||
grid_string (str): Grid string | ||
|
||
Returns: | ||
Self: New instance of GridProj which can then be converted to any format | ||
""" | ||
crs = CRS.from_proj4(proj_string) | ||
grid_args = { | ||
key[1:]: float(value) | ||
for key, value in ( | ||
pair.split('=') for pair in grid_string.split(' ') | ||
) | ||
} | ||
return GridProj(crs, | ||
grid_args['lat_ll'], grid_args['lon_ll'], | ||
int(grid_args['w']), int(grid_args['h']), | ||
grid_args['dx'], grid_args['dy']) | ||
|
||
def map_proj_to_pixel(self, x, y) -> Tuple[float, float]: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It would be best of all _to_pixel functions supported rounding, so add the same logic to map_proj_to_pixel() as map_geo_to_pixel() has. It might be good to create a private function to round so that both _to_pixel functions can call the same code. |
||
"""Map projection to a pixel""" | ||
return self.map_geo_to_pixel(*self._trans.transform(x, y)) # pylint: disable=not-an-iterable | ||
|
||
def map_pixel_to_proj(self, x: float, y: float) -> Tuple[Any, Any]: | ||
"""Map pixel to a projection""" | ||
return self._trans.transform(*self.map_pixel_to_geo(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_pixel_to_geo(self, x: float, y: float) -> Tuple[float, float]: | ||
"""Map pixel to geographic coordinates""" | ||
return x * self._dx + self._x_offset, y * self._dy + self._y_offset | ||
|
||
def map_geo_to_proj(self, x: float, y: float) -> Tuple[Any, Any]: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Nice to have: the return type of a tuple of Any is less helpful than Tuple[Union[int, float], Union[int, float]]. I think that is the correct syntax. |
||
"""Map geographic coordinates to a projection""" | ||
return self._trans.transform(x, y, direction=TransformDirection.INVERSE) | ||
|
||
def map_geo_to_pixel(self, x: float, y: float) -> Tuple[int, int]: | ||
"""Map geographic coordinates to a pixel | ||
|
||
Pixels are only identified by whole numbers when graphically rendered, | ||
so the transformed numerical values (floats) can be safely rounded to ints | ||
""" | ||
return round((x - self._x_offset) / self._dx), round((y - self._y_offset) / self._dy) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @Geary-Layne it's safe to assume that pixels can be rounded integers, correct? From my general understanding of how browsers actually paint onto a display, a fraction of a pixel is not possible. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I ask because unit testing uncovered this strange precision behavior, before I made the change here to round the resulting pixel. I don't think this extra precision is meaningful (or comprehendible) for a browser client, but I could be wrong: geo_x, geo_y = self.map_pixel_to_geo(0, 1)
pixel_x, pixel_y = self.map_geo_to_pixel(geo_x, geo_y)
print(pixel_x, pixel_y) # (0.0 0.9999999999951331) There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. You are correct that for some uses, an non integer pixel value is not meaningful and is problematic. There are however some use case for this code (now that it is a general tool and not specific to the slice operator) where you would want floating point numbers. Also, sometime the use would dictate round and sometime floor. In you test code the precision error is totally executable, and you should us approx(). When I'd implement code like this class in Java and python in the past, I supported an option argument with _to_pixel() function. Something like: There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I don't feel like Enum is overkill here. It makes the comparison easy, although I didn't find an elegant way to directly map, and invoke, Enum members to functions. Simple if statements seemed sufficient for now, and default to no rounding as previously defined in the data-access-service version of GridProj. |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,133 @@ | ||
"""Test suite for grid_proj.py""" | ||
# -------------------------------------------------------------------------------- | ||
# Created on Wed Aug 2 2023 | ||
# | ||
# Copyright (c) 2023 Colorado State University. All rights reserved. (1) | ||
# | ||
# Contributors: | ||
# Mackenzie Grimes (1) | ||
# | ||
# -------------------------------------------------------------------------------- | ||
# pylint: disable=missing-function-docstring,redefined-outer-name,invalid-name,protected-access | ||
|
||
from typing import Tuple | ||
|
||
from pytest import fixture, approx | ||
|
||
from idsse.common.grid_proj import GridProj | ||
|
||
# example data | ||
EXAMPLE_PROJ_SPEC = '+proj=lcc +lat_0=25.0 +lon_0=-95.0 +lat_1=25.0 +r=6371200' | ||
PROJ_SPEC_WITH_OFFSET = ( | ||
'+proj=lcc +lat_0=25.0 +lon_0=-95.0 +lat_1=25.0 +x_0=3271.152832031251 ' | ||
'+y_0=263.7934687500001 +r=6371.2, +units=km' | ||
) | ||
EXAMPLE_GRID_SPEC = "+dx=2.539703 +dy=2.539703 +w=2345 +h=1597 +lat_ll=19.229 +lon_ll=-126.2766" | ||
|
||
EXAMPLE_PIXELS = [ | ||
(0, 0), | ||
(0, 1), | ||
(2000, 1500) | ||
] | ||
EXAMPLE_PROJECTIONS = [ | ||
(-126.2766, 19.229), | ||
(-126.27660549554, 19.229022224607), | ||
(-126.23803580508, 19.272773488997) | ||
] | ||
EXAMPLE_GEOS = [ | ||
(-3275807.350733, -260554.63043505), | ||
(-3275807.350733, -260552.09073205), | ||
(-3270727.944733, -256745.07593505) | ||
] | ||
|
||
|
||
# utility to roughly compare tuples of floats with less floating point precision | ||
def approx_tuple(values: Tuple[float, float]) -> Tuple: | ||
return (approx(values[0]), approx(values[1])) | ||
|
||
|
||
# fixtures | ||
@fixture | ||
def grid_proj() -> GridProj: | ||
return GridProj.from_proj_grid_spec(EXAMPLE_PROJ_SPEC, EXAMPLE_GRID_SPEC) | ||
|
||
|
||
# test class methods | ||
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._h == 1597 | ||
assert grid_proj._w == 2345 | ||
assert grid_proj._dx == approx(2.539703) | ||
assert grid_proj._dy == grid_proj._dx | ||
|
||
t = grid_proj._trans | ||
assert t.source_crs is not None and t.source_crs.type_name == 'Geographic 2D CRS' | ||
assert t.target_crs is not None and t.target_crs.type_name == 'Projected CRS' | ||
|
||
|
||
def test_from_proj_grid_spec_with_offset(): | ||
proj_with_offset = GridProj.from_proj_grid_spec(PROJ_SPEC_WITH_OFFSET, EXAMPLE_GRID_SPEC) | ||
assert proj_with_offset._x_offset == approx(-3272536.1979) | ||
assert proj_with_offset._y_offset == approx(-260290.8369) | ||
|
||
|
||
# transformation methods testing | ||
def test_map_proj_to_pixel(grid_proj: GridProj): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should include testing rounding, once you add rounding option to map_proj_to_pixel(). |
||
for index, proj in enumerate(EXAMPLE_PROJECTIONS): | ||
pixel_x, pixel_y = grid_proj.map_proj_to_pixel(*proj) | ||
assert (pixel_x, pixel_y) == EXAMPLE_PIXELS[index] | ||
|
||
|
||
def test_map_proj_to_geo(grid_proj: GridProj): | ||
for index, proj in enumerate(EXAMPLE_PROJECTIONS): | ||
geo_x, geo_y = grid_proj.map_proj_to_geo(*proj) | ||
assert (geo_x, geo_y) == approx_tuple(EXAMPLE_GEOS[index]) | ||
|
||
|
||
def test_map_pixel_to_geo(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]) | ||
|
||
|
||
def test_map_pixel_to_proj(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]) | ||
|
||
|
||
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_geo_to_pixel(grid_proj: GridProj): | ||
for index, geo in enumerate(EXAMPLE_GEOS): | ||
pixel_x, pixel_y = grid_proj.map_geo_to_pixel(*geo) | ||
assert (pixel_x, pixel_y) == EXAMPLE_PIXELS[index] | ||
|
||
|
||
def test_compound_tranformations_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) | ||
|
||
# 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 | ||
|
||
# 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) | ||
assert (pixel_x, pixel_y) == 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) | ||
|
||
# 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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I found two real comments and a couple nice to haves.