Copied from the private repo

georadii/rsp.py

@@ -193,7 +193,7 @@ def date_check(self, date):
 				message = 'Error [RSP_arcsix]: date %s is not avaiable.'
 				raise OSError(message)
 			else:
-				print('Camera images for [%s]' % self._flights[date]['description'])
+				print('RSP images for [%s]' % self._flights[date]['description'])
 
 	def load_rsp_h5(self, start_time, end_time, location='.'):
 		self.h5_allfiles = self.locate_allh5(location=location)
@@ -219,6 +219,7 @@ def load_rsp_h5(self, start_time, end_time, location='.'):
 				wvls    = h5_file['Data']['Wavelength'][...][0:nband]
 
 		glat_arr, glon_arr = np.array([], dtype=float), np.array([], dtype=float)
+		vza_arr, vaa_arr = np.array([], dtype=float), np.array([], dtype=float)
 		plat_arr, plon_arr = np.array([], dtype=float), np.array([], dtype=float)
 		time_arr = np.array([], dtype='datetime64')
 		intens_arr = np.zeros((0, 9), dtype=float)
@@ -229,24 +230,34 @@ def load_rsp_h5(self, start_time, end_time, location='.'):
 				nband   = h5_file['dim_Bands'][...].shape[0]
 				ground_latitude  = h5_file['Geometry']['Ground_Latitude'][...][0:2, 0:nscan, 0:nsector]
 				ground_longitude = h5_file['Geometry']['Ground_Longitude'][...][0:2, 0:nscan, 0:nsector]
+				viewing_azimuth  = h5_file['Geometry']['Viewing_Azimuth'][...][0:2, 0:nscan, 0:nsector]
+				viewing_zenith   = 180. - h5_file['Geometry']['Viewing_Zenith'][...][0:2, 0:nscan, 0:nsector]
 				platform_latitude  = h5_file['Platform']['Platform_Latitude'][...][0:nscan]
 				platform_longitude = h5_file['Platform']['Platform_Longitude'][...][0:nscan]
 				# time_scan        = h5_file['Platform']['Fraction_of_Day'][...][0:nscan]
 				time_pixel       = h5_file['Geometry']['Measurement_Time'][...][0:2, 0:nscan, 0:nsector]
+				solar_cons       = h5_file['Calibration']['Solar_Constant'][...][0:nband]
 				intensity_1   = h5_file['Data']['Intensity_1'][...][0:nscan, 0:nsector, 0:nband]
-				intensity_2   = h5_file['Data']['Intensity_2'][...][0:nscan, 0:nsector, 0:nband]            
-				intensity_avg = (intensity_1 + intensity_2) / 2
+				intensity_2   = h5_file['Data']['Intensity_2'][...][0:nscan, 0:nsector, 0:nband]
+				intensity_avg = (intensity_1 + intensity_2) / 2.
+				radiance      = intensity_avg*solar_cons[np.newaxis, np.newaxis, 0:nband]/(np.pi*1000.)
 				glat_arr = np.append(glat_arr, ground_latitude[0, 0:nscan, 0:nsector].flatten())
 				glon_arr = np.append(glon_arr, ground_longitude[0, 0:nscan, 0:nsector].flatten())
-				plat_arr = np.append(plat_arr, platform_latitude[0:nscan].flatten())
-				plon_arr = np.append(plon_arr, platform_longitude[0:nscan].flatten())
+				vza_arr  = np.append(vza_arr, viewing_zenith[0, 0:nscan, 0:nsector].flatten())
+				vaa_arr  = np.append(vaa_arr, viewing_azimuth[0, 0:nscan, 0:nsector].flatten())
+				plat_arr = np.append(plat_arr, np.repeat(platform_latitude[0:nscan], nsector).flatten())
+				plon_arr = np.append(plon_arr, np.repeat(platform_longitude[0:nscan], nsector).flatten())
 				# time_arr = np.append(time_arr, np.array([datetime.strptime(start_time_str.split(' ')[0] + ' 00:00', '%Y-%m-%d %H:%M') + datetime.timedelta(days=fofd) for fofd in time_scan[0:nscan]]))
 				time_arr = np.append(time_arr, Time(time_pixel[0, 0:nscan, 0:nsector].flatten(), format='mjd').to_datetime())
-				intens_arr = np.append(intens_arr, intensity_avg[0:nscan, 0:nsector, 0:nband].reshape(nscan*nsector, nband), axis=0)
+				intens_arr = np.append(intens_arr, radiance[0:nscan, 0:nsector, 0:nband].reshape(nscan*nsector, nband), axis=0)
 		ipixels = np.where((st_dt + self.t_offset < time_arr) & (time_arr < en_dt + self.t_offset))[0]
 		colpix  = np.zeros((len(ipixels), 1, 3))
 		latpix  = np.zeros((len(ipixels), 1))
 		lonpix  = np.zeros((len(ipixels), 1))
