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
55 changes: 51 additions & 4 deletions hvplot/tests/testutil.py
Original file line number Diff line number Diff line change
Expand Up @@ -217,6 +217,14 @@ def test_proj_to_cartopy(self):

assert isinstance(crs, self.ccrs.CRS)

def test_proj_to_cartopy_wkt_string(self):
from ..util import proj_to_cartopy
crs = proj_to_cartopy('GEOGCRS["unnamed",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],DERIVINGCONVERSION["unknown",METHOD["PROJ ob_tran o_proj=latlon"],PARAMETER["o_lon_p",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lat_p",37.5,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["lon_0",357.5,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]') # noqa: E501

assert isinstance(crs, self.ccrs.RotatedPole)
assert crs.proj4_params["lon_0"] == 357.5
assert crs.proj4_params["o_lat_p"] == 37.5


class TestDynamicArgs(TestCase):

Expand Down Expand Up @@ -283,20 +291,59 @@ def test_check_crs():

@pytest.mark.parametrize("input", [
"+init=epsg:26911",
])
def test_process_crs(input):
pytest.importorskip("pyproj")
ccrs = pytest.importorskip("cartopy.crs")
crs = process_crs(input)
assert isinstance(crs, ccrs.CRS)


def test_process_crs_pyproj_crs():
pyproj = pytest.importorskip("pyproj")
ccrs = pytest.importorskip("cartopy.crs")
crs = process_crs(pyproj.CRS.from_epsg(4326))
assert isinstance(crs, ccrs.PlateCarree)


def test_process_crs_pyproj_proj():
pyproj = pytest.importorskip("pyproj")
ccrs = pytest.importorskip("cartopy.crs")
crs = process_crs(pyproj.Proj(init='epsg:4326'))
assert isinstance(crs, ccrs.PlateCarree)


@pytest.mark.parametrize("input", [
"4326",
4326,
"epsg:4326",
"EPSG: 4326",
"+init=epsg:4326",
'GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]',
ahuang11 marked this conversation as resolved.
Show resolved Hide resolved
], ids=lambda x: str(x)[:20])
def test_process_crs_platecarree(input):
pytest.importorskip("pyproj")
ccrs = pytest.importorskip("cartopy.crs")
crs = process_crs(input)
assert isinstance(crs, ccrs.PlateCarree)


@pytest.mark.parametrize("input", [
"3857",
3857,
"epsg:3857",
"EPSG: 3857",
"+init=epsg:3857",
'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.

'PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["Popular Visualisation Pseudo-Mercator",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Web mapping and visualisation."],AREA["World between 85.06°S and 85.06°N."],BBOX[-85.06,-180,85.06,180]],ID["EPSG",3857]]'
# Created with pyproj.CRS("EPSG:4326").to_wkt()
ahuang11 marked this conversation as resolved.
Show resolved Hide resolved
'PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["Popular Visualisation Pseudo-Mercator",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Web mapping and visualisation."],AREA["World between 85.06°S and 85.06°N."],BBOX[-85.06,-180,85.06,180]],ID["EPSG",3857]]',
ahuang11 marked this conversation as resolved.
Show resolved Hide resolved
], ids=lambda x: str(x)[:20])
def test_process_crs(input):
def test_process_crs_mercator(input):
pytest.importorskip("pyproj")
ccrs = pytest.importorskip("cartopy.crs")
crs = process_crs(input)
assert isinstance(crs, ccrs.Mercator)

assert isinstance(crs, ccrs.CRS)

def test_process_crs_rasterio():
pytest.importorskip("pyproj")
Expand Down
61 changes: 34 additions & 27 deletions hvplot/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ def check_crs(crs):

