Why can't I set data_crs using ccrs.CRS(), how can I do it if I already know the coordinate system of the data.

[Curallin]

I want to use the data under the utm projection but using the latitude/longitude scale, and for that I wrote the function to visualize the dem. But it seems to work only for data with EPSG=4326 and not the utm projection (EPSG=3857).How can I figure it out? The function is below:

import rioxarray  as rxr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib import rcParams

config = {'font.family': ['Arial','Microsoft Yahei'],'font.size':18}
rcParams.update(config)

def visualize_single_band_data(DataArray,label='Elevation(m)',map_crs=ccrs.PlateCarree(),data_crs=ccrs.PlateCarree(),cmap='terrain'):
    bounds = DataArray.rio.bounds()
    extent = bounds[0],bounds[2],bounds[1],bounds[3]
    nodata = DataArray._FillValue
    DataArray = DataArray.where(DataArray!= nodata)
    fig,ax = plt.subplots(figsize=(12,9),subplot_kw={'projection':map_crs})
    ax.set_extent(extent,data_crs)
    img = ax.imshow(DataArray.squeeze(), origin='upper', extent=extent, transform=data_crs,cmap=cmap)
    ax.gridlines(draw_labels={'top':'x','left':'y'}, linewidth=1, color='gray', alpha=0.5, linestyle='--')
    cbar = plt.colorbar(img,ax=ax,orientation='horizontal',pad=0.05)
    cbar.set_label(label=label)
    plt.show()

It works properly under the following conditions:

visualize_single_band_data(cn_dem)

But when I use data under utm projection, it causes errors.

visualize_single_band_data(cn_dem_utm,data_crs=cn_dem_utm.rio.crs)

[greglucas]

Does it work if you use a Projection subclass like the message suggests? Also it helps if you copy/paste raw text rather than screenshots of the issues and make a minimal reproducer that we can run locally without all of the extras.

data_crs=ccrs.Projection(cn_dem_utm.rio.crs), or if you can make the epsg directly I think that might work ccrs.epsg(3857)

[Curallin]

@greglucas The former answer may not work,but the latter one does work in utm projection.However, it does nor work for epsg=4326. Next, I will show you in details. The former solution(data_crs=ccrs.Projection(cn_dem_utm.rio.crs)):

import rioxarray  as rxr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib import rcParams

config = {'font.family': ['Arial','Microsoft Yahei'],'font.size':18}
rcParams.update(config)

def visualize_single_band_data(DataArray,label='Elevation(m)',map_crs=ccrs.PlateCarree(),cmap='terrain'):
    data_crs=ccrs.Projection(DataArray.rio.crs)
    bounds = DataArray.rio.bounds()
    extent = bounds[0],bounds[2],bounds[1],bounds[3]
    nodata = DataArray._FillValue 
    DataArray = DataArray.where(DataArray!= nodata)
    fig,ax = plt.subplots(figsize=(12,9),subplot_kw={'projection':map_crs})
    ax.set_extent(extent,data_crs)
    img = ax.imshow(DataArray.squeeze(), origin='upper', extent=extent, transform=data_crs,cmap=cmap)
    ax.gridlines(draw_labels={'top':'x','left':'y'}, linewidth=1, color='gray', alpha=0.5, linestyle='--')
    cbar = plt.colorbar(img,ax=ax,orientation='horizontal',pad=0.05)
    cbar.set_label(label=label)
    plt.show()
visualize_single_band_data(cn_dem)

The latter one:

def visualize_single_band_data(DataArray,label='Elevation(m)',map_crs=ccrs.PlateCarree(),data_crs=ccrs.PlateCarree(),cmap='terrain'):
    epsg = DataArray.rio.crs.to_epsg()
    data_crs = ccrs.epsg(epsg)
    bounds = DataArray.rio.bounds()
    extent = bounds[0],bounds[2],bounds[1],bounds[3]
    nodata = DataArray._FillValue 
    DataArray = DataArray.where(DataArray!= nodata)
    fig,ax = plt.subplots(figsize=(12,9),subplot_kw={'projection':map_crs})
    ax.set_extent(extent,data_crs)
    img = ax.imshow(DataArray.squeeze(), origin='upper', extent=extent, transform=data_crs,cmap=cmap)
    ax.gridlines(draw_labels={'top':'x','left':'y'}, linewidth=1, color='gray', alpha=0.5, linestyle='--')
    cbar = plt.colorbar(img,ax=ax,orientation='horizontal',pad=0.05)
    cbar.set_label(label=label)
    plt.show()
visualize_single_band_data(cn_dem_utm) # same in cn_dem_utm

Well,for utm it does work

However,when I use cn_dem which epsg is 4326, it causes errors.

There seems to be a need for a generalized approach to this problem

[greglucas]

Does this page provide any help for moving between the CRS'? https://pyproj4.github.io/pyproj/stable/crs_compatibility.html#converting-from-rasterio-crs-crs-to-pyproj-crs-crs

I agree that this would be good to add some examples and better interoperability between all of the various packages.

xref previous discussions: https://github.com/SciTools/cartopy/issues/923

[Curallin]

def visualize_single_band_data(DataArray,label='Elevation(m)',map_crs=ccrs.PlateCarree(),cmap='terrain'):
    epsg = cn_dem.rio.crs.to_epsg()
    data_crs = CRS.from_epsg(epsg)
    # data_crs_str = DataArray.rio.crs.to_proj4()
    # data_crs = ccrs.Projection(data_crs_str)
    bounds = DataArray.rio.bounds()
    extent = bounds[0],bounds[2],bounds[1],bounds[3]
    nodata = DataArray._FillValue 
    DataArray = DataArray.where(DataArray!= nodata)
    fig,ax = plt.subplots(figsize=(12,9),subplot_kw={'projection':map_crs})
    ax.set_extent(extent,data_crs)
    img = ax.imshow(DataArray.squeeze(), origin='upper', extent=extent, transform=data_crs,cmap=cmap)
    ax.gridlines(draw_labels={'top':'x','left':'y'}, linewidth=1, color='gray', alpha=0.5, linestyle='--')
    cbar = plt.colorbar(img,ax=ax,orientation='horizontal',pad=0.05)
    cbar.set_label(label=label)
    plt.show()

But using proj4 does work for both conditions,but I get a FutureWarning:

d:\PyEnv\envs\dp\Lib\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)

I don't know if I'll be able to use it after the version update