+		platpix  = np.zeros((len(ipixels), 1))
+		plonpix  = np.zeros((len(ipixels), 1))
+		vzapix  = np.zeros((len(ipixels), 1))
+		vaapix  = np.zeros((len(ipixels), 1))
 		timepix = np.zeros((len(ipixels), 1), dtype='datetime64')
 		# print(ipixels)
 		# print(wvls)
@@ -256,15 +267,21 @@ def load_rsp_h5(self, start_time, end_time, location='.'):
 		colpix[:, 0, 2]  = intens_arr[ipixels, 1]
 		latpix[:, 0]     = glat_arr[ipixels]
 		lonpix[:, 0]     = glon_arr[ipixels]
+		vzapix[:, 0]     = vza_arr[ipixels]
+		vaapix[:, 0]     = vaa_arr[ipixels]
+		platpix[:, 0]    = plat_arr[ipixels]
+		plonpix[:, 0]    = plon_arr[ipixels]
 		# timepix[:, 0]    = time_arr[ipixels]
-		img = {'data': colpix, 'type': 'radiance', 'unit': 'radiance'}
-		latlon_meta = {'longeo': lonpix, 'latgeo': latpix}#, 'timepix': timepix}
+		img = {'data': colpix, 'type': 'radiance', 'unit': 'radiance', 'wavelength': wvls[np.array([3, 2, 1])]}
+		# latlon_meta = {'longeo': lonpix, 'latgeo': latpix}#, 'timepix': timepix}
+		latlon_meta = {'longeo': lonpix, 'latgeo': latpix, 'vza': vzapix, 'vaa': vaapix, 'lonplat': plonpix, 'latplat': platpix}#, 'timepix': timepix}
 		return img, latlon_meta
 
 	def locate_allh5(self, location='.'):
 		yyyy, mm, dd = self.date.split('-')
 		if True:
-			path = os.path.join(location, 'ARCSIX-RSP-L1C_P3B_%s%s%s_R01/*.h5' % (yyyy, mm, dd))
+			# path = os.path.join(location, 'ARCSIX-RSP-L1C_P3B_%s%s%s_R01/*.h5' % (yyyy, mm, dd))
+			path = os.path.join(location, 'RSP1-L1C_P3_%s%s%s_R*/*.h5' % (yyyy, mm, dd)) #RSP1-L1C_P3_20240605_R02
 			# path = location
 		h5_allfiles = glob.glob(path, recursive=True)
 		if len(h5_allfiles) == 0:
@@ -275,10 +292,12 @@ def locate_allh5(self, location='.'):
 
 	def find_files_in_range(self, h5_files, start_time, end_time):
 		file_timestamps = []
-		print(h5_files, start_time, end_time)
+		# print(h5_files, start_time, end_time)
 		for file in h5_files:
-			timestamp_str = os.path.basename(file).split('_')[2]
-			timestamp = datetime.datetime.strptime(timestamp_str, "%Y%m%d%H%M%S")
+			# timestamp_str = os.path.basename(file).split('_')[2]
+			# timestamp = datetime.datetime.strptime(timestamp_str, "%Y%m%d%H%M%S")
+			timestamp_str = os.path.basename(file).split('_')[2] #RSP1-P3_L1C-RSPCOL-CollocatedRadiances_20240605T113433Z_V002-20241114T193424Z.h5
+			timestamp = datetime.datetime.strptime(timestamp_str, "%Y%m%dT%H%M%SZ")
 			file_timestamps.append((file, timestamp))
 		files_all = sorted([file for file, timestamp in file_timestamps])
 		files_before_start = sorted([file for file, timestamp in file_timestamps if start_time > timestamp])

georadii/util.py

@@ -4,6 +4,11 @@
 import h5py
 from astropy.io import fits
 import cv2
+import datetime
+import netCDF4 as nc
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+from matplotlib.ticker import FixedLocator
 
 # Read housekeeping file (HDF5 binary data format)
 def read_hsk_camp2ex(hsk_filename):
@@ -20,6 +25,10 @@ def read_hsk_camp2ex(hsk_filename):
                 'pit': h5f['pitch_angle'][...]  }
     return hsk_data
 
