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

Fix handling wkt strings and add pyproj.CRS support #1139

Merged
merged 18 commits into from
Sep 26, 2023
Merged

Conversation

ahuang11
Copy link
Collaborator

@ahuang11 ahuang11 commented Sep 14, 2023

Closes #1137

The issue was that rioxarray gets loaded after running the second cell, which triggered the first if clause, but we didn't handle WKT strings properly; this PR fixes that.

    def _process_crs(self, data, crs):
        """Given crs as proj4 string, data.attr, or cartopy.crs return cartopy.crs
        """
        if hasattr(data, 'rio') and data.rio.crs is not None:
            # if data is a rioxarray
            _crs = data.rio.crs.to_wkt()
            print("RIO")
        else:
            # get the proj string: either the value of data.attrs[crs] or crs itself
            _crs = getattr(data, 'attrs', {}).get(crs or 'crs', crs)
            print("NOT RIO")

Supersedes contaminated #1138

Note the notebook cell execution number
image

@hoxbro
Copy link
Member

hoxbro commented Sep 15, 2023

Thank you for finding the cause of the confusing behavior. Your changes stops the error, but I still find the behavior not to be ideal, e.g.,

import cartopy.crs as ccrs
import rioxarray
import xarray as xr
from hvplot.util import process_crs

crs = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
ds = xr.Dataset(attrs={"crs": crs})

a1 = process_crs(ds.rio.crs.to_proj4())
a2 = process_crs(ds.attrs["crs"])

type(a1), type(a2)

Outputs: (cartopy.crs.PlateCarree, cartopy.crs.RotatedPole)

After looking into it more, it seems these lines are the problem:

hvplot/hvplot/util.py

Lines 130 to 131 in c049c7b

if proj_is_latlong(proj):
return ccrs.PlateCarree()

Looking at this: ds.rio.crs.to_wkt() == crs.to_wkt() will return True, and if I comment out those lines I get type(a1), type(a2) to be (cartopy.crs.RotatedPole, cartopy.crs.RotatedPole)

@ahuang11
Copy link
Collaborator Author

ahuang11 commented Sep 15, 2023

Good catch! I believe if proj_is_latlong is making a wrong assumption that only PlateCarree is lat lon; rotated pole is also in latlong.

    A rotated latitude/longitude projected coordinate system
    with cylindrical topology and projected distance.

Also, I think to_wkt is preferred over to_proj4
https://proj.org/en/9.3/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems

hvplot/util.py Outdated Show resolved Hide resolved
hvplot/util.py Show resolved Hide resolved
hvplot/util.py Outdated Show resolved Hide resolved
hvplot/util.py Outdated Show resolved Hide resolved
hvplot/util.py Outdated Show resolved Hide resolved
@ahuang11 ahuang11 changed the title Fix handling wkt strings Fix handling wkt strings and add pyproj.CRS support Sep 21, 2023
hvplot/util.py Outdated Show resolved Hide resolved
Copy link
Member

@hoxbro hoxbro left a comment

Choose a reason for hiding this comment

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

It looks good to me.

One small thing is we have identical functions in GeoViews. Can you update them and test that we don't break anything in that CI, too, with these changes?

I think it is out-of-scope for this PR, but we are getting some warnings related to these tests. (Already there before this PR)

hvplot/tests/testgeo.py::TestProjections::test_plot_with_crs_as_nonexistent_attr_str
hvplot/tests/testutil.py::TestGeoUtil::test_proj_to_cartopy
hvplot/tests/testutil.py::test_check_crs
hvplot/tests/testutil.py::test_check_crs
hvplot/tests/testutil.py::test_process_crs[+init=epsg:26911]
hvplot/tests/testutil.py::test_process_crs_pyproj_proj
hvplot/tests/testutil.py::test_process_crs_platecarree[+init=epsg:4326]
hvplot/tests/testutil.py::test_process_crs_mercator[+init=epsg:3857]
  /Users/runner/miniconda3/envs/test-environment/lib/python3.11/site-packages/pyproj/crs/crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
    in_crs_string = _prepare_from_proj_string(in_crs_string)

'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]',
# Created with pyproj.CRS("EPSG:3857").to_wkt()
Copy link
Member

Choose a reason for hiding this comment

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

Could be nice to still have this comment.

hvplot/tests/testutil.py Outdated Show resolved Hide resolved
Co-authored-by: Simon Høxbro Hansen <simon.hansen@me.com>
hvplot/tests/testutil.py Outdated Show resolved Hide resolved
@ahuang11 ahuang11 merged commit e4a65f0 into main Sep 26, 2023
7 checks passed
@ahuang11 ahuang11 deleted the fix_rio_crs branch September 26, 2023 20:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Weird behavior with crs vectorfield
3 participants