I am trying to plot a topographic raster in Cartopy. I have downloaded some sample GeoTIFF data from this database: https://zenodo.org/record/3940482. I then import the data and metadata using the Python GDAL library:
from osgeo import gdal
data_object = gdal.Open(path_geotiff)
data_array = data_object.ReadAsArray()
data_transform = data_object.GetGeoTransform()
proj_wkt = data_object.GetProjection()
I need to get the projection as a Cartopy object, which I achieve using the PyProj library for an intermediate step:
from pyproj.crs import CRS
proj_crs = CRS.from_wkt(proj_wkt)
import cartopy.crs as ccrs
proj_cartopy = ccrs.Projection(proj_crs)
This seems to work, and indicates that the projection is one of those supported by Cartopy, Albers Equal Area:
>>> print(proj_cartopy)
PROJCRS["Albers",BASEGEOGCRS["NAD83",DATUM["North American Datum 1983",ELLIPSOID["GRS 1980",6378137,298.257222101004,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4269]],CONVERSION["unnamed",METHOD["Albers Equal Area",ID["EPSG",9822]],PARAMETER["Latitude of false origin",23,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",-106,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",29.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",45.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]
I then try to create a GeoAxes object:
import matplotlib.pyplot as plt
subplot_kw = dict(projection = proj_cartopy)
fig, ax = plt.subplots(subplot_kw = subplot_kw)
However, I receive a NotImplemented error:
File "/Users/my_username/opt/anaconda3/envs/carto/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py", line 1598, in _boundary path, = cpatch.geos_to_path(self.projection.boundary)
File "/Users/my_username/opt/anaconda3/envs/carto/lib/python3.9/site-packages/cartopy/crs.py", line 679, in boundary
raise NotImplementedError
Please help me to understand this error (I know I could manually enter the projection, for example using the EPSG code, but I prefer to keep my current automatic method if possible).
CodePudding user response：
From what I understand, this is something that is not fully supported yet.