I am working with GOES-16 satellite imagery which come in netcdf files. I want to crop the file so that it only shows Puerto Rico and the outline. I can get it to crop, but latitude and longitude won't show up, neither will the coastline of Puerto Rico. Below is my sample code.
import netCDF4 as nc
from netCDF4 import Dataset
import cartopy as ccrs
import cartopy.mpl.ticker as cticker
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
#from mpl_toolkits.basemap import Basemap
import metpy
import numpy as np
import cv2
from cv2 import remap
import xarray
fn = 'OR_ABI-L2-ACMC-M6_G16_s20211730001172_e20211730003545_c20211730004314.nc'
ds = xarray.open_dataset(fn)
nc = nc.Dataset(fn)
variable = 'BCM'
bcm = ds['BCM'][:].data
bcm = np.clip(bcm, 0, 1)
geo_extent = nc.variables['geospatial_lat_lon_extent']
print(geo_extent)
#set latitude and longitude values from metadata
min_lon = float(geo_extent.geospatial_westbound_longitude)
max_lon = float(geo_extent.geospatial_eastbound_longitude)
min_lat = float(geo_extent.geospatial_southbound_latitude)
max_lat = float(geo_extent.geospatial_northbound_latitude)
dat = ds.metpy.parse_cf('BCM')
geos = dat.metpy.cartopy_crs
x = dat.x
y = dat.y
print(x.min())
fig = plt.figure(figsize=(12, 15))
ax = fig.add_subplot(1, 1, 1, projection=geos)
ax.set_extent([-68.07, -64.65, 16.7, 19.79], crs=geos)
ax.imshow(bcm, origin='upper', extent=(min_lon, max_lon, min_lat, max_lat), transform=geos)
ax.coastlines(resolution='10m', color='black', linewidth=1)
ax.add_feature(ccrs.cartopy.feature.COASTLINE)
plt.title('Puerto Rico', loc='left', fontweight='bold', fontsize=15)
plt.show()
Here is a picture of the map. [1]: https://i.sstatic.net/Y0qbY.png
I figured it out. I was using the incorrect syntax.
dat = ds.metpy.parse_cf('BCM')
geos = dat.metpy.cartopy_crs
x = dat.x
y = dat.y
pc = ccrs.crs.PlateCarree()
fig = plt.figure(figsize=(12, 15))
ax = fig.add_subplot(1, 1, 1, projection = pc)
ax.set_extent([-68.07, -64.65, 16.7, 19.79], crs = pc)
ax.imshow(bcm, origin='upper', extent=(x.min(), x.max(), y.min(), y.max()), transform=geos, interpolation = 'none')
ax.coastlines(resolution='10m', color='black', linewidth=1)
ax.add_feature(ccrs.cartopy.feature.COASTLINE)
ax.set_xticks((-68.07, -64.65), minor = 'True', crs= pc)
ax.set_yticks((16.7, 19.79), minor= 'True', crs= pc)
plt.title('Puerto Rico Clear-Sky Mask', loc='center', fontweight='bold', fontsize=15)
plt.show()