Returns
-------
A valid crs if possible, otherwise None
A valid crs if possible, otherwise None.
"""
import pyproj

Expand All @@ -79,8 +79,11 @@ class Dummy:
out = pyproj.Proj(crs.to_wkt(), preserve_units=True)
elif isinstance(crs, dict) or isinstance(crs, str):
if isinstance(crs, str):
# quick fix for https://github.com/pyproj4/pyproj/issues/345
crs = crs.replace(' ', '').replace('+', ' +')
try:
crs = pyproj.CRS.from_wkt(crs)
except RuntimeError:
# quick fix for https://github.com/pyproj4/pyproj/issues/345
crs = crs.replace(' ', '').replace('+', ' +')
try:
out = pyproj.Proj(crs, preserve_units=True)
except RuntimeError:
Expand Down Expand Up @@ -124,10 +127,10 @@ def proj_to_cartopy(proj):
except ImportError:
has_gdal = False

proj = check_crs(proj)

if proj_is_latlong(proj):
return ccrs.PlateCarree()
input_proj = proj
proj = check_crs(input_proj)
if proj is None:
raise ValueError(f"Invalid proj projection {input_proj!r}")

srs = proj.srs
if has_gdal:
Expand Down Expand Up @@ -167,19 +170,23 @@ def proj_to_cartopy(proj):
except:
pass
if k == 'proj':
if v == 'tmerc':
if v == "longlat":
cl = ccrs.PlateCarree
elif v == 'tmerc':
cl = ccrs.TransverseMercator
kw_proj['approx'] = True
if v == 'lcc':
elif v == 'lcc':
cl = ccrs.LambertConformal
if v == 'merc':
elif v == 'merc':
cl = ccrs.Mercator
if v == 'utm':
elif v == 'utm':
cl = ccrs.UTM
if v == 'stere':
elif v == 'stere':
cl = ccrs.Stereographic
if v == 'ob_tran':
elif v == 'ob_tran':
cl = ccrs.RotatedPole
else:
raise NotImplementedError(f'Unknown projection {v}')
if k in km_proj:
if k == 'zone':
v = int(v)
Expand Down Expand Up @@ -228,7 +235,9 @@ def process_crs(crs):
1. EPSG codes: Defined as string of the form "EPSG: {code}" or an integer
2. proj.4 string: Defined as string of the form "{proj.4 string}"
3. cartopy.crs.CRS instance
4. None defaults to crs.PlateCaree
3. pyproj.Proj or pyproj.CRS instance
4. WKT string: Defined as string of the form "{WKT string}"
5. None defaults to crs.PlateCaree
"""
missing = []
try:
Expand All @@ -248,29 +257,27 @@ def process_crs(crs):

if crs is None:
return ccrs.PlateCarree()
elif isinstance(crs, ccrs.CRS):
return crs
elif isinstance(crs, pyproj.CRS):
crs = crs.to_wkt()

errors = []
if isinstance(crs, str) and crs.lower().startswith('epsg'):
if isinstance(crs, (str, int)): # epsg codes
try:
crs = crs[5:].lstrip().rstrip()
return ccrs.epsg(crs)
crs = pyproj.CRS.from_epsg(crs).to_wkt()
ahuang11 marked this conversation as resolved.
Show resolved Hide resolved
except Exception as e:
errors.append(e)
if isinstance(crs, int):
try:
return ccrs.epsg(crs)
except Exception as e:
crs = str(crs)
errors.append(e)
if isinstance(crs, (str, pyproj.Proj)):
if isinstance(crs, (str, pyproj.Proj)): # proj4/wkt strings
try:
return proj_to_cartopy(crs)
except Exception as e:
errors.append(e)
if isinstance(crs, ccrs.CRS):
return crs

raise ValueError("Projection must be defined as a EPSG code, proj4 string, cartopy CRS or pyproj.Proj.") from Exception(*errors)
raise ValueError(
"Projection must be defined as a EPSG code, proj4 string, "
"WKT string, cartopy CRS, pyproj.Proj, or pyproj.CRS."
) from Exception(*errors)


def is_list_like(obj):
Expand Down
Loading