+def jd_to_doy(jd, year):
+    jan1_jd = 1721425.5 + 365 * (year - 1) + int((year - 1) / 4)
+    return jd - jan1_jd + 1
+
 def read_hsk_arcsix(hsk_filename):
     print('Reading {}'.format(hsk_filename))
     h5f = h5py.File(hsk_filename, 'r')
@@ -51,74 +60,37 @@ def reshapeFITS(img):
     return newimg
 
 # Read image file (Fits file format)
-def read_fits(fits_filename, flipud=False, fliplr=True, mask_fits_filename=None):
+def read_fits(fits_filename, flipud=False, fliplr=True, mask_fits_filename=None, header_only=False):
 	if not os.path.exists(fits_filename):
 		print('Error: {} not found.'.format(fits_filename))
 		sys.exit()
 	print('Reading {}'.format(fits_filename))
 	handle = fits.open(fits_filename)
-	fimg = reshapeFITS(handle[0].data)
-	if flipud:
-		fimg = np.flipud(fimg)
-	if fliplr:
-		fimg = np.fliplr(fimg)
-	fheader = handle[0].header
-
-	if mask_fits_filename is not None:
-		print('Reading {}'.format(mask_fits_filename))
-		handle_msk = fits.open(mask_fits_filename)
-		fmsk = reshapeFITS(handle_msk[0].data)
+	if not header_only:
+		fimg = reshapeFITS(handle[0].data)
 		if flipud:
-			fmsk = np.flipud(fmsk)
+			fimg = np.flipud(fimg)
 		if fliplr:
-			fmsk = np.fliplr(fmsk)
-		fheader_msk = handle_msk[0].header
+			fimg = np.fliplr(fimg)
+
+		if mask_fits_filename is not None:
+			print('Reading {}'.format(mask_fits_filename))
+			handle_msk = fits.open(mask_fits_filename)
+			fmsk = reshapeFITS(handle_msk[0].data)
+			if flipud:
+				fmsk = np.flipud(fmsk)
+			if fliplr:
+				fmsk = np.fliplr(fmsk)
+			fheader_msk = handle_msk[0].header
 
-		fimg[fmsk > 0.] = np.nan
+			fimg[fmsk > 0.] = np.nan
+	else:
+		fimg = None
+
+	fheader = handle[0].header
 
 	return fimg, fheader
 
-# # Calculate Direct-Cosine-Matrix to convert NED(North-East-Down) coordinate to
-# #  camera coordinate
-# # Usage: d(Cam) = R*d(NED)
-# def R_NED2Cam(roll, pitch, yaw): #roll, pitch, yaw [radian]
-# 	cr = np.cos(roll)
-# 	sr = np.sin(roll)
-# 	cp = np.cos(pitch)
-# 	sp = np.sin(pitch)
-# 	cy = np.cos(yaw)
-# 	sy = np.sin(yaw)
-# 	Rr = np.matrix([[ 1.,  0.,  0.],
-# 					[ 0.,  cr,  sr],
-# 					[ 0., -sr,  cr]])
-# 	Rp = np.matrix([[ cp,  0., -sp],
-# 					[ 0.,  1.,  0.],
-# 					[ sp,  0.,  cp]])
-# 	Ry = np.matrix([[ cy,  sy,  0.],
-# 					[-sy,  cy,  0.],
-# 					[ 0.,  0.,  1.]])
-# 	R = np.matmul(Rr, np.matmul(Rp, Ry))
-# 	return R
-
-# def new_get_image(self, tile):
-# 	import six
-# 	from PIL import Image
-# 	if six.PY3:
-# 		from urllib.request import urlopen, Request
-# 	else:
-# 		from urllib2 import urlopen
-# 	url = self._image_url(tile)  # added by H.C. Winsemius
-# 	req = Request(url) # added by H.C. Winsemius
-# 	req.add_header('User-agent', 'your bot 0.1')
-# 	# fh = urlopen(url)  # removed by H.C. Winsemius
-# 	fh = urlopen(req)
-# 	im_data = six.BytesIO(fh.read())
-# 	fh.close()
-# 	img = Image.open(im_data)
-# 	img = img.convert(self.desired_tile_form)
-
-# 	return img, self.tileextent(tile), 'lower'
-
 def find_center(xs, ys):
 	A = np.vstack((xs, ys, np.ones((len(xs))))).T
 	v = -(xs ** 2 + ys ** 2)
@@ -220,3 +192,421 @@ def matched_points_to_latlon(m1, m2, lon_xx, lat_yy):
     latm2 = (1. - rym2)*lat_yy[iym2, ixm2] + rym2*lat_yy[iym2 + 1, ixm2]
     return lonm1, latm1, lonm2, latm2
 
