-
Notifications
You must be signed in to change notification settings - Fork 389
Added winkel tripel projection class #2442
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
base: main
Are you sure you want to change the base?
Conversation
|
|
greglucas
left a comment
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.
Nice, this looks like a great start.
Regarding the tests: Probably add a test similar to one of the other crs's here containing some transform points and initialization expectations: https://github.com/SciTools/cartopy/tree/main/lib/cartopy/tests/crs
Add the projection to the parameterization list here and add a comparison image.
cartopy/lib/cartopy/tests/mpl/test_mpl_integration.py
Lines 161 to 200 in d57fa13
| @pytest.mark.parametrize('proj', [ | |
| ccrs.Aitoff, | |
| ccrs.EckertI, | |
| ccrs.EckertII, | |
| ccrs.EckertIII, | |
| ccrs.EckertIV, | |
| ccrs.EckertV, | |
| ccrs.EckertVI, | |
| ccrs.EqualEarth, | |
| ccrs.Gnomonic, | |
| ccrs.Hammer, | |
| pytest.param((ccrs.InterruptedGoodeHomolosine, dict(emphasis='land')), | |
| id='InterruptedGoodeHomolosine'), | |
| ccrs.LambertCylindrical, | |
| pytest.param((ccrs.Mercator, dict(min_latitude=-85, max_latitude=85)), | |
| id='Mercator'), | |
| ccrs.Miller, | |
| ccrs.Mollweide, | |
| ccrs.NorthPolarStereo, | |
| ccrs.Orthographic, | |
| pytest.param((ccrs.OSGB, dict(approx=True)), id='OSGB'), | |
| ccrs.PlateCarree, | |
| ccrs.Robinson, | |
| pytest.param((ccrs.RotatedPole, | |
| dict(pole_latitude=45, pole_longitude=180)), | |
| id='RotatedPole'), | |
| ccrs.Stereographic, | |
| ccrs.SouthPolarStereo, | |
| pytest.param((ccrs.TransverseMercator, dict(approx=True)), | |
| id='TransverseMercator'), | |
| pytest.param( | |
| (ccrs.ObliqueMercator, dict(azimuth=0.)), id='ObliqueMercator_default' | |
| ), | |
| pytest.param( | |
| (ccrs.ObliqueMercator, dict(azimuth=90., central_latitude=-22)), | |
| id='ObliqueMercator_rotated', | |
| ), | |
| ]) | |
| @pytest.mark.mpl_image_compare(style='mpl20') | |
| def test_global_map(proj): |
A PR you can look at for updating one of the other projections is #1561
lib/cartopy/crs.py
Outdated
| globe = globe or Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS) | ||
| proj4_params = [('proj', 'wintri'), | ||
| ('lon_0', central_longitude), | ||
| ('lat_0', 0.0)] |
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.
Do we need to pass in the lat_0 parameter even if 0 is the default, and does that even work, it looks like it is lat_1 in the documentation: https://proj.org/en/9.4/operations/projections/wintri.html
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.
lat_0 is required and lat_1 doesn't work. Maybe this is a typo in proj? I'm also not very familiar with that library so maybe I'm hacking a parameter.
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 don't think lat_0 is required when I tried this out locally. It looks like it may actually be ignored within proj. lat_1 does work for me and is the first standard parallel. For now, I would just remove the lat_0 from this list. In addition, all of our other _WarpedRectangularProjection subclasses also accept the false_easting and false_northing, so could you pipe those through the initializer as well for consistency with the other classes?
Analysis:
All of these produce a pyproj CRS, but there are differences in the WKT which I think the lat_0 indicates isn't being forwarded because it has a "REMARK" at the end. You can also look at the first standard parallel and the first and second are identical, whereas the third (with lat_1) does have the expected 5.
import pyproj
pyproj.CRS.from_proj4("proj=+wintri").to_wkt()
pyproj.CRS.from_proj4("proj=+wintri +lat_0=5").to_wkt()
pyproj.CRS.from_proj4("proj=+wintri +lat_1=5").to_wkt()'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Winkel Tripel"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Winkel Tripel"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]],REMARK["PROJ CRS string: +proj=wintri +lat_0=5"]]'
'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Winkel Tripel"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Latitude of 1st standard parallel",5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
| return result | ||
|
|
||
|
|
||
| class WinkelTripel(_WarpedRectangularProjection): |
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.
Do we need to do any of the nan handling here like we are doing in the Robinson projection? i.e. should this subclass Robinson instead if it is similar to that?
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 was able to generate maps containing maps with no problem. Are there other tests for nans that I would need to do?
Co-authored-by: Greg Lucas <[email protected]>
|
ping @AgentOxygen did you want to try and get this in still? It looks like it needs tests and a possible update to the keyword argument. |
|
Apologies for the delay, I just finished up my master's thesis! I've added the test and an image for comparison. I noticed several of the tests failed, though I'm not sure how to proceed. What are the next steps I should take to get this integrated. |
lib/cartopy/crs.py
Outdated
| globe = globe or Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS) | ||
| proj4_params = [('proj', 'wintri'), | ||
| ('lon_0', central_longitude), | ||
| ('lat_0', 0.0)] |
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 don't think lat_0 is required when I tried this out locally. It looks like it may actually be ignored within proj. lat_1 does work for me and is the first standard parallel. For now, I would just remove the lat_0 from this list. In addition, all of our other _WarpedRectangularProjection subclasses also accept the false_easting and false_northing, so could you pipe those through the initializer as well for consistency with the other classes?
Analysis:
All of these produce a pyproj CRS, but there are differences in the WKT which I think the lat_0 indicates isn't being forwarded because it has a "REMARK" at the end. You can also look at the first standard parallel and the first and second are identical, whereas the third (with lat_1) does have the expected 5.
import pyproj
pyproj.CRS.from_proj4("proj=+wintri").to_wkt()
pyproj.CRS.from_proj4("proj=+wintri +lat_0=5").to_wkt()
pyproj.CRS.from_proj4("proj=+wintri +lat_1=5").to_wkt()'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Winkel Tripel"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Winkel Tripel"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]],REMARK["PROJ CRS string: +proj=wintri +lat_0=5"]]'
'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Winkel Tripel"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Latitude of 1st standard parallel",5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
lib/cartopy/crs.py
Outdated
| ('lat_0', 0.0)] | ||
|
|
||
| super().__init__(proj4_params, | ||
| central_longitude, |
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.
Should be indented the same amount as the parentheses, see the other classes in this file for an example.
🎉 Congratulations, @AgentOxygen! |
Removing unused central latitude
Fixing formatting of super() call
|
@AgentOxygen, it looks like the image comparison is failing now. I pushed the formatting fix, but otherwise this looks good to go to me once you can update the images to the expected content. |
Rationale
Winkel-Tripel is a popular projection that would be used by multiple users (my group at UT and folks at NCAR). It is also the primary projection for the National Geographic Society.
Implications
I have created a new CRS projection object without modifying any code by utilizing existing class objects and the pyproj backend.
Checklist
This is my first pull request! I am open to input how on to properly contribute.