CONUS Albers CRS works in hvplot.pandas, but not in hvplot.xarray?

hvPlot
1

I’ve got a geotiff that is geolocated correctly in QGIS (so I know the CRS is correct), but when I load it with rioxarray into xarray I’m struggling to get it visualized correctly with Hvplot.

I tried using both the Rasterio CRS and converting to a Cartopy CRS, but hvplot is unhappy with both of them.

This is a reproducible example:

import xarray as xr
import hvplot.xarray
import fsspec

fs = fsspec.filesystem('s3', anon=True, endpoint_url='https://usgs.osn.mghpcc.org')
ds = xr.open_dataset(fs.open('esip/rsignell/testing/1646955_2487285.tiff'), engine='rasterio').squeeze('band', drop=True)
rasterio_crs = ds['band_data'].rio.crs
ds['band_data'].hvplot(rasterize=True, crs=rasterio_crs)

which fails with:

Exception: (CRSError('Invalid projection: ...
ValueError: Projection must be defined as a EPSG code, proj4 string, WKT string, cartopy CRS instance, cartopy CRS class name string, pyproj.Proj, or pyproj.CRS.

which was strange, because I was pretty sure that I’ve used a rasterio CRS before. But I went ahead and converted the Rasterio CRS to a Cartopy CRS using:

from pyproj.crs import CRS as PyprojCRS
import cartopy.crs as ccrs

pyproj_crs = PyprojCRS.from_user_input(rasterio_crs)
cartopy_crs = ccrs.CRS(pyproj_crs.to_wkt())

but when I tried:

ds['band_data'].hvplot(rasterize=True, crs=cartopy_crs)

I got the same error:

ValueError: Projection must be defined as a EPSG code, proj4 string, WKT string, cartopy CRS instance, cartopy CRS class name string, pyproj.Proj, or pyproj.CRS.

The problem doesn’t seem to be with rasterio, since if I use rasterio to reproject it works fine:

image
image1182×659 141 KB

Using da.hvplot(crs=...) should work, right?

What am I missing??

Full reproducible notebook is here.

2

Maybe you can try pyproj?

ValueError: Projection must be defined as a EPSG code, proj4 string, WKT string, cartopy CRS instance, cartopy CRS class name string, pyproj.Proj, or pyproj.CRS.

3

Ah would that it be that easy! :smile:
Nope, same error. :frowning_face:
I’ve updated the reproducible notebook to show that.

The projection used here is actually this EPSG:5070

4

Hmm, the plot thickens. I tried a geotiff image from a different dataset on AWS and the rasterio CRS worked just fine:

import fsspec
import xarray as xr
import hvplot.xarray

fs = fsspec.filesystem('s3', anon=True)
url = 'sentinel-1-global-coherence-earthbigdata/data/tiles/N00E029/N00E029_winter_vv_rmse.tif'
ds = xr.open_dataset(fs.open(url), engine='rasterio').squeeze('band', drop=True)
crs = ds['band_data'].rio.crs
ds['band_data'].hvplot(rasterize=True, crs=crs)

image

So why would this rasterio crs work:

CRS.from_wkt(‘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”]],AXIS[“Latitude”,NORTH],AXIS[“Longitude”,EAST],AUTHORITY[“EPSG”,“4326”]]’)

and not this one:

CRS.from_wkt(‘PROJCS[“Projection = Albers Conical Equal Area”,GEOGCS[“WGS 84”,DATUM[“WGS_1984”,SPHEROID[“WGS 84”,6378137,298.257223563,AUTHORITY[“EPSG”,“7030”]],AUTHORITY[“EPSG”,“6326”]],PRIMEM[“Greenwich”,0],UNIT[“degree”,0.0174532925199433,AUTHORITY[“EPSG”,“9122”]],AUTHORITY[“EPSG”,“4326”]],PROJECTION[“Albers_Conic_Equal_Area”],PARAMETER[“latitude_of_center”,23],PARAMETER[“longitude_of_center”,-96],PARAMETER[“standard_parallel_1”,29.5],PARAMETER[“standard_parallel_2”,45.5],PARAMETER[“false_easting”,0],PARAMETER[“false_northing”,0],UNIT[“metre”,1,AUTHORITY[“EPSG”,“9001”]],AXIS[“Easting”,EAST],AXIS[“Northing”,NORTH]]’)

5

Update. I had the bright idea to just create the CRS from scratch using Cartopy:

import cartopy.crs as ccrs
albers_crs = ccrs.AlbersEqualArea(
    central_longitude=-96,
    central_latitude=23,
    standard_parallels=(29.5, 45.5)
)

and typing albers_crs in the notebook results in:

image
image657×530 51.4 KB

so I was excited, and thought this must then work:

ds['band_data'].hvplot(rasterize=True, crs=albers_crs)

but NO it doesn’t. What the fudge? Give same error as before!

6

So my question now is why doesn’t this simple code work:

image
image1147×843 70.2 KB

Reproducible notebook here
Are there only some cartopy CRS specifications that hvplot/geoviews can handle?

7

Okay, WTF (what the fudge!): if I convert some of the points to a dataframe, it works fine:

import hvplot.pandas
df = ds['band_data'][::100,::100].to_dataframe()
df.hvplot.points(
    x='x', 
    y='y',
    hover_cols=['point_id', 'description'],
    crs=albers_crs,
    tiles='OSM'
)

produces:
image

Why doesn’t hvplot.xarray work?

Reproducible notebook

8

Have you tried upgrading your packages?

hvplot 0.11.3
It seems to be running for me, but not finishing

import fsspec
import xarray as xr
import hvplot.xarray

fs = fsspec.filesystem('s3', anon=True)
url = 'sentinel-1-global-coherence-earthbigdata/data/tiles/N00E029/N00E029_winter_vv_rmse.tif'
ds = xr.open_dataset(fs.open(url), engine='rasterio').squeeze('band', drop=True)
crs = ds['band_data'].rio.crs
ds['band_data'].hvplot(rasterize=True, crs=crs)
9

That global coherence notebook is the one that is working.

It’s this one that isn’t working:
conus_albers_issue.ipynb

10

Thanks I think the offending lines are:


        if hasattr(data, 'rio') and data.rio.crs is not None:
            # if data is a rioxarray
            _crs = data.rio.crs.to_wkt()

I’m looking into it now.

11

Okay PR is up; if you can try it out on your datasets, that’d be appreciated!

1 Like