import earthaccess # needed to discover and download TEMPO data
import netCDF4 as nc # needed to read TEMPO data
import os
import sys
import copy
import platform
from subprocess import Popen
import shutil
from shapely.geometry import Point, Polygon # needed to search a point within a polygon
from scipy.interpolate import griddata # needed to interpolate TEMPO data to the point of interest
from scipy import stats # needed for linear regression analysis
import requests # needed to search for and download Pandora data
import codecs # needed to read Pandora data
import numpy as np
import h5py # needed to read DSCOVR_EPIC_L2_TO3 files
import matplotlib.pyplot as plt # needed to plot the resulting time series
from urllib.request import urlopen, Request # needed to search for and download Pandora data
from pathlib import Path # needed to check whether a needed data file is already downloaded
from datetime import datetime, timedelta # needed to work with time in plotting time seriesTEMPO UV Aerosol Index from L2_O3TOT
This notebook illustrates comparison of TEMPO ultra-violet aerosol index (UVAI) against DSCOVR EPIC UVAI as well as comparison of DSCOVR EPIC AOD against AERONET.
Since AERONET data files cannot be downloaded programmatically, a user is supposed to upload an AERONET AOD level 1.5 or 2.0 data file and provide its name as a constant. Since TEMPO spatial coverage is regional and limited to North America, it is user’s responsibilty to select AERONET station within TEMPO’s field of regard (FOR). If the selected station is outside FOR, no TEMPO time series will be generated.
The user is allowed to choose the time period of interest by providing start and end dates in the form YYYYMMDD. Please be aware, that if the selecte period of interest is outside of available time span of one of the sensors, corresponding time series will not be generated.
TEMPO and DSCOVR data are downloaded on-the-fly with earthaccess library.
TEMPO and DSCOVR L2 AER data are interpolated to the location of the selected AERONET station.
1 Importing necessary libraries
2 Defining functions
2.1 function to read DSCOVR AER data files
function read_epic_l2_AER reads DSCOVR_EPIC_L2_AER product file given by its fname and returns arrays of 2D latitudes, longitudes, UVAI, and AOD along with their fill values, time, and wavelengths at which AOD are retrieved
def read_epic_l2_AER(fname):
aod_name = '/HDFEOS/SWATHS/Aerosol NearUV Swath/Data Fields/FinalAerosolOpticalDepth'
uvai_name = '/HDFEOS/SWATHS/Aerosol NearUV Swath/Data Fields/UVAerosolIndex'
lat_name = '/HDFEOS/SWATHS/Aerosol NearUV Swath/Geolocation Fields/Latitude'
lon_name = '/HDFEOS/SWATHS/Aerosol NearUV Swath/Geolocation Fields/Longitude'
wl_name = '/HDFEOS/SWATHS/Aerosol NearUV Swath/Data Fields/Wavelength'
try:
f = h5py.File(fname, "r" )
item = f[aod_name]
aod3D = np.array(item[:])
fv_aod = item.fillvalue
item = f[uvai_name]
uvai2D = np.array(item[:])
fv_uvai = item.fillvalue
item = f[lat_name]
lat2D = np.array(item[:])
fv_lat = item.fillvalue
item = f[lon_name]
lon2D = np.array(item[:])
fv_lon = item.fillvalue
item = f[wl_name]
wl = np.array(item[:])
fv_wl = item.fillvalue
f.close()
# getting time from the granule's filename
fname_split = fname.split('_')
timestamp = fname_split[-2]
yyyy= int(timestamp[0 : 4])
mm = int(timestamp[4 : 6])
dd = int(timestamp[6 : 8])
hh = int(timestamp[8 : 10])
mn = int(timestamp[10 : 12])
ss = int(timestamp[12 : 14])
except:
print("Unable to find or read hdf5 input granule file ", fname)
aod3D = 0.
fv_aod = 0.
uvai2D = 0.
fv_uvai = 0.
lat2D = 0.
fv_lat = 0.
lon2D = 0.
fv_lon = 0.
wl = 0.
fv_wl = 0.
yyyy = 0.
mm = 0.
dd = 0.
hh = 0.
mn = 0.
ss = 0.
return aod3D, fv_aod, uvai2D, fv_uvai, lat2D, fv_lat, lon2D, fv_lon\
, wl, fv_wl, yyyy, mm, dd, hh, mn, ss2.2 function to read AERONET data files
Function read_aeronet_mw takes aeronet file name, a set of wavelengths, and start and end dates of the timeframe of interest as arguments and returns a timeseries array of AODs at these wavelengths.
If the wavelength is not in the list of wavelengths of the data file, or if there are no retrieval at it, empty array is returned.
def read_aeronet_mw(fname, wl, start_date, end_date):
wln = np.array(copy.deepcopy(sorted(wl)))
print(wln)
f_aero = open(fname, 'r')
while True:
line = f_aero.readline()
if not line: break
if 'Date' in line:
header = line.split(',')
nheader = len(header)
break
# finding positions of necessary information in the header
# date and time are always 1st and 2nd, so no need to search for these names in the header
date_pos = 0
time_pos = 1
# wavelengths first
nwl = len(wln)
aod_pos = np.empty([nwl], dtype = int)
aod_pos[:] = -1
for iwl in range(nwl):
aod_str = str('AOD_%dnm' %wln[iwl])
print(aod_str)
for i in range(nheader):
if header[i] == aod_str: aod_pos[iwl] = i
wln = wln[aod_pos >= 0]
nwl = len(wln)
aod_pos = aod_pos[aod_pos >= 0]
print(aod_pos)
print(wln)
if nwl == 0:
print('wavelengths not found')
aod = np.empty([0])
wln = []
date_time = np.empty([0, 2])
lat = -999.
lon = -999.
AERONET_name = ''
return AERONET_name, lat, lon, wln, date_time, aod
# AERONET Site name
name_pos = -1
name_str = 'AERONET_Site_Name'
for i in range(nheader):
if header[i] == name_str: name_pos = i
# AERONET Site latitude
lat_pos = -1
lat_str = 'Site_Latitude(Degrees)'
for i in range(nheader):
if header[i] == lat_str: lat_pos = i
# AERONET Site longitude
lon_pos = -1
lon_str = 'Site_Longitude(Degrees)'
for i in range(nheader):
if header[i] == lon_str: lon_pos = i
aod = np.empty([0, nwl])
date_time = np.empty([0, 2])
# reading the 1st line of data, get name, lat, lon,
# and if aod is valid append arrays aod and date_time
line = f_aero.readline()
values = line.split(',')
AERONET_name = values[name_pos]
lat = float(values[lat_pos])
lon = float(values[lon_pos])
aod_loc = np.empty([nwl])
date_loc = values[date_pos]
time_loc = values[time_pos]
dd = int(date_loc[0:2])
mm = int(date_loc[3:5])
yyyy = int(date_loc[6:10])
date_stamp = yyyy*10000 + mm*100 + dd
if date_stamp >= start_date and date_stamp <= end_date:
for iwl in range(nwl): aod_loc[iwl] = float(values[aod_pos[iwl]])
aod = np.append(aod, [aod_loc], axis = 0)
date_time = np.append(date_time, [[date_loc, time_loc]], axis = 0)
# reading all other lines of data,
# if aod is valid append arrays aod and date_time
while True:
line = f_aero.readline()
if not line: break
values = line.split(',')
aod_loc = np.empty([nwl])
date_loc = values[date_pos]
time_loc = values[time_pos]
dd = int(date_loc[0:2])
mm = int(date_loc[3:5])
yyyy = int(date_loc[6:10])
date_stamp = yyyy*10000 + mm*100 + dd
if date_stamp < start_date or date_stamp > end_date: continue
for iwl in range(nwl): aod_loc[iwl] = float(values[aod_pos[iwl]])
aod = np.append(aod, [aod_loc], axis = 0)
date_time = np.append(date_time, [[date_loc, time_loc]], axis = 0)
f_aero.close()
return AERONET_name, lat, lon, wln, date_time, aod2.3 functions to work with TEMPO
2.3.1 function to read UVAI from TEMPO_O3TOT_L2
function read_TEMPO_O3TOT_L2_UVAI reads the following arrays from the TEMPO L2 O3TOT product TEMPO_O3TOT_L2_V01(2):
uv_aerosol_index;
quality_flag;
and returns respective fields along with their fill values and coordinates of the pixels.
If one requested variables cannot be read, all returned variables are zeroed.
def read_TEMPO_O3TOT_L2_UVAI(fn):
var_name = 'uv_aerosol_index'
var_QF_name = 'quality_flag'
try:
ds = nc.Dataset(fn)
prod = ds.groups['product'] # this opens group product, /product, as prod
var = prod.variables[var_name] # this reads variable column_amount_o3 from prod (group product, /product)
uvai = np.array(var)
uvai_fv = var.getncattr('_FillValue')
var_QF = prod.variables[var_QF_name] # this reads variable column_amount_o3 from prod (group product, /product)
uvai_QF = np.array(var_QF)
# there is no fill value for the quality flag.
# Once it is available in the next version of the product,
# un-comment the line below and add fv_QF to the return line.
# fv_QF = var_QF.getncattr('_FillValue')
geo = ds.groups['geolocation'] # this opens group geolocation, /geolocation, as geo
lat = np.array(geo.variables['latitude']) # this reads variable latitude from geo (geolocation group, /geolocation) into a numpy array
lon = np.array(geo.variables['longitude']) # this reads variable longitude from geo (geolocation group, /geolocation) into a numpy array
fv_geo = geo.variables['latitude'].getncattr('_FillValue')
# it appeared that garbage values of latitudes and longitudes in the L2 files
# are 9.969209968386869E36 while fill value is -1.2676506E30
# (after deeper search it was found that actual value in the file is -1.2676506002282294E30).
# For this reason, fv_geo is set to 9.96921E36 to make the code working.
# Once the problem is resolved and garbage values of latitudes and longitudes
# equal to their fill value, the line below must be removed.
fv_geo = 9.969209968386869E36
time = np.array(geo.variables['time'] )# this reads variable longitude from geo (geolocation group, /geolocation) into a numpy array
ds.close()
except:
print('variable '+var_name+' cannot be read in file '+fn)
lat = 0.
lon = 0.
time = 0.
fv_geo = 0.
uvai = 0.
uvai_QF = 0.
fv_uvai = 0.
# fv_QF = -999
prod_unit = ''
return lat, lon, fv_geo, time, uvai, uvai_QF, uvai_fv2.3.2 function creating TEMPO O3 granule polygon
def TEMPO_L2_polygon(lat, lon, fv_geo):
nx = lon.shape[0]
ny = lon.shape[1]
print('granule has %3d scanlines by %4d pixels' %(nx, ny))
dpos = np.empty([0,2])
x_ind = np.empty([nx, ny], dtype = int) # creating array in x indices
for ix in range(nx): x_ind[ix, :] = ix # populating array in x indices
y_ind = np.empty([nx, ny], dtype = int)
for iy in range(ny): y_ind[:, iy] = iy # populating array in x indices
mask = (lon[ix, iy] != fv_geo)&(lat[ix, iy] != fv_geo)
if len(lon[mask]) == 0:
print('the granule is empty - no meaningful positions')
return dpos
# right boundary
r_m = min(x_ind[mask].flatten())
local_mask = (lon[r_m, :] != fv_geo)&(lat[r_m, :] != fv_geo)
r_b = np.stack((lon[r_m, local_mask], lat[r_m, local_mask])).T
# left boundary
l_m = max(x_ind[mask].flatten())
local_mask = (lon[l_m, :] != fv_geo)&(lat[l_m, :] != fv_geo)
l_b = np.stack((lon[l_m, local_mask], lat[l_m, local_mask])).T
#top and bottom boundaries
t_b = np.empty([0,2])
b_b = np.empty([0,2])
for ix in range(r_m + 1, l_m):
local_mask = (lon[ix, :] != fv_geo)&(lat[ix, :] != fv_geo)
local_y_ind = y_ind[ix, local_mask]
y_ind_top = min(local_y_ind)
y_ind_bottom = max(local_y_ind)
t_b = np.append(t_b, [[lon[ix, y_ind_top], lat[ix, y_ind_top]]], axis=0)
b_b = np.append(b_b, [[lon[ix, y_ind_bottom], lat[ix, y_ind_bottom]]], axis=0)
# combining right, top, left, and bottom boundaries together, going along the combined boundary counterclockwise
dpos = np.append(dpos, r_b[ : :-1, :], axis=0) # this adds right boundary, counterclockwise
dpos = np.append(dpos, t_b, axis=0) # this adds top boundary, counterclockwise
dpos = np.append(dpos, l_b, axis=0) # this adds left boundary, counterclockwise
dpos = np.append(dpos, b_b[ : :-1, :], axis=0) # this adds bottom boundary, counterclockwise
print('polygon shape: ',dpos.shape)
return dposMain code begins here
3 Establishing access to EarthData
User needs to create an account at https://www.earthdata.nasa.gov/ Function earthaccess.login prompts for EarthData login and password.
# logging in
auth = earthaccess.login(strategy="interactive", persist=True)
# Creating local directory
homeDir = os.path.expanduser("~") + os.sep
with open(homeDir + '.dodsrc', 'w') as file:
file.write('HTTP.COOKIEJAR={}.urs_cookies\n'.format(homeDir))
file.write('HTTP.NETRC={}.netrc'.format(homeDir))
file.close()
print('Saved .dodsrc to:', homeDir)
# Set appropriate permissions for Linux/macOS
if platform.system() != "Windows":
Popen('chmod og-rw ~/.netrc', shell=True)
else:
# Copy dodsrc to working directory in Windows
shutil.copy2(homeDir + '.dodsrc', os.getcwd())
print('Copied .dodsrc to:', os.getcwd())Saved .dodsrc to: /home/jovyan/4 Select timeframe of interest
the timeframe is going to be common for all instruments
print('enter period of interest, start and end dates, in the form YYYYMMDD')
datestamp_ini = input('enter start date of interest ')
datestamp_fin = input('enter end date of interest ')
start_date = int(datestamp_ini)
end_date = int(datestamp_fin)
yyyy_ini = start_date//10000
mm_ini = (start_date//100 - yyyy_ini*100)
dd_ini = (start_date - yyyy_ini*10000 - mm_ini*100)
yyyy_fin = end_date//10000
mm_fin = (end_date//100 - yyyy_fin*100)
dd_fin = (end_date - yyyy_fin*10000 - mm_fin*100)
print(yyyy_ini, mm_ini, dd_ini, yyyy_fin, mm_fin, dd_fin)
date_start = str('%4.4i-%2.2i-%2.2i 00:00:00' %(yyyy_ini, mm_ini, dd_ini))
date_end = str('%4.4i-%2.2i-%2.2i 23:59:59' %(yyyy_fin, mm_fin, dd_fin))enter period of interest, start and end dates, in the form YYYYMMDDenter start date of interest 20230805
enter end date of interest 202308052023 8 5 2023 8 54.1 calculate datetime objects for the beginning and end of the timeframe of interest
using datetime library
yyyy = yyyy_ini
mm = mm_ini
dd = dd_ini
hh = 0
mn = 0
ss = 0
us = 0 # microseconds
dt_ini = datetime(yyyy, mm, dd, hh, mn, ss, us)
yyyy = yyyy_fin
mm = mm_fin
dd = dd_fin
hh = 23
mn = 59
ss = 59
us = 999999 # microseconds
dt_fin = datetime(yyyy, mm, dd, hh, mn, ss, us) # this is time 1 second before the end of the timeframe of interest5 work with AERONET data
As of now, April 24, 2024, AERONET data cannot be downloaded programmatically. User must upload an AERONET data file using AERONET download tool here https://aeronet.gsfc.nasa.gov/new_web/webtool_aod_v3.html. The first line below is the name of the AERONET file used in this example.
fn = '20230801_20241231_CCNY.lev15'
wl = [500, 340, 380]5.1 read AERONET data
AERONET_name, lat, lon, wln, date_time, aod = read_aeronet_mw(fn, wl, start_date, end_date)
POI_lat = lat
POI_lon = lon
POI_name = AERONET_name
nt = len(date_time)
print(wln)
print(len(aod))
print('name %s, latitude %9.4f, longitude %10.4f' %(AERONET_name, lat, lon))
nwl = len(wln)
for iwl in range(nwl): print(wln[iwl], len(aod[aod[:,iwl]>0, iwl]))[340 380 500]
AOD_340nm
AOD_380nm
AOD_500nm
[25 24 18]
[340 380 500]
[340 380 500]
22
name CCNY, latitude 40.8213, longitude -73.9490
340 22
380 22
500 225.2 create timeseries of AERONET AODs
calculate time difference, dy_loc, in fractional days between AERONET observations and the beginning of the timeframe of interest, write date and time as yyyy, mm, dd, hh, mn, ss, dt_loc along with wavelengths and corresponding AODs to a file
out_Q_AERONET = 'AOD_AERONET'
fout = open(out_Q_AERONET+'_'+datestamp_ini+'_'+datestamp_fin+'_'\
+POI_name+'_'+str('%08.4fN_%08.4fW.txt' %(POI_lat, -POI_lon)), 'w')
fout.write('timeseries of '+out_Q_AERONET+' at '+POI_name+' '+str('%08.4fN %08.4fW' %(POI_lat, -POI_lon))+'\n')
fout.write('yyyy mm dd hh mn ss time,days wl,nm AOD wl,nm AOD wl,nm AOD\n')
for it in range(nt):
d_t = date_time[it]
yyyy = int(d_t[0][6:10])
mm = int(d_t[0][3:5])
dd = int(d_t[0][0:2])
hh = int(d_t[1][0:2])
mn = int(d_t[1][3:5])
ss = int(d_t[1][6:8])
us = 0 # microseconds
dt_AERONET = datetime(yyyy, mm, dd, hh, mn, ss, us)
dt_loc = (dt_AERONET - dt_ini).total_seconds()/86400
fout.write('%4i %2.2i %2.2i %2.2i %2.2i %2.2i %9.6f' %(yyyy, mm, dd, hh, mn, ss, dt_loc))
for wl, tau in zip(wln, aod[it]): fout.write(' %5.0f %6.3f' %(wl, max(tau, -1.)))
fout.write('\n')
print(d_t, yyyy, mm, dd, hh, mn, ss, dt_loc, aod[it])
fout.close()['05:08:2023' '12:21:54'] 2023 8 5 12 21 54 0.5152083333333334 [0.494809 0.458947 0.333314]
['05:08:2023' '12:37:38'] 2023 8 5 12 37 38 0.5261342592592593 [0.460212 0.425241 0.306128]
['05:08:2023' '12:39:36'] 2023 8 5 12 39 36 0.5275 [0.468725 0.432512 0.311574]
['05:08:2023' '12:41:35'] 2023 8 5 12 41 35 0.5288773148148148 [0.484136 0.453612 0.324803]
['05:08:2023' '12:43:34'] 2023 8 5 12 43 34 0.5302546296296297 [0.488803 0.454636 0.328589]
['05:08:2023' '12:45:32'] 2023 8 5 12 45 32 0.5316203703703704 [0.490196 0.453847 0.335888]
['05:08:2023' '12:54:52'] 2023 8 5 12 54 52 0.5381018518518519 [0.478464 0.44467 0.325635]
['05:08:2023' '13:10:01'] 2023 8 5 13 10 1 0.5486226851851852 [0.494555 0.460215 0.333986]
['05:08:2023' '13:19:13'] 2023 8 5 13 19 13 0.5550115740740741 [0.486189 0.452743 0.330567]
['05:08:2023' '13:35:01'] 2023 8 5 13 35 1 0.5659837962962962 [0.418951 0.386741 0.276587]
['05:08:2023' '13:44:19'] 2023 8 5 13 44 19 0.5724421296296296 [0.404995 0.374061 0.266198]
['05:08:2023' '13:59:33'] 2023 8 5 13 59 33 0.5830208333333333 [0.427243 0.394882 0.280744]
['05:08:2023' '14:09:10'] 2023 8 5 14 9 10 0.5896990740740741 [0.400944 0.369637 0.262378]
['05:08:2023' '14:24:59'] 2023 8 5 14 24 59 0.6006828703703704 [0.402717 0.373348 0.269201]
['05:08:2023' '14:29:02'] 2023 8 5 14 29 2 0.6034953703703704 [0.376444 0.346455 0.245482]
['05:08:2023' '15:32:29'] 2023 8 5 15 32 29 0.6475578703703704 [0.326602 0.30022 0.208331]
['05:08:2023' '15:42:29'] 2023 8 5 15 42 29 0.6545023148148148 [0.35099 0.325715 0.23699 ]
['05:08:2023' '15:47:29'] 2023 8 5 15 47 29 0.657974537037037 [0.328565 0.302978 0.216409]
['05:08:2023' '16:02:30'] 2023 8 5 16 2 30 0.6684027777777778 [0.334116 0.308406 0.220808]
['05:08:2023' '16:04:48'] 2023 8 5 16 4 48 0.67 [0.334563 0.308769 0.221642]
['05:08:2023' '16:14:27'] 2023 8 5 16 14 27 0.6767013888888889 [0.405008 0.380736 0.294151]
['05:08:2023' '16:30:17'] 2023 8 5 16 30 17 0.6876967592592592 [0.356773 0.32694 0.246413]6 work with DSCOVR data
6.1 search for DSCOVR EPIC granules
falling into the timeframe of interest and 1 degree box surrounding the AERONET location
short_name = 'DSCOVR_EPIC_L2_AER' # collection name to search for in the EarthData
bbox = (POI_lon - 0.5, POI_lat - 0.5, POI_lon + 0.5, POI_lat + 0.5)
POI_results_EPIC = earthaccess.search_data(short_name = short_name\
, temporal = (date_start, date_end)\
, bounding_box = bbox)
n_EPIC = len(POI_results_EPIC)
print('total number of DSCOVR EPIC L2_AER granules found for POI', POI_name\
, '\nwithin period of interes between', date_start, 'and', date_end, '\nis', n_EPIC)Granules found: 11
total number of DSCOVR EPIC L2_AER granules found for POI CCNY
within period of interes between 2023-08-05 00:00:00 and 2023-08-05 23:59:59
is 116.2 ensure all granules have download links
lines below before the call of earthaccess.download() check whether all found granules have download links. granules without links are removed from the list of search results without this step, those granules crash the call of earthaccess.download()
granule_links_EPIC = []
POI_results_EPIC_bad = []
for result in POI_results_EPIC:
try:
granule_links_EPIC.append(result['umm']['RelatedUrls'][0]['URL'])
except:
POI_results_EPIC_bad.append(result)
for granule_link in granule_links_EPIC: print(granule_link)
for result in POI_results_EPIC_bad: POI_results_EPIC.remove(result)https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805135123_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805114028_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805223501_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805191839_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805170745_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805145650_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805181312_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805212934_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805160217_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805202406_03.he5
https://asdc.larc.nasa.gov/data/DSCOVR/EPIC/L2_AER_03/2023/08/DSCOVR_EPIC_L2_AER_03_20230805124555_03.he56.3 download DSCOVR EPIC granules using earthaccess library
check whether all downloads were successful and try another time to get the failed granules if the second attempt was not successful, remove failed granules from the list
downloaded_files = earthaccess.download(POI_results_EPIC, local_path='.',)
# Checking whether all DSCOVR EPIC data files have been downloaded
for granule_link in granule_links_EPIC:
EPIC_fname = granule_link.split('/')[-1]
# check if file exists in the local directory
if not os.path.exists(EPIC_fname):
print(EPIC_fname, 'does not exist in local directory')
# repeat attempt to download
downloaded_files = earthaccess.download(granule_link,
local_path='.')
# if file still does not exist in the directory, remove its link from the list of links
if not os.path.exists(EPIC_fname): granule_links_EPIC.remove(granule_link) Getting 11 granules, approx download size: 0.0 GBFile DSCOVR_EPIC_L2_AER_03_20230805135123_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805223501_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805114028_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805170745_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805145650_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805191839_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805181312_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805212934_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805160217_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805202406_03.he5 already downloaded
File DSCOVR_EPIC_L2_AER_03_20230805124555_03.he5 already downloaded6.4 read DSCOVR EPIC granules and compile timeseries of AODs and UVAI
# parameter geo_deviation defines maximum difference between lats and lons of EPIC pixels and that of AERONET
geo_deviation = 0.075
pp = np.array([POI_lon, POI_lat])
out_Q_EPIC = 'AOD_UVAI_EPIC'
fout = open(out_Q_EPIC+'_'+datestamp_ini+'_'+datestamp_fin+'_'\
+POI_name+'_'+str('%08.4fN_%08.4fW.txt' %(POI_lat, -POI_lon)), 'w')
fout.write('timeseries of '+out_Q_EPIC+' at '+POI_name+' '+str('%08.4fN %08.4fW' %(POI_lat, -POI_lon))+'\n')
fout.write('yyyy mm dd hh mn ss time,days wl,nm AOD wl,nm AOD wl,nm AOD UVAI\n')
for granule_link in sorted(granule_links_EPIC):
EPIC_fname = granule_link.split('/')[-1]
print(EPIC_fname)
aod3D, fv_aod, uvai2D, fv_uvai, lat2D, fv_lat, lon2D, fv_lon\
, wl, fv_wl, yyyy, mm, dd, hh, mn, ss = read_epic_l2_AER(EPIC_fname)
dt_EPIC = datetime(yyyy, mm, dd, hh, mn, ss)
dt_loc = (dt_EPIC - dt_ini).total_seconds()/86400
print(yyyy, mm, dd, hh, mn, ss, (dt_EPIC - dt_ini).total_seconds()/86400)
fout.write('%4i %2.2i %2.2i %2.2i %2.2i %2.2i %9.6f' %(yyyy, mm, dd, hh, mn, ss, dt_loc))
# check whether POI is in the granule. If not - move to the next granule
mask_geo = (lat2D < POI_lat+geo_deviation)&(lat2D > POI_lat-geo_deviation)\
&(lon2D < POI_lon+geo_deviation)&(lon2D > POI_lon-geo_deviation)
nwl_EPIC = len(wl)
for iwl in range(nwl_EPIC):
mask = mask_geo&(aod3D[:, :, iwl] > 0.)
lat_loc = lat2D[mask]
lon_loc = lon2D[mask]
aod_loc = aod3D[mask, iwl]
n_loc = len(aod_loc)
# if n_loc < 1: continue
if n_loc < 1: prod_loc = -0.999
else:
points = np.empty([0,2])
ff = np.empty(0)
for i in range(n_loc):
points = np.append(points, [[lon_loc[i], lat_loc[i]]], axis=0)
ff = np.append(ff, aod_loc[i])
try:
[prod_loc] = griddata(points, ff, pp, method='linear', fill_value=-0.999, rescale=False)
except:
try:
prod_loc = np.mean(ff)
except: prod_loc = -0.999
# if prod_loc == -999.: continue
fout.write(' %5.0f %6.3f' %(wl[iwl], prod_loc))
print(wl[iwl], prod_loc)
mask = mask_geo&(uvai2D != fv_uvai)
lat_loc = lat2D[mask]
lon_loc = lon2D[mask]
uvai_loc = uvai2D[mask]
n_loc = len(aod_loc)
# if n_loc < 1: continue
if n_loc < 1: prod_loc = -99.
else:
points = np.empty([0,2])
ff = np.empty(0)
for i in range(n_loc):
points = np.append(points, [[lon_loc[i], lat_loc[i]]], axis=0)
ff = np.append(ff, uvai_loc[i])
try:
[prod_loc] = griddata(points, ff, pp, method='linear', fill_value=-99., rescale=False)
except:
try:
prod_loc = np.mean(ff)
except: prod_loc = -99.
fout.write(' %6.2f\n' %prod_loc)
print('UVAI', prod_loc)
fout.close()DSCOVR_EPIC_L2_AER_03_20230805114028_03.he5
2023 8 5 11 40 28 0.4864351851851852
340.0 -0.999
388.0 -0.999
500.0 -0.999
UVAI -99.0
DSCOVR_EPIC_L2_AER_03_20230805124555_03.he5
2023 8 5 12 45 55 0.531886574074074
340.0 -0.999
388.0 -0.999
500.0 -0.999
UVAI -99.0
DSCOVR_EPIC_L2_AER_03_20230805135123_03.he5
2023 8 5 13 51 23 0.5773495370370371
340.0 0.658370852470398
388.0 0.6070993542671204
500.0 0.5208560228347778
UVAI 1.0490291118621826
DSCOVR_EPIC_L2_AER_03_20230805145650_03.he5
2023 8 5 14 56 50 0.622800925925926
340.0 1.0037826597690582
388.0 0.925611823797226
500.0 0.7941211760044098
UVAI 1.269827425479889
DSCOVR_EPIC_L2_AER_03_20230805160217_03.he5
2023 8 5 16 2 17 0.6682523148148148
340.0 -0.999
388.0 -0.999
500.0 -0.999
UVAI -99.0
DSCOVR_EPIC_L2_AER_03_20230805170745_03.he5
2023 8 5 17 7 45 0.7137152777777778
340.0 2.0829079262483576
388.0 1.6970317278376417
500.0 1.0809463897588194
UVAI -0.14213168194983128
DSCOVR_EPIC_L2_AER_03_20230805181312_03.he5
2023 8 5 18 13 12 0.7591666666666667
340.0 1.6575589177292698
388.0 1.528474366786192
500.0 1.3113422907254944
UVAI 0.9618683264479624
DSCOVR_EPIC_L2_AER_03_20230805191839_03.he5
2023 8 5 19 18 39 0.8046180555555555
340.0 1.46158707100744
388.0 1.347764097820558
500.0 1.1563033843281367
UVAI 1.5943157504218373
DSCOVR_EPIC_L2_AER_03_20230805202406_03.he5
2023 8 5 20 24 6 0.8500694444444444
340.0 -0.999
388.0 -0.999
500.0 -0.999
UVAI -99.0
DSCOVR_EPIC_L2_AER_03_20230805212934_03.he5
2023 8 5 21 29 34 0.8955324074074074
340.0 4.511193513870239
388.0 4.159878492355347
500.0 3.568934202194214
UVAI 3.3587218523025513
DSCOVR_EPIC_L2_AER_03_20230805223501_03.he5
2023 8 5 22 35 1 0.9409837962962962
340.0 -0.999
388.0 -0.999
500.0 -0.999
UVAI -99.07 Work with TEMPO data
7.1 Search for TEMPO data files within 0.5 degree range around the POI
(position of the AERONET station) falling into the timeframe of interest
# Setting TEMPO name constants
short_name = 'TEMPO_O3TOT_L2' # collection name to search for in the EarthData
out_Q = 'UVAI_TEMPO' # name of the output quantity with unit
version = 'V03'
# Searching TEMPO data files within 0.5 degree range around the POI
bbox = (POI_lon - 0.5, POI_lat - 0.5, POI_lon + 0.5, POI_lat + 0.5)
POI_results = earthaccess.search_data(short_name = short_name\
, version = version\
, temporal = (date_start, date_end)\
, bounding_box = bbox)
n_gr = len(POI_results)
print('total number of TEMPO version ', version,' granules found for POI', POI_name, \
'\nwithin period of interes between', date_start, 'and', date_end, ' is', n_gr)
if n_gr == 0:
print('program terminated')
sys.exit()
# Print links to the granules.
granule_links = []
for result in POI_results: granule_links.append(result['umm']['RelatedUrls'][0]['URL'])
for granule_link in granule_links: print(granule_link)Granules found: 11
total number of TEMPO version V03 granules found for POI CCNY
within period of interes between 2023-08-05 00:00:00 and 2023-08-05 23:59:59 is 11
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T122445Z_S001G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T132716Z_S002G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T142947Z_S003G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T153218Z_S004G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T163449Z_S005G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T173720Z_S006G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T183951Z_S007G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T194222Z_S008G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T204453Z_S009G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T214724Z_S010G03.nc
https://data.asdc.earthdata.nasa.gov/asdc-prod-protected/TEMPO/TEMPO_O3TOT_L2_V03/2023.08.05/TEMPO_O3TOT_L2_V03_20230805T224955Z_S011G03.nc7.2 download TEMPO files
check whether all downloads were successful, try to download failed granules yet another time. if second attempt fails, remove those granules from the list.
# Downloading TEMPO data files
downloaded_files = earthaccess.download(
POI_results,
local_path='.')
# Checking whether all TEMPO data files have been downloaded
for granule_link in granule_links:
TEMPO_fname = granule_link.split('/')[-1]
# check if file exists in the local directory
if not os.path.exists(TEMPO_fname):
print(TEMPO_fname, 'does not exist in local directory')
# repeat attempt to download
downloaded_files = earthaccess.download(granule_link,
local_path='.')
# if file still does not exist in the directory, remove its link from the list of links
if not os.path.exists(TEMPO_fname): granule_links.remove(granule_link) Getting 11 granules, approx download size: 1.14 GB
Accessing cloud dataset using dataset endpoint credentials: https://data.asdc.earthdata.nasa.gov/s3credentials
Downloaded: TEMPO_O3TOT_L2_V03_20230805T122445Z_S001G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T132716Z_S002G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T142947Z_S003G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T153218Z_S004G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T163449Z_S005G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T173720Z_S006G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T183951Z_S007G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T194222Z_S008G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T204453Z_S009G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T214724Z_S010G03.nc
Downloaded: TEMPO_O3TOT_L2_V03_20230805T224955Z_S011G03.nc7.3 compile TEMPO UVAI timeseries
days = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
fout = open(out_Q+'_'+datestamp_ini+'_'+datestamp_fin+'_'\
+POI_name+'_'+str('%08.4fN_%08.4fW.txt' %(POI_lat, -POI_lon)), 'w')
fout.write('timeseries of '+out_Q+' at '+POI_name+' '+str('%08.4fN %08.4fW' %(POI_lat, -POI_lon))+'\n')
fout.write('yyyy mm dd hh mn ss time,days UVAI\n')
pp = np.array([POI_lon, POI_lat])
p = Point(pp)
for granule_link in sorted(granule_links):
last_slash_ind = granule_link.rfind('/')
fname = granule_link[last_slash_ind+1 : ]
print('\ngranule ', fname)
lat, lon, fv_geo, time, uvai, uvai_QF, uvai_fv = read_TEMPO_O3TOT_L2_UVAI(fname)
if isinstance(lat, float): continue
polygon = TEMPO_L2_polygon(lat, lon, fv_geo) # create granule polygon
coords_poly = list(polygon)
poly = Polygon(coords_poly)
nx = lon.shape[0]
ny = lon.shape[1]
# getting time from the granule filename.
# This time corresponds to the 1st element of array time above; that is GPS time in seconds
Tind = fname.rfind('T')
yyyy= int(fname[Tind-8 : Tind-4])
mm = int(fname[Tind-4 : Tind-2])
dd = int(fname[Tind-2 : Tind])
hh = int(fname[Tind+1 : Tind+3])
mn = int(fname[Tind+3 : Tind+5])
ss = float(fname[Tind+5 : Tind+7])
# check whether POI is in the granule. If not - move to the next granule
if not p.within(poly): continue
print('point', p, 'is in granule polygon' )
POI_found = False
for ix in range(nx-1):
for iy in range(ny-1):
if lon[ix, iy] == fv_geo: continue
if lat[ix, iy] == fv_geo: continue
if lon[ix, iy+1] == fv_geo: continue
if lat[ix, iy+1] == fv_geo: continue
if lon[ix+1, iy+1] == fv_geo: continue
if lat[ix+1, iy+1] == fv_geo: continue
if lon[ix+1, iy] == fv_geo: continue
if lat[ix+1, iy] == fv_geo: continue
coords_poly_loc = [[lon[ix, iy], lat[ix, iy]], [lon[ix, iy+1], lat[ix, iy+1]] \
, [lon[ix+1, iy+1], lat[ix+1, iy+1]], [lon[ix+1, iy], lat[ix+1, iy]]]
poly_loc = Polygon(coords_poly_loc)
if p.within(poly_loc):
print('scanl pixel latitude longitude UVAI UVAI_QF')
for scl in range(ix, ix+2, 1):
for pix in range(iy, iy+2, 1):
print(" %3d %4d %9.6f %10.6f %5.1f %6i"\
%(scl, pix, lat[scl, pix], lon[scl, pix]\
, uvai[scl, pix], uvai_QF[scl, pix]))
POI_found = True
print('POI', POI_name, 'at', POI_lat,'N ', POI_lon, 'E found')
uvai_loc = np.array([uvai[ix, iy],\
uvai[ix, iy+1],\
uvai[ix+1, iy+1],\
uvai[ix+1, iy]])
uvai_QF_loc = np.array([uvai_QF[ix, iy],\
uvai_QF[ix, iy+1],\
uvai_QF[ix+1, iy+1],\
uvai_QF[ix+1, iy]])
lat_loc = np.array([lat[ix, iy], lat[ix, iy+1],\
lat[ix+1, iy+1], lat[ix+1, iy]])
lon_loc = np.array([lon[ix, iy], lon[ix, iy+1],\
lon[ix+1, iy+1], lon[ix+1, iy]])
# mask_noFV = (uvai_loc != uvai_fv)&\
# (uvai_QF_loc == 0)
mask_noFV = (uvai_loc != uvai_fv)
points_noFV = np.column_stack((lon_loc[mask_noFV], lat_loc[mask_noFV]))
ff_noFV = uvai_loc[mask_noFV]
points = np.empty([0,2])
ff = np.empty(0)
if ff_noFV.shape[0] == 0:
continue
elif ff_noFV.shape[0] < 4:
uvai_noFV = np.mean(ff_noFV)
elif ff_noFV.shape[0] == 4:
uvai_noFV = griddata(points_noFV, ff_noFV, pp,\
method='linear', fill_value=-99., rescale=False)[0]
if uvai_noFV == -99.: continue
# handling time first:
delta_t = (time[ix+1] + time[ix])*0.5 - time[0]
ss = ss + delta_t
if ss >= 60.:
delta_mn = int(ss/60.)
ss = ss - 60.*delta_mn
mn = mn + delta_mn
if mn >= 60:
mn = mn - 60
hh = hh + 1
if hh == 24:
hh = hh - 24
dd = dd + 1
day_month = days[mm]
if (yyyy//4)*4 == yyyy and mm == 2: day_month = day_month + 1
if dd > day_month:
dd = 1
mm = mm + 1
if mm > 12:
mm = 1
yyyy = yyyy + 1
int_ss = int(ss)
us = int((ss - int_ss)*10**6)
dt_TEMPO = datetime(yyyy, mm, dd, hh, mn, int_ss, us)
dt_loc = (dt_TEMPO - dt_ini).total_seconds()/86400
fout.write(str('%4.4i %2.2i %2.2i %2.2i %2.2i %2.2i %9.6f %10.3e '\
%(yyyy, mm, dd, hh, mn, ss, dt_loc, uvai_noFV)))
fout.write(str('%9.4fN %9.4fW %10.3e '\
%(lat[ix, iy], -lon[ix, iy], uvai[ix, iy])))
fout.write(str('%9.4fN %9.4fW %10.3e '\
%(lat[ix, iy+1], -lon[ix, iy+1], uvai[ix, iy+1])))
fout.write(str('%9.4fN %9.4fW %10.3e '\
%(lat[ix+1, iy+1], -lon[ix+1, iy+1], uvai[ix+1, iy+1])))
fout.write(str('%9.4fN %9.4fW %10.3e\n'\
%(lat[ix+1, iy], -lon[ix+1, iy], uvai[ix+1, iy])))
break
if POI_found: break
fout.close()
granule TEMPO_O3TOT_L2_V03_20230805T122445Z_S001G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
59 661 40.834648 -73.943382 -0.9 8
59 662 40.813667 -73.949684 -0.8 8
60 661 40.831947 -74.001534 -0.8 8
60 662 40.810905 -74.007851 -0.8 8
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T132716Z_S002G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
59 660 40.834454 -73.943253 -1.8 13
59 661 40.813469 -73.949562 -1.4 13
60 660 40.833084 -74.000923 -1.8 13
60 661 40.812099 -74.007202 -1.5 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T142947Z_S003G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
58 660 40.837505 -73.885612 -1.5 13
58 661 40.816528 -73.891937 -1.3 13
59 660 40.836674 -73.944603 -2.4 13
59 661 40.815800 -73.950844 -2.3 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T153218Z_S004G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
58 660 40.821766 -73.911758 -2.4 13
58 661 40.800797 -73.918060 -2.4 13
59 660 40.821194 -73.968063 -2.6 13
59 661 40.800110 -73.974419 -2.5 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T163449Z_S005G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
58 660 40.824287 -73.898758 -1.1 13
58 661 40.803333 -73.905052 -1.5 13
59 660 40.823364 -73.957039 -1.6 13
59 661 40.802170 -73.963463 -1.7 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T173720Z_S006G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
58 660 40.831722 -73.893791 -1.5 13
58 661 40.810757 -73.900108 -1.5 13
59 660 40.829990 -73.952919 -2.1 13
59 661 40.808990 -73.959221 -1.9 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T183951Z_S007G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
58 659 40.843174 -73.891533 -1.4 13
58 660 40.822109 -73.897896 -1.5 13
59 659 40.841095 -73.952232 -1.2 13
59 660 40.820091 -73.958534 -1.6 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T194222Z_S008G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
59 660 40.837013 -73.913071 -1.0 8
59 661 40.815960 -73.919426 -0.9 13
60 660 40.835320 -73.974220 -1.1 8
60 661 40.814457 -73.980431 -0.8 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T204453Z_S009G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
59 660 40.825069 -73.947708 -1.1 8
59 661 40.803844 -73.954155 -0.8 8
60 660 40.823502 -74.007164 -0.9 8
60 661 40.802341 -74.013550 -0.8 13
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T214724Z_S010G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
59 660 40.833599 -73.914108 -1.0 8
59 661 40.812668 -73.920380 -1.1 8
60 660 40.830616 -73.972809 -1.2 8
60 661 40.810017 -73.978859 -0.9 8
POI CCNY at 40.821313 N -73.949036 E found
granule TEMPO_O3TOT_L2_V03_20230805T224955Z_S011G03.nc
granule has 123 scanlines by 2048 pixels
polygon shape: (4338, 2)
point POINT (-73.949036 40.821313) is in granule polygon
scanl pixel latitude longitude UVAI UVAI_QF
59 661 40.826336 -73.932693 -1.0 8
59 662 40.805405 -73.938957 -1.0 8
60 661 40.823360 -73.993324 -1.2 8
60 662 40.802219 -73.999695 -1.3 8
POI CCNY at 40.821313 N -73.949036 E found8 reading back timeseries written in the files and plot the result
8.1 AERONET
fin = open(out_Q_AERONET+'_'+datestamp_ini+'_'+datestamp_fin+'_'\
+POI_name+'_'+str('%08.4fN_%08.4fW.txt' %(POI_lat, -POI_lon)), 'r')
data_lines = fin.readlines()
fin.close()
timeseries_AERONET_340 = np.empty([0, 2])
timeseries_AERONET_380 = np.empty([0, 2])
timeseries_AERONET_500 = np.empty([0, 2])
for data_line in data_lines[2 : ]:
split = data_line.split()
tt = float(split[6])
aod_340 = float(split[8])
aod_380 = float(split[10])
aod_500 = float(split[12])
if aod_340 > 0.: timeseries_AERONET_340 = np.append(timeseries_AERONET_340, [[tt, aod_340]], axis = 0)
if aod_380 > 0.: timeseries_AERONET_380 = np.append(timeseries_AERONET_380, [[tt, aod_380]], axis = 0)
if aod_500 > 0.: timeseries_AERONET_500 = np.append(timeseries_AERONET_500, [[tt, aod_500]], axis = 0)8.2 DSCOVR EPIC
fin = open(out_Q_EPIC+'_'+datestamp_ini+'_'+datestamp_fin+'_'\
+POI_name+'_'+str('%08.4fN_%08.4fW.txt' %(POI_lat, -POI_lon)), 'r')
data_lines = fin.readlines()
fin.close()
timeseries_EPIC_340 = np.empty([0, 2])
timeseries_EPIC_388 = np.empty([0, 2])
timeseries_EPIC_500 = np.empty([0, 2])
timeseries_EPIC_UVAI = np.empty([0, 2])
for data_line in data_lines[2 : ]:
split = data_line.split()
tt = float(split[6])
aod_EPIC_340 = float(split[8])
aod_EPIC_388 = float(split[10])
aod_EPIC_500 = float(split[12])
EPIC_UVAI = float(split[13])
if aod_EPIC_340 > 0.: timeseries_EPIC_340 = np.append(timeseries_EPIC_340, [[tt, aod_EPIC_340]], axis = 0)
if aod_EPIC_388 > 0.: timeseries_EPIC_388 = np.append(timeseries_EPIC_388, [[tt, aod_EPIC_388]], axis = 0)
if aod_EPIC_500 > 0.: timeseries_EPIC_500 = np.append(timeseries_EPIC_500, [[tt, aod_EPIC_500]], axis = 0)
if EPIC_UVAI > -99.: timeseries_EPIC_UVAI = np.append(timeseries_EPIC_UVAI, [[tt, EPIC_UVAI]], axis = 0)# third, TEMPO UVAI8.3 TEMPO
fin = open(out_Q+'_'+datestamp_ini+'_'+datestamp_fin+'_'\
+POI_name+'_'+str('%08.4fN_%08.4fW.txt' %(POI_lat, -POI_lon)), 'r')
data_lines = fin.readlines()
fin.close()
timeseries_TEMPO_UVAI = np.empty([0, 2])
for data_line in data_lines[2 : ]:
split = data_line.split()
tt = float(split[6])
TEMPO_UVAI = float(split[7])
if TEMPO_UVAI > -99.: timeseries_TEMPO_UVAI = np.append(timeseries_TEMPO_UVAI, [[tt, TEMPO_UVAI]], axis = 0)8.4 plot Aerosol Optical Depth (AOD) from AERONET and DSCOVR EPIC
plot_title = 'AOD_'+datestamp_ini+'_'+datestamp_fin+'\n'+POI_name
img_name = 'AOD_'+datestamp_ini+'_'+datestamp_fin+'_'+POI_name+'.jpg'
plt.plot(timeseries_AERONET_340[:, 0], timeseries_AERONET_340[:, 1],\
label = "AOD AERONET at 340nm", linestyle='None', markersize=2,\
marker='s', markerfacecolor='r', markeredgecolor='r')
plt.plot(timeseries_AERONET_380[:, 0], timeseries_AERONET_380[:, 1],\
label = "AOD AERONET at 380nm", linestyle='None', markersize=2,\
marker='o', markerfacecolor='r', markeredgecolor='r')
plt.plot(timeseries_AERONET_500[:, 0], timeseries_AERONET_500[:, 1],\
label = "AOD AERONET at 500nm", linestyle='None', markersize=2,\
marker='D', markerfacecolor='r', markeredgecolor='r')
plt.plot(timeseries_EPIC_340[:, 0], timeseries_EPIC_340[:, 1],\
label = "AOD EPIC at 340nm", linestyle='None', markersize=2,\
marker='s', markerfacecolor='c', markeredgecolor='c')
plt.plot(timeseries_EPIC_388[:, 0], timeseries_EPIC_388[:, 1],\
label = "AOD EPIC at 388nm", linestyle='None', markersize=2,\
marker='o', markerfacecolor='c', markeredgecolor='c')
plt.plot(timeseries_EPIC_500[:, 0], timeseries_EPIC_500[:, 1],\
label = "AOD EPIC at 500nm", linestyle='None', markersize=2,\
marker='D', markerfacecolor='c', markeredgecolor='c')
# Set the range of x-axis
l_lim = 0.
u_lim = ((dt_fin - dt_ini).total_seconds() + 1.)/86400.
plt.xlim(l_lim, u_lim)
# Set the range of y-axis
l_lim = 0.
u_lim = 2.
plt.ylim(l_lim, u_lim)
plt.xlabel(r'GMT, day from beginning of '+datestamp_ini, fontsize=12)
plt.ylabel(r'AOD', fontsize=12)
plt.legend(loc='lower left')
plt.title(plot_title+str(', %08.4fN %08.4fW' %(POI_lat, -POI_lon)))
plt.savefig(img_name, format='jpg', dpi=300)
plt.close()8.5 plot UV Aerosol Index
plot_title = 'UVAI_'+datestamp_ini+'_'+datestamp_fin+'\n'+POI_name
img_name = 'UVAI_'+datestamp_ini+'_'+datestamp_fin+'_'+POI_name+'.jpg'
plt.plot(timeseries_EPIC_UVAI[:, 0], timeseries_EPIC_UVAI[:, 1],\
label = "UVAI EPIC", linestyle='None', markersize=2,\
marker='o', markerfacecolor='c', markeredgecolor='c')
plt.plot(timeseries_TEMPO_UVAI[:, 0], timeseries_TEMPO_UVAI[:, 1],\
label = "UVAI TEMPO", linestyle='None', markersize=2,\
marker='o', markerfacecolor='m', markeredgecolor='m')
# Set the range of x-axis
l_lim = 0.
u_lim = ((dt_fin - dt_ini).total_seconds() + 1.)/86400.
plt.xlim(l_lim, u_lim)
# Set the range of y-axis
l_lim = -3.
u_lim = 3.
plt.ylim(l_lim, u_lim)
plt.xlabel(r'GMT, day from beginning of '+datestamp_ini, fontsize=12)
plt.ylabel(r'UV Aerosol Index', fontsize=12)
plt.legend(loc='lower left')
plt.title(plot_title+str(', %08.4fN %08.4fW' %(POI_lat, -POI_lon)))
plt.savefig(img_name, format='jpg', dpi=300)
plt.close()