Skip to content

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. #2510

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

Open
Curallin opened this issue Mar 11, 2025 · 4 comments

Comments

@Curallin
Copy link

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)

Image

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)

Image

@greglucas
Copy link
Contributor

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
Copy link
Author

@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)

Image

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

Image

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

Image

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

@greglucas
Copy link
Contributor

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: #923

@Curallin
Copy link
Author

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()

Image

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

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

No branches or pull requests

2 participants