I have a ABI-L1b radiance file in the *.nc file format. I'm trying to cut out a selection of this file with a shapely polygon. I can plot both the netCDF4 file and polygon but I am unsure of how to cut out the data. Here is the code and a plot
#!/usr/bin/python
# See https://makersportal.com/blog/2018/11/25/goes-r-satellite-latitude-and-longitude-grid-projection-algorithm
from netCDF4 import Dataset
import matplotlib as mpl
mpl.use("TkAgg")
import pyproj
from functools import partial
from shapely.geometry import Point
from shapely.geometry import Polygon
from shapely.ops import transform
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import numpy as np
import os
# Site of the fire center
l = {"latitude": 40.109, "longitude": -122.789}
# Drawing km circle around what we're interested in
def geodesic_point_buffer(lat, lon, km):
proj_wgs84 = pyproj.Proj(init="epsg:4326")
# Azimuthal equidistant projection
aeqd_proj = "+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0"
project = partial(
pyproj.transform, pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)), proj_wgs84
)
buf = Point(0, 0).buffer(km * 1000) # distance in meters
return transform(project, buf).exterior.coords[:]
def lat_lon_reproj(nc_folder, nc_indx):
os.chdir(nc_folder)
full_direc = os.listdir()
nc_files = [ii for ii in full_direc if ii.endswith(".nc")]
g16_data_file = nc_files[nc_indx] # select .nc file
print(nc_files[nc_indx]) # print file name
# designate dataset
g16nc = Dataset(g16_data_file, "r")
var_names = [ii for ii in g16nc.variables]
var_name = var_names[0]
try:
band_id = g16nc.variables["band_id"][:]
band_id = " (Band: {},".format(band_id[0])
band_wavelength = g16nc.variables["band_wavelength"]
band_wavelength_units = band_wavelength.units
band_wavelength_units = "{})".format(band_wavelength_units)
band_wavelength = " {0:.2f} ".format(band_wavelength[:][0])
print("Band ID: {}".format(band_id))
print("Band Wavelength: {} {}".format(band_wavelength, band_wavelength_units))
except:
band_id = ""
band_wavelength = ""
band_wavelength_units = ""
# GOES-R projection info and retrieving relevant constants
proj_info = g16nc.variables["goes_imager_projection"]
lon_origin = proj_info.longitude_of_projection_origin
H = proj_info.perspective_point_height + proj_info.semi_major_axis
r_eq = proj_info.semi_major_axis
r_pol = proj_info.semi_minor_axis
# grid info
lat_rad_1d = g16nc.variables["x"][:]
lon_rad_1d = g16nc.variables["y"][:]
# data info
data = g16nc.variables[var_name][:]
data_units = g16nc.variables[var_name].units
data_time_grab = ((g16nc.time_coverage_end).replace("T", " ")).replace("Z", "")
data_long_name = g16nc.variables[var_name].long_name
# close file when finished
g16nc.close()
g16nc = None
# create meshgrid filled with radian angles
lat_rad, lon_rad = np.meshgrid(lat_rad_1d, lon_rad_1d)
# lat/lon calc routine from satellite radian angle vectors
lambda_0 = (lon_origin * np.pi) / 180.0
a_var = np.power(np.sin(lat_rad), 2.0) + (
np.power(np.cos(lat_rad), 2.0)
* (
np.power(np.cos(lon_rad), 2.0)
+ (((r_eq * r_eq) / (r_pol * r_pol)) * np.power(np.sin(lon_rad), 2.0))
)
)
b_var = -2.0 * H * np.cos(lat_rad) * np.cos(lon_rad)
c_var = (H ** 2.0) - (r_eq ** 2.0)
r_s = (-1.0 * b_var - np.sqrt((b_var ** 2) - (4.0 * a_var * c_var))) / (2.0 * a_var)
s_x = r_s * np.cos(lat_rad) * np.cos(lon_rad)
s_y = -r_s * np.sin(lat_rad)
s_z = r_s * np.cos(lat_rad) * np.sin(lon_rad)
lat = (180.0 / np.pi) * (
np.arctan(
((r_eq * r_eq) / (r_pol * r_pol))
* ((s_z / np.sqrt(((H - s_x) * (H - s_x)) + (s_y * s_y))))
)
)
lon = (lambda_0 - np.arctan(s_y / (H - s_x))) * (180.0 / np.pi)
return (
lon,
lat,
data,
data_units,
data_time_grab,
data_long_name,
band_id,
band_wavelength,
band_wavelength_units,
var_name,
)
nc_folder = "/path/to/ABI-L1b-radiances" # define folder where .nc files are located
nc_folder = "/home/n/Documents/Research/Fire-Prediction/data/2019/248/00"
file_indx = (
0 # be sure to pick the correct file. Make sure the file is not too big either,
)
# main data grab from function above
(
lon,
lat,
data,
data_units,
data_time_grab,
data_long_name,
band_id,
band_wavelength,
band_units,
var_name,
) = lat_lon_reproj(nc_folder, file_indx)
flat_lon = lon.flatten()
flat_lat = lat.flatten()
flat_data = data.flatten()
b = geodesic_point_buffer(l["latitude"], l["longitude"], 5)
b_x, b_y = zip(*b)
poly_km = Polygon(b)
bbox = [np.min(lon), np.min(lat), np.max(lon), np.max(lat)] # set bounds for plotting
# figure routine for visualization
fig = plt.figure(figsize=(9, 4), dpi=200)
n_add = 0
m = Basemap(
resolution="i",
projection="cyl",
area_thresh=50000,
llcrnrlon=l["longitude"] - 1,
llcrnrlat=l["latitude"] - 1,
urcrnrlon=l["longitude"] + 1,
urcrnrlat=l["latitude"] + 1,
)
m.drawcoastlines(linewidth=0.5)
m.drawcountries(linewidth=0.25)
m.pcolormesh(lon.data, lat.data, data, latlon=True)
parallels = np.linspace(np.min(lat), np.max(lat), 5)
m.drawparallels(parallels, labels=[True, False, False, False])
meridians = np.linspace(np.min(lon), np.max(lon), 5)
m.drawmeridians(meridians, labels=[False, False, False, True])
cb = m.colorbar()
data_units = ((data_units.replace("-", "^{-")).replace("1", "1}")).replace("2", "2}")
plt.rc("text", usetex=True)
cb.set_label(r"%s $ \left[ \mathrm{%s} \right] $" % (var_name, data_units))
plt.title(
"{0}{2}{3}{4} on {1}".format(
data_long_name, data_time_grab, band_id, band_wavelength, band_units
)
)
plt.tight_layout()
m_bx, m_by = m(b_x, b_y)
m.plot(b_x, b_y, c="r", linewidth=0.1)
plt.show()
I get the following image:
So I need to cut out the data in the red circle from the image. I've tried looking into GDAL python cut geotiff image with geojson file. So I convert the .nc file to .tif using gtranslate then using rasterio intersect the geometry with the .tif file. However, I keep getting that the geometry does not intersect the raster file. Am I going about this the wrong way?
Here is a dropbox link to the .nc file: https://www.dropbox.com/s/nyj6zpdvubgqppb/glm.nc?dl=0
