Skip to content
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

WIP: Updating to PROJ 6+ proj.h C API #1252

Closed
wants to merge 2 commits into from
Closed

Conversation

snowman2
Copy link
Contributor

@snowman2 snowman2 commented Dec 31, 2018

This likely has bugs and whatnot at this stage, but it gets the concepts of changes needed out there for discussion.

Rationale

#1140

@SciTools-assistant SciTools-assistant added the Blocked: CLA needed See https://scitools.org.uk. Submit the form at: https://scitools.org.uk/cla/v4/form label Dec 31, 2018
@snowman2 snowman2 force-pushed the proj6 branch 3 times, most recently from 28fba75 to 6dbff9a Compare January 2, 2019 02:24
@snowman2
Copy link
Contributor Author

snowman2 commented Jan 2, 2019

Current issues in CRS tests is the error:

proj_create_operations: Source and target ellipsoid do not belong to the same celestial body

@snowman2
Copy link
Contributor Author

Also, this will mean that python 2.7 won't work on windows due to c++11 requirement for Proj 6.0.0

@dopplershift
Copy link
Contributor

I think we were already close to dropping 2.7, so this is probably as good a reason as any. Might want to have one more release before we merge this...but then we'd have to be careful about what proj C library is used to build. Ick. @pelson ?

@snowman2
Copy link
Contributor Author

snowman2 commented Apr 2, 2019

I wouldn't recommend merging this just yet. It was created before 6.0.0 was released in preparation for the release and the API has changed since then.

@SciTools-assistant SciTools-assistant removed the Blocked: CLA needed See https://scitools.org.uk. Submit the form at: https://scitools.org.uk/cla/v4/form label Apr 10, 2019
setup.py Show resolved Hide resolved
setup.py Outdated Show resolved Hide resolved
@snowman2 snowman2 force-pushed the proj6 branch 2 times, most recently from 923db9f to a333bba Compare April 10, 2019 01:19
@snowman2
Copy link
Contributor Author

Well, at least the errors I mentioned earlier are showing in the build now :)

@snowman2
Copy link
Contributor Author

I see a lot of failures due to the Source and target ellipsoid do not belong to the same celestial body error.

More information about this issue:

Essentially, the tests are using invalid transformations in multiple places if I understand correctly.