+def calc_viewing_angles(lon, lat, lon_air, lat_air, alt_air):
+    from pyproj import Geod as geod
+    g = geod(ellps='WGS84')
+    lon_air_tmp = lon_air*np.ones_like(lon)
+    lat_air_tmp = lat_air*np.ones_like(lat)
+    fwd_az1, fwd_az2, dist = g.inv(lon_air_tmp, lat_air_tmp, lon, lat)
+    theta = np.arctan(dist/alt_air)*180./np.pi
+    phi   = fwd_az1
+    phi[phi < 0.] += 360.
+    return theta, phi # vza, vaa (deg)
+
+def get_spec_resp(spec_resp_txt, instrument='cam'):
+    if instrument == 'cam':
+        nch = 3
+    elif instrument == 'rsp':
+        nch = 9
+    wavelengths = []
+    response = [[] for _ in range(nch)]
+    with open(spec_resp_txt, "r") as file:
+        lines = file.readlines()
+    for line in lines[1:]: # Skip the header
+        values = line.strip().split()
+        if len(values) == nch + 1:
+            wavelengths.append(float(values[0]))
+            for i in range(nch):
+                response[i].append(float(values[i + 1]))
+    wavelengths = np.array(wavelengths)
+    spec_resp = np.zeros((nch, len(wavelengths)))
+    for i in range(nch):
+        spec_resp[i, :] = np.array(response[i])
+        spec_resp[i, :] /= np.trapz(spec_resp[i, :], x=wavelengths)
+        # print('Center wavelength: %9.2f nm' % np.trapz(wavelengths*spec_resp[i, :], x=wavelengths))
+    return wavelengths, spec_resp, nch
+
+def get_ssfr_fluxdn(f_RA, fhsk, start_time, end_time, flux_source='ssfr', instrument='cam', spec_resp_txt="./response_ARCSIX.txt"):
+    ### Open SSFR file ###
+    with h5py.File(f_RA, 'r') as f:
+        tmhr_all = f['tmhr'][...]
+        if flux_source == 'ssfr':
+            rad_all = f['zen/flux'][...]
+            wvl = f['zen/wvl'][...] # nm
+            tmhr_ssfr = f['tmhr'][...]
+        else:
+            msg = 'Unknown flux source: %s' % flux_source
+            raise ValueError(msg)
+
+    hsk_data = read_hsk_arcsix(fhsk)
+    tmhr_hsk = hsk_data['hrs']
+    pit = np.interp(tmhr_ssfr, tmhr_hsk, hsk_data['pit'])
+    rol = np.interp(tmhr_ssfr, tmhr_hsk, hsk_data['rol'])
+
+    level = np.sqrt(pit**2. + rol**2.) < 2.5 # check if the aircraft is not tilted too much
+
+    st_dt = datetime.datetime.strptime(start_time, '%H:%M:%S')
+    en_dt = datetime.datetime.strptime(end_time,   '%H:%M:%S')
+
+    st_hr = st_dt.hour + st_dt.minute/60. + st_dt.second/3600.
+    en_hr = en_dt.hour + en_dt.minute/60. + en_dt.second/3600.
+
+    target = level & (st_hr <= tmhr_all) & (tmhr_all < en_hr) # time and attitude filter
+    tmhr = tmhr_all[target]
+    rad  = rad_all[target]
+
+    ### Open camera spectral response ###
+    wavelengths, spec_resp, nch = get_spec_resp(spec_resp_txt, instrument=instrument)
+
+    # Match the SSFR data to the camera spectral response function wavelength increments
+    flux_ssfr = np.zeros((len(tmhr), len(wavelengths)))
+    for it in range(len(tmhr)):
+        flux_ssfr[it, :] = np.interp(wavelengths, wvl, rad[it, :])
+
+    # Calculate the flux down corresponding to each camera channel
+    f_resp = np.zeros((nch, len(tmhr)))
+    for ich in range(nch):
+        f_resp[ich, :] = np.trapezoid(spec_resp[ich, :][np.newaxis, :]*flux_ssfr, x=wavelengths, axis=1)
+
+    # Average temporally
+    flux_down = np.nanmean(f_resp[:, :], axis=1)
+    for ich in range(nch):
+        print('Flux down (ch=%d): %9.4f (W/m^2/nm)' % (ich, flux_down[ich]))
+
+    return flux_down
+
+def write_surface_grid_to_nc(output_ncfile, datout, lonxx, latyy, vzagrid, vaagrid, ncount, date, tact, reflectance=None, flxdn=None, wvlc=None, sza=None, saa=None):
+    print('Writing to %s' % (output_ncfile))
+    with nc.Dataset(output_ncfile, "w", format="NETCDF4") as ncfile:
+        # Define dimensions
+        lat_dim = ncfile.createDimension("lat", datout.shape[0])
+        lon_dim = ncfile.createDimension("lon", datout.shape[1])
+        ch_dim  = ncfile.createDimension("ch",  datout.shape[2] if len(datout.shape) == 3 else 1)
+
+        # Define variables
+        lon_var = ncfile.createVariable("longitude", "f4", ("lat", "lon"))
+        lat_var = ncfile.createVariable("latitude", "f4", ("lat", "lon"))
+        vza_var = ncfile.createVariable("vza", "f4", ("lat", "lon"))
+        vaa_var = ncfile.createVariable("vaa", "f4", ("lat", "lon"))
+        count_var = ncfile.createVariable("count", "i4", ("lat", "lon"))
+        radiance_var = ncfile.createVariable("radiance", "f4", ("lat", "lon", 'ch'), fill_value=np.nan)
+        if reflectance is not None:
+            reflectance_var = ncfile.createVariable("hdrf", "f4", ("lat", "lon", 'ch'), fill_value=np.nan)
+        if flxdn is not None:
+            down_flux_var = ncfile.createVariable("fluxdn", "f4", ('ch'))
+        if wvlc is not None:
+            center_wvl_var = ncfile.createVariable("wvlc", "f4", ('ch'))
+        if sza is not None:
+            sza_var = ncfile.createVariable("sza", "f4")
+        if saa is not None:
+            saa_var = ncfile.createVariable("saa", "f4")
+        yyyy_var = ncfile.createVariable("year", "i4")
+        mm_var = ncfile.createVariable("month", "i4")
+        dd_var = ncfile.createVariable("day", "i4")
+        time_var = ncfile.createVariable("time", "f4")
+
+        # Add attributes
+        lon_var.long_name = "Longitude"
+        lon_var.units = "degree"
+        lon_var.description = "Longitude at the surface level."
+
+        lat_var.long_name = "Latitude"
+        lat_var.units = "degree"
+        lat_var.description = "Latitude at the surface level."
+
+        vza_var.long_name = "Viewing zenith angle"
+        vza_var.units = "degree"
+        vza_var.description = "Viewing zenith angle at the surface level."
+
+        vaa_var.long_name = "Viewing azimuth angle"
+        vaa_var.units = "degree"
+        vaa_var.description = "Viewing azimuth angle at the surface level."
+
+        count_var.long_name = "Count"
+        count_var.units = "dimensionless"
+        count_var.description = "Number of image pixels averaged over to generate radiance for a given grid."
+
+        radiance_var.long_name = "Radiance"
+        radiance_var.units = "W/m^2/nm/sr"
+        radiance_var.description = "Radiance as a function of viewing zenith and azimuth angles."
+
+        if reflectance is not None:
+            reflectance_var.long_name = "HDRF"
+            reflectance_var.units = "dimensionless"
+            reflectance_var.description = "HDRF (reflectance) as a function of viewing zenith and azimuth angles."
+
+        if flxdn is not None:
+            down_flux_var.long_name = "Downward irradiance"
+            down_flux_var.units = "W/m^2/nm"
+            down_flux_var.description = "Downward irradiance averaged over the flight leg. The camera response function is applied."
+
+        if wvlc is not None:
+            center_wvl_var.long_name = "Center wavelength"
+            center_wvl_var.units = "nm"
+            center_wvl_var.description = "Center wavelength of the camera's spectral response function."
+
+        if sza is not None:
+            sza_var.long_name = "Solar Zenith Angle"
+            sza_var.units = "degrees"
+            sza_var.description = "Solar zenith angle averaged over the flight leg."
+
+        if saa is not None:
+            saa_var.long_name = "Solar Azimuth Angle"
+            saa_var.units = "degrees"
+            saa_var.description = "Solar azimuth angle averaged over the flight leg."
+
+        yyyy_var.long_name = "Year"
+        yyyy_var.units = "dimensionless"
+        yyyy_var.description = "Year of the flight"
+
+        mm_var.long_name = "Month"
+        mm_var.units = "dimensionless"
+        mm_var.description = "Month of the flight"
+
+        dd_var.long_name = "Day"
+        dd_var.units = "dimensionless"
+        dd_var.description = "Day of the flight"
+
+        time_var.long_name = "Time"
+        time_var.units = "decimal hour"
+        time_var.description = "Time at which the camera data was obtained. It may contain some offset due to time sync error."
+
+        ncfile.title = "Gridded image"
+        ncfile.description = "Hemispherical radiance derived from a airborne downward hemispherical camera on %s at %s UTC." % (date, tact.strftime("%H%M%S"))
+
+        # Assign data to variables
+        lon_var[:, :] = lonxx
+        lat_var[:, :] = latyy
+        vza_var[:, :] = vzagrid
+        vaa_var[:, :] = vaagrid
+        count_var[:, :] = ncount
+        radiance_var[:, :, :] = datout if len(datout.shape) == 3 else datout[:, :, np.newaxis]
+        if reflectance is not None:
+            reflectance_var[:, :, :] = reflectance if len(reflectance.shape) == 3 else reflectance[:, :, np.newaxis]
+        if flxdn is not None:
+            down_flux_var[:] = flxdn
+        if wvlc is not None:
+            center_wvl_var[:] = wvlc
+        if sza is not None:
+            sza_var.assignValue(sza)
+        if saa is not None:
+            saa_var.assignValue(saa)
+        yyyy_var.assignValue(date.split('-')[0])
+        mm_var.assignValue(date.split('-')[1])
+        dd_var.assignValue(date.split('-')[2])
+        time_var.assignValue(tact.hour + tact.minute/60. + tact.second/3600. + tact.microsecond/(3600*1e6))
+
+def write_angular_grid_to_nc(output_ncfile, datout, zen_yy, razi_xx, date, start_time, end_time, reflectance=None, azi_xx=None, sza=None, saa=None, flxdn=None, wvlc=None, nimg=None, alt=None):
+    print('Writing to %s' % (output_ncfile))
+    with nc.Dataset(output_ncfile, "w", format="NETCDF4") as ncfile:
+        # Define dimensions
+        zenith_dim = ncfile.createDimension("zenith", datout.shape[0])
+        azimuth_dim = ncfile.createDimension("azimuth", datout.shape[1])
+        ch_dim = ncfile.createDimension("ch", datout.shape[2] if len(datout.shape) == 3 else 1)
+
+        # Define variables
+        zenith_var = ncfile.createVariable("vza", "f4", ("zenith", "azimuth"))
+        relative_azimuth_var = ncfile.createVariable("raa", "f4", ("zenith", "azimuth"))
+        if azi_xx is not None:
+            azimuth_var = ncfile.createVariable("vaa", "f4", ("zenith", "azimuth"))
+        radiance_var = ncfile.createVariable("radiance", "f4", ("zenith", "azimuth", 'ch'), fill_value=np.nan)
+        if reflectance is not None:
+            reflectance_var = ncfile.createVariable("hdrf", "f4", ("zenith", "azimuth", 'ch'), fill_value=np.nan)
+        if flxdn is not None:
+            down_flux_var = ncfile.createVariable("fluxdn", "f4", ('ch'))
+        if wvlc is not None:
+            center_wvl_var = ncfile.createVariable("wvlc", "f4", ('ch'))
+        if sza is not None:
+            sza_var = ncfile.createVariable("solar_zenith_angle", "f4")
+        if saa is not None:
+            saa_var = ncfile.createVariable("solar_azimuth_angle", "f4")
+        if alt is not None:
+            alt_var = ncfile.createVariable("alt", "f4")
+        yyyy_var = ncfile.createVariable("year", "i4")
+        mm_var = ncfile.createVariable("month", "i4")
+        dd_var = ncfile.createVariable("day", "i4")
+        st_time_var = ncfile.createVariable("start time", "f4")
+        en_time_var = ncfile.createVariable("end time", "f4")
+        if nimg is not None:
+            nimg_var = ncfile.createVariable("image number", "i4")
+
+        # Add attributes
+        zenith_var.long_name = "Viewing Zenith Angle"
+        zenith_var.units = "degrees"
+        zenith_var.description = "Angle from the vertical axis to the line of sight. Nadir = 0 degrees."
+
+        relative_azimuth_var.long_name = "Relative Azimuth Angle"
+        relative_azimuth_var.units = "degrees"
+        relative_azimuth_var.description = "Horizontal angle between the line of sight and the principal plane. Sunglint should be at 0 degrees."
+
+        if azi_xx is not None:
+            azimuth_var.long_name = "Viewing Azimuth Angle"
+            azimuth_var.units = "degrees"
+            azimuth_var.description = "Angle from the north direction to the line of sight. North = 0 degrees."
+
+        radiance_var.long_name = "Radiance"
+        radiance_var.units = "W/m^2/nm/sr"
+        radiance_var.description = "Radiance as a function of viewing zenith and azimuth angles."
+
+        if reflectance is not None:
+            reflectance_var.long_name = "HDRF"
+            reflectance_var.units = "dimensionless"
+            reflectance_var.description = "HDRF (reflectance) as a function of viewing zenith and azimuth angles."
+
+        if flxdn is not None:
+            down_flux_var.long_name = "Downward irradiance"
+            down_flux_var.units = "W/m^2/nm"
+            down_flux_var.description = "Downward irradiance averaged over the flight leg. The camera response function is applied."
+
+        if wvlc is not None:
+            center_wvl_var.long_name = "Center wavelength"
+            center_wvl_var.units = "nm"
+            center_wvl_var.description = "Center wavelength of the camera's spectral response function."
+
+        if sza is not None:
+            sza_var.long_name = "Solar Zenith Angle"
+            sza_var.units = "degrees"
+            sza_var.description = "Solar zenith angle averaged over the flight leg."
+
+        if saa is not None:
+            saa_var.long_name = "Solar Azimuth Angle"
+            saa_var.units = "degrees"
+            saa_var.description = "Solar azimuth angle averaged over the flight leg."
+
+        if alt is not None:
+            alt_var.long_name = "Altitude"
+            alt_var.units = "m"
+            alt_var.description = "Altitude of the aircraft averaged over the flight leg."
+
+        yyyy_var.long_name = "Year"
+        yyyy_var.units = "dimensionless"
+        yyyy_var.description = "Year of the flight"
+
+        mm_var.long_name = "Month"
+        mm_var.units = "dimensionless"
+        mm_var.description = "Month of the flight"
+
+        dd_var.long_name = "Day"
+        dd_var.units = "dimensionless"
+        dd_var.description = "Day of the flight"
+
+        st_time_var.long_name = "Start time"
+        st_time_var.units = "decimal hour"
+        st_time_var.description = "Beginning of the time frame over which the camera data was averaged."
+
+        en_time_var.long_name = "End time"
+        en_time_var.units = "decimal hour"
+        en_time_var.description = "End of the time frame over which the camera data was averaged."
+
+        if nimg is not None:
+            nimg_var.long_name = "Image number"
+            nimg_var.units = "dimensionless"
+            nimg_var.description = "Number of images averaged to produce radiance/HDRF."
+
+        ncfile.title = "HDRF Data"
+        ncfile.description = "HDRF (reflectance) values derived from %d airborne downward hemispherical camera images on %s between %s - %s UTC." % (nimg, date, start_time, end_time)
+
+        # Assign data to variables
+        zenith_var[:, :] = zen_yy
+        relative_azimuth_var[:, :] = razi_xx % 360.
+        if azi_xx is not None:
+            azimuth_var[:, :] = azi_xx % 360.
+        radiance_var[:, :, :] = datout if len(datout.shape) == 3 else datout[:, :, np.newaxis]
+        if reflectance is not None:
+            reflectance_var[:, :, :] = reflectance if len(reflectance.shape) == 3 else reflectance[:, :, np.newaxis]
+        if flxdn is not None:
+            down_flux_var[:] = flxdn
+        if wvlc is not None:
+            center_wvl_var[:] = wvlc
+        if sza is not None:
+            sza_var.assignValue(sza)
+        if saa is not None:
+            saa_var.assignValue(saa)
+        if alt is not None:
+            alt_var.assignValue(alt)
+        yyyy_var.assignValue(date.split('-')[0])
+        mm_var.assignValue(date.split('-')[1])
+        dd_var.assignValue(date.split('-')[2])
+        st_time_var.assignValue(start_time)
+        en_time_var.assignValue(end_time)
+        if nimg is not None:
+            nimg_var.assignValue(nimg)
+
+def plot_surface_and_angular_grid_image(fn_out, lon_xx, lat_yy, ref_surface, rel_azimuth, zenith, ref_angular, tact):
+    fig  = plt.figure(figsize=(8, 4.5))
+
+    cartopy_proj = ccrs.Orthographic(central_longitude=np.nanmean(lon_xx), central_latitude=np.nanmean(lat_yy),)
+    ax = fig.add_subplot(121, projection=cartopy_proj)
+    xmin, xmax = np.nanmin(lon_xx[~np.isnan(ref_surface)]), np.nanmax(lon_xx[~np.isnan(ref_surface)])
+    ymin, ymax = np.nanmin(lat_yy[~np.isnan(ref_surface)]), np.nanmax(lat_yy[~np.isnan(ref_surface)])
+    cp = ax.pcolormesh(lon_xx, lat_yy, ref_surface, transform=ccrs.PlateCarree(), zorder=10)
+    ax.set_extent([xmin, xmax, ymin, ymax])
+    g1 = ax.gridlines(lw=0.5, color='gray', draw_labels=True, ls='-')
+    g1.xlocator = FixedLocator(np.arange(-180, 180.1, 0.5*10.**(np.round(np.log10(np.abs(xmax - xmin))))))
+    g1.ylocator = FixedLocator(np.arange(-90.0, 89.9, 0.5*10.**(np.round(np.log10(np.abs(ymax - ymin))))))
+    g1.top_labels = False
+    g1.right_labels = False
+    ax.annotate('Surface grid reflectance', (0.01, 1.05), xycoords='axes fraction')
+    cbar = fig.colorbar(cp, ax=ax, orientation='horizontal')
+    cbar.set_label('Reflectance')
+
+    ax2 = fig.add_subplot(122, projection='polar')
+    ax2.set_theta_direction(-1)
+    ax2.set_theta_offset(np.pi/2.)
+    cw = ax2.pcolormesh(rel_azimuth*np.pi/180., zenith, ref_angular)
+    ax2.set_rticks([30., 60., 90.])
+    ax2.annotate('Angular reflectance', (-0.10, 1.05), xycoords='axes fraction')
+    cbar2 = fig.colorbar(cw, ax=ax2, orientation='horizontal')
+    cbar2.set_label('Reflectance')
+
+    fig.suptitle('Camera image: ' + datetime.datetime.strftime(tact, '%Y-%m-%d %H:%M:%S.%f'), fontsize=9)
+    fig.tight_layout()
+
+    fig.savefig(fn_out, dpi=300)
+
+def plot_angular_grid_rad_and_ref(fn_out, rel_azimuth, zenith, rad, ref, sza, flxdn, start_time, end_time=None, nimg=None, alt=None, meta={}):
+    meta1 = meta.get('rad', {})
+    meta2 = meta.get('ref', {})
+    vmin1 = meta1.get('vmin', None)
+    vmax1 = meta1.get('vmax', None)
+    cmap1 = meta1.get('cmap', 'viridis')
+    vmin2 = meta2.get('vmin', None)
+    vmax2 = meta2.get('vmax', None)
+    cmap2 = meta2.get('cmap', 'viridis')
+
+    fig  = plt.figure(figsize=(8, 4.5))
+
+    ax1 = fig.add_subplot(121, projection='polar')
+    ax1.set_theta_direction(-1)
+    ax1.set_theta_offset(np.pi/2.)
+    cw = ax1.pcolormesh(rel_azimuth*np.pi/180., zenith, rad, vmin=vmin1, vmax=vmax1, cmap=cmap1)
+    ax1.set_rticks([30., 60., 90.])
+    ax1.annotate('Radiance', (-0.10, 1.05), xycoords='axes fraction')
+    cbar1 = fig.colorbar(cw, ax=ax1, orientation='horizontal')
+    cbar1.set_label(r'Radiance $ \rm (W/m^2/nm/sr) $')
+
+    ax1.annotate('Average SZA = %5.2f deg' % (sza), (0.01, -0.19), xycoords='axes fraction', fontsize=9)
+
+    ax2 = fig.add_subplot(122, projection='polar')
+    ax2.set_theta_direction(-1)
+    ax2.set_theta_offset(np.pi/2.)
+    cw = ax2.pcolormesh(rel_azimuth*np.pi/180., zenith, ref, vmin=vmin2, vmax=vmax2, cmap=cmap2)
+    ax2.set_rticks([30., 60., 90.])
+    ax2.annotate('Reflectance', (-0.10, 1.05), xycoords='axes fraction')
+    if alt is not None: ax2.annotate('Alt: %7.1f m' % alt, ( 0.70, -0.07), xycoords='axes fraction', fontsize=9)
+    cbar2 = fig.colorbar(cw, ax=ax2, orientation='horizontal')
+    cbar2.set_label('Reflectance')
+
+    ax2.annotate(r'Average Downward flux = %6.3f $ \rm (W/m^2/nm) $' % (flxdn), (0.01, -0.19), xycoords='axes fraction', fontsize=9)
+
+    if end_time is not None:
+        if nimg is not None:
+            fig.suptitle('Average of %d images between %s - %s ' % (nimg, start_time, end_time), fontsize=9)
+        else:
+            fig.suptitle('Images averaged between %s - %s ' % (start_time, end_time), fontsize=9)
+    else:
+        fig.suptitle('Image at %s' % (start_time), fontsize=9)
+
+    fig.tight_layout()
+
+    fig.savefig(fn_out, dpi=300)