2

I read a shapefile containing geometries that I want to plot using Cartopy. My challenge is finding the correct CRS projection to use for this. The .prj file contains this information:

PROJCS["NAD_1983_Albers",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-154.0],PARAMETER["Standard_Parallel_1",55.0],PARAMETER["Standard_Parallel_2",65.0],PARAMETER["Latitude_Of_Origin",50.0],UNIT["Meter",1.0]]

When I plot the shapefile in Google Earth it looks correct, Google Earth with correct projection

but when I plot using Cartopy it is off Cartopy with wrong projection

I obviously use the incorrect projection (crs.AlbersEqualArea assumed from PROJECTION["Albers"]), but I am not sure how to identify the correct one from the .prj file.

How can I extract the necessary information and plot the shapefile using Cartopy without having to guess the projection? I would welcome any suggestions.

The shapefile is found here if anyone would take a look and the code for making the map is here:

import geopandas as gdp

def create_colors(N):
   color = iter(cm.tab20b(np.linspace(0,1,N)))
   return [next(color) for c in range(N)]

ax.coastlines(resolution="10m", linewidth=0.6, color="black", alpha=0.8, zorder=4)
ax.add_feature(cpf.BORDERS, linestyle=':',alpha=0.4)
ax.add_feature(cpf.LAND, color="orange")
extent = [-182, -100, 40, 65]

projection = ccrs.AlbersEqualArea(central_longitude=-155,
                        central_latitude=50,
                        standard_parallels=(0, 80))

ax.set_extent(extent, crs=ccrs.PlateCarree())

shdf = gpd.read_file(lme_file)

LMES = ["630","640","650","620","610"]
counter = 0
for LME_NUMBER in shdf['NMFS_AREA']:

    shdf_sel = shdf[shdf['NMFS_AREA']==LME_NUMBER]
    ax.add_geometries(shdf_sel['geometry'],
                          projection,
                          facecolor='LightGray',
                          edgecolor='k')

plt.show()

I managed to get the projection correct by manually adding the optional parameters to Albers Equal Area projection in Cartopy.

projection = ccrs.AlbersEqualArea(central_longitude=-154, central_latitude=50, false_easting=0, false_northing=0)

This gives a nice result like this: enter image description here

I still would like to know if there are any way I could have extracted this information automatically from the shapefile without having to manually look into the .prj file, but regardless this works.

1

This turned out to be pretty straightforward when usin GeoPandas. I first converted the shapefile to a standard projection (EPSG:4326):

import geopandas as gpd

shdf = gpd.read_file("shapefile")
shdf = shdf.to_crs("EPSG:4326")

Next, I used the PlateCarreeas my projection in Cartopy. I could then add my shapefile geometries using:

ax = plt.figure(figsize=(16,10)).gca(projection=ccrs.PlateCarree())
ax.coastlines(resolution=res, linewidth=0.6, color="black", alpha=0.8, zorder=4)
ax.add_feature(cpf.BORDERS, linestyle=':', alpha=0.4, zorder=2)
ax.add_feature(cpf.LAND, color="lightgrey", zorder=2)
projection = ccrs.PlateCarree(central_longitude=0) 
ax.add_geometries(shdf.geometry,
                  projection,
                  facecolor="red",
                  edgecolor='k')

Reading the documentation of GeoPandas was very useful.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.