@@ -1188,7 +1188,7 @@ def __init__(self, central_longitude=-96.0, central_latitude=39.0,
lons[1:-1] = np.linspace(central_longitude - 180 + 0.001,
central_longitude + 180 - 0.001, n)

points = self.transform_points(PlateCarree(), lons, lats)
points = self.transform_points(PlateCarree(globe=self.globe), lons, lats)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E501 line too long (81 > 79 characters)

@@ -1266,7 +1266,7 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0,
lon = central_longitude + 180
sign = np.sign(central_latitude) or 1
lat = -central_latitude + sign * 0.01
x, max_y = self.transform_point(lon, lat, PlateCarree())
x, max_y = self.transform_point(lon, lat, PlateCarree(globe=self.globe))

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E501 line too long (80 > 79 characters)

@@ -190,8 +190,8 @@ def test_globe(self):
rugby_moll = ccrs.Mollweide(globe=rugby_globe)
footy_moll = ccrs.Mollweide(globe=footy_globe)

rugby_pt = rugby_moll.transform_point(10, 10, ccrs.Geodetic())
footy_pt = footy_moll.transform_point(10, 10, ccrs.Geodetic())
rugby_pt = rugby_moll.transform_point(10, 10, ccrs.Geodetic(globe=rugby_globe))

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E501 line too long (87 > 79 characters)

rugby_pt = rugby_moll.transform_point(10, 10, ccrs.Geodetic())
footy_pt = footy_moll.transform_point(10, 10, ccrs.Geodetic())
rugby_pt = rugby_moll.transform_point(10, 10, ccrs.Geodetic(globe=rugby_globe))
footy_pt = footy_moll.transform_point(10, 10, ccrs.Geodetic(globe=footy_globe))

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E501 line too long (87 > 79 characters)

@@ -504,7 +504,7 @@ def format_coord(self, x, y):
A string formatted for the Matplotlib GUI status bar.

"""
lon, lat = ccrs.Geodetic().transform_point(x, y, self.projection)
lon, lat = ccrs.Geodetic(self.projection.globe).transform_point(x, y, self.projection)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E501 line too long (94 > 79 characters)

@@ -583,7 +583,7 @@ def tissot(self, rad_km=500, lons=None, lats=None, n_samples=80, **kwargs):
n_samples=n_samples)
geoms.append(sgeom.Polygon(circle))

feature = cartopy.feature.ShapelyFeature(geoms, ccrs.Geodetic(),
feature = cartopy.feature.ShapelyFeature(geoms, ccrs.Geodetic(self.projection.globe),

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E501 line too long (93 > 79 characters)

@snowman2
Copy link
Contributor Author

Question for @pelson @dopplershift @QuLogic
Do y'all know the reasoning behind the logic for the default globe for PlateCaree?

globe = Globe(semimajor_axis=math.degrees(1))

>>> from cartopy.crs import PlateCarree
>>> pc = PlateCarree()
>>> pc.proj4_init
'+proj=eqc +ellps=WGS84 +a=57.29577951308232 +lon_0=0.0 +no_defs'

The semimajor_axis varies quite a bit from any other value I have seen.

Also, from what I can tell, it seems that meters are the units expected for the semi major axis:

>>> from pyproj import CRS
>>> crs = CRS('+proj=eqc +ellps=WGS84 +a=57.29577951308232 +lon_0=0.0 +no_defs')
>>> crs.ellipsoid
ELLIPSOID["unknown",57.2957795130823,298.257223563,
    LENGTHUNIT["metre",1,
        ID["EPSG",9001]]]
>>> crs.ellipsoid.semi_major_metre
57.29577951308232

@snowman2
Copy link
Contributor Author

Another data point on the plate carree ellipsoid base on https://proj.org/operations/projections/eqc.html:

>>> from pyproj import CRS
>>> cc = CRS.from_epsg(32622)
>>> cc.ellipsoid
ELLIPSOID["WGS 84",6378137,298.257223563,
    LENGTHUNIT["metre",1],
    ID["EPSG",7030]]
>>> cc.ellipsoid.semi_major_metre
6378137.0
>>> ceqc = CRS("+proj=eqc +lat_ts=0")
>>> ceqc.ellipsoid
ELLIPSOID["WGS 84",6378137,298.257223563,
    LENGTHUNIT["metre",1],
    ID["EPSG",7030]]
>>> ceqc.ellipsoid.semi_major_metre
6378137.0

It seems 6378137.0 could be a good default?

@dopplershift
Copy link
Contributor

@snowman2 That particular choice of ellipsoid as the default for PlateCarree is a hack clever solution to allow passing lats and lons as Cartesian coordinates to PlateCarree, which is really just an equirectangular projection. That choice of Globe also means you get unexpected behavior when changing it (i.e. it no longer takes lons/lats). See #709.

@snowman2
Copy link
Contributor Author

Thanks for the explanation - that make sense. Well, it seems this is no longer allowed in PROJ 6+, so it will need to be updated before upgrading.

Based on #1164, I think it would make sense to use latlon. In pyproj we default to assume all input is degrees unless you toggle the radians flag and do conversions accordingly as needed by PROJ:
https://github.com/pyproj4/pyproj/blob/cecbdf4a516bf031c6e8fab0cf7188148dafc9f4/pyproj/_transformer.pyx#L164-L176

Is this something you would be interested in?

@dopplershift
Copy link
Contributor

@snowman2 I think what we need is:

  1. Make a new lat/lon projection that takes degrees
  2. Add an EquidistantCylindrical projection that works like everything else with regards to globes
  3. Deprecate PlateCarree, or at least its current behavior

@snowman2
Copy link
Contributor Author

snowman2 commented Jul 6, 2019

@dopplershift, I think something just clicked in my brain. Is the goal to create something that converts from the "WGS84" or "EPGS:4326" to a global projected coordinate system? Or to convert from the geodetic CRS to the projected CRS. Either way, I think a new base class is needed. This class would be similar to the pyproj.Transformer, however it would default to have the from_crs always be either "WGS84" or the geodetic CRS based on the projected CRS depending on what you want (See: http://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter). This would enable the globe parameter to be supported as well open up the option for users to use additional Transformer classes that convert from latlon to a projected coordinate system.

@snowman2
Copy link
Contributor Author

snowman2 commented Jul 6, 2019

An idea of different global projected coordinate systems could include:
EPSG:3857 (Pseudo web mercator)
EPSG:3395 "WGS 84 / World Mercator"

Ideas from: https://lists.osgeo.org/pipermail/proj/2019-June/008653.html

@snowman2
Copy link
Contributor Author

Currently upgrading to use the proj.h API in PROJ 6+ (an soon 7) has many changes/challenges. @dopplershift has expressed interest in potentially using pyproj 2 itself within cartopy to replace the direct need for building against the proj.h API. Is this something that the cartopy project is interested in?

@greglucas
Copy link
Contributor

I'm 👍 on using pyproj2, since you've already built in all of the new Transformer hooks. Let's see how others feel about this too though. Would this remove the need for Cartopy to use pyepsg too? There was a bit of a hiccup recently when epsg.io went down and not depending on an external website would be nice.

Cartopy is currently still supporting some pretty old Python versions, so we would have to drop some older Python versions to make this work it looks like.

I think that all of this should actually be thought about from a broader perspective of what packages should handle what operations in a future big release. It looks like pygeos and Shapely will be merging their GEOS library work together in the (hopefully not too distant) future. There was a discussion about moving the GEOS optimizations out of Cartopy. Maybe this goes in parallel with that update to start moving some of the Cartopy Cython optimization code into other upstream libraries where it makes sense. (i.e. code in trace.pyx to shapely/pygeos, and _crs.pyx into pyproj2)

@snowman2
Copy link
Contributor Author

Let's see how others feel about this too though.

Should I open up a separate issue?

Would this remove the need for Cartopy to use pyepsg too?

Yes, it will remove that.

@dopplershift
Copy link
Contributor

@snowman2 A separate issue might make sense, since we're veering off the topic of this PR. I'll go ahead and open it with my thoughts.

@dopplershift
Copy link
Contributor

See #1477.

@snowman2
Copy link
Contributor Author

Going to close in favor of #1477.

@snowman2 snowman2 closed this Feb 25, 2020
@snowman2 snowman2 changed the title WIP: updating to proj.4 6.0.0 WIP: Updating to PROJ 6+ proj.h C API Jun 23, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants