diff --git a/CHANGELOG.md b/CHANGELOG.md index 33c94b09..9e8f311a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added python package for post-processing ObsFcstAna output. + ### Changed - Specify only ntasks_model for SLURM resource request. diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/Get_ObsFcstAna_stats.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Get_ObsFcstAna_stats.py new file mode 100644 index 00000000..cf7e61a0 --- /dev/null +++ b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Get_ObsFcstAna_stats.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 + +""" +Sample script for post-processing GEOSldas ObsFcstAna output into data assimilation diagnostics. +First, compute and store monthly sums and sums of squares and cross-products of raw ObsFcstAna output. +Data assimilation diagnostics ("stats") such as the mean and std-dev of the observation-minus-forecast +residuals can then be diagnosed quickly from these intermediate "sums" files. +Sample script optionally computes and plots: +- Maps of long-term data assimilation diagnostics (see also Plot_stats_maps.py). +- Monthly time series of spatially averaged data assimilation diagnostics (see also Plot_stats_timeseries.py). + +Usage: + 1. Edit "user_config.py" with experiments information. + 2. Run this script as follows (on Discover): + $ module load python/GEOSpyD + $ ./Get_ObsFcstAna_stats.py + + # Background execution: + $ nohup ./Get_ObsFcstAna_stats.py > out.log & + +Authors: Q. Liu, R. Reichle, A. Fox +Last Modified: May 2025 +""" + +import sys; sys.path.append('../../shared/python/') +import warnings; warnings.filterwarnings("ignore") +import os + +import numpy as np + +from datetime import timedelta +from postproc_ObsFcstAna import postproc_ObsFcstAna +from user_config import get_config # user-defined config info + +# --- +# +# If the script is run in the background, uncomment the following lines to see the redirected +# standard output in the out.log file immediately. When the lines are commented out, the redirected +# standard output will not appear in the out.log file until the job has completed. +# If the script is run in the foreground, the lines must be commented out. +# +#import io +#sys.stdout = io.TextIOWrapper(open(sys.stdout.fileno(), 'wb', 0), write_through=True) +#sys.stderr = io.TextIOWrapper(open(sys.stderr.fileno(), 'wb', 0), write_through=True) +# +# --- + +def main(): + + config = get_config() # edit user-defined inputs in user_config.py + + exp_list = config['exp_list'] + start_time = config['start_time'] + end_time = config['end_time'] + sum_path = config['sum_path'] + out_path = config['out_path'] + + # ------------------------------------------------------------------------------------ + # Postprocess raw ObsFcstAna output data into monthly sums + + # Initialize the postprocessing object + postproc = postproc_ObsFcstAna(exp_list, start_time, end_time, sum_path=sum_path) + + # Compute and save monthly sums + postproc.save_monthly_sums() + + # -------------------------------------------------------------------------------------- + # Optionally compute long-term temporal/spatial statistics and create plots. + # The plotting scripts can also run standalone using the individual Plot_stats_*.py scripts, + # as long as the monthly sum files are available. + + plot_maps = False + plot_timeseries = False + + if plot_maps: # Compute long-term temporal stats and plot maps + + stats_file = out_path + 'temporal_stats_'+exp_list[0]['exptag']+ '_'+start_time.strftime('%Y%m%d')+'_'+ \ + (end_time+timedelta(days=-1)).strftime('%Y%m%d')+'.nc4' + + # temporal_stats is a dictionary that contains all mean/variances fields for computing long-term O-F/O-A stats + # each field has the dimension [N_tile, N_species] + + temporal_stats = postproc.calc_temporal_stats_from_sums(write_to_nc=True, fout_stats=stats_file) + + # Example to plot some O-F maps + from Plot_stats_maps import plot_OmF_maps + plot_OmF_maps(postproc, temporal_stats, fig_path=out_path ) + + if plot_timeseries: # Compute spatially averaged stats and plot monthly time series of stats + + from Plot_stats_timeseries import Plot_monthly_OmF_bars + Plot_monthly_OmF_bars(postproc, fig_path=out_path) + +if __name__ == '__main__': + main() + +# ====================== EOF ========================================================= diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_maps.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_maps.py new file mode 100644 index 00000000..ed44a8db --- /dev/null +++ b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_maps.py @@ -0,0 +1,201 @@ +#!/usr/bin/env python3 + +""" +Sample script for plotting maps of long-term data assimilation diagnostics. +Computes Nobs-weighted avg of each metric across all species. +Requires saved files with monthly sums (see Get_ObsFcstAna_stat.py). +Stats of *normalized* OmFs are approximated! +""" + +import sys; sys.path.append('../../shared/python/') +import warnings; warnings.filterwarnings("ignore") + +import numpy as np +import matplotlib.pyplot as plt + +from datetime import timedelta + +from remap_1d_to_2d import remap_1d_to_2d +from plot import plotMap +from EASEv2 import EASEv2_ind2latlon + +def plot_OmF_maps(postproc_obj, stats, fig_path='./'): + + start_time = postproc_obj.start_time + end_time = postproc_obj.end_time + domain = postproc_obj.domain + tc = postproc_obj.tilecoord + + # Sample of final compuation of selected diagnostic metrics + + Nmin = 20 + + # Then compute metrics of O-F, O-A, etc. based on above computed + N_data = stats['N_data'] + O_mean = stats['obs_mean'] + # mean(x-y) = E[x] - E[y] + OmF_mean = stats['obs_mean'] - stats['fcst_mean'] + OmA_mean = stats['obs_mean'] - stats['ana_mean'] + # var(x-y) = var(x) + var(y) - 2cov(x,y) + # cov(x,y) = E[xy] - E[x]E[y] + OmF_stdv = np.sqrt(stats['obs_variance'] + stats['fcst_variance'] - \ + 2 * (stats['oxf_mean'] - stats['obs_mean']*stats['fcst_mean'])) + + OmA_stdv = np.sqrt(stats['obs_variance'] + stats['ana_variance'] - \ + 2 * (stats['oxa_mean'] - stats['obs_mean']*stats['ana_mean'])) + + # ***************************************************************************************** + # The time series mean and std-dev of the *normalized* OmF computed here are APPROXIMATED! + # ***************************************************************************************** + # Here, we first compute the stats of the OmF time series and then normalize using + # the time-avg "obsvar" and "fcstvar" values. + # Since "fcstvar" changes with time, the OmF values should be normalized at each time + # step (as in the older matlab scripts), and then the time series stats can be computed. + # To compute the exact stats with this python package, the sum and sum-of-squares of + # the normalized OmF values would need to be added into the sums files. + # + OmF_norm_mean = OmF_mean / np.sqrt(stats['obsvar_mean'] + stats['fcstvar_mean']) # APPROXIMATED stat! + OmF_norm_stdv = np.sqrt(OmF_stdv**2 / (stats['obsvar_mean'] + stats['fcstvar_mean']) ) # APPROXIMATED stat! + + # Mask out data points with insufficent observations using the Nmin threshold + # Do NOT apply to N_data + OmF_mean[ N_data < Nmin] = np.nan + OmF_stdv[ N_data < Nmin] = np.nan + OmF_norm_mean[N_data < Nmin] = np.nan + OmF_norm_stdv[N_data < Nmin] = np.nan + OmA_mean[ N_data < Nmin] = np.nan + OmA_stdv[ N_data < Nmin] = np.nan + N_data[ N_data < Nmin] = 0 + + # Compute Nobs-weighted avg of each metric across all species. + # Typically used for SMAP Tb_h/h from asc and desc overpasses, + # or ASCAT soil moisture from Metop-A/B/C. + # DOES NOT MAKE SENSE IF, SAY, SPECIES HAVE DIFFERENT UNITS! + Nobs_data = np.nansum( N_data, axis=1) + OmF_mean = np.nansum(OmF_mean *N_data, axis=1)/Nobs_data + OmF_stdv = np.nansum(OmF_stdv *N_data, axis=1)/Nobs_data + OmF_norm_mean = np.nansum(OmF_norm_mean*N_data, axis=1)/Nobs_data + OmF_norm_stdv = np.nansum(OmF_norm_stdv*N_data, axis=1)/Nobs_data + OmA_mean = np.nansum(OmA_mean *N_data, axis=1)/Nobs_data + OmA_stdv = np.nansum(OmA_stdv *N_data, axis=1)/Nobs_data + + # Plotting + exptag = postproc_obj.exptag + + fig, axes = plt.subplots(2,2, figsize=(18,10)) + plt.rcParams.update({'font.size':14}) + + for i in np.arange(2): + for j in np.arange(2): + units = '[k]' + if i == 0 and j == 0: + tile_data = Nobs_data + # crange is [cmin, cmax] + crange =[0, np.ceil((end_time-start_time).days/150)*300] + colormap = plt.get_cmap('jet',20) + title_txt = exptag + ' Tb Nobs '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m') + units = '[-]' + if i == 0 and j ==1: + tile_data = OmF_mean + crange =[-3, 3] + colormap = plt.get_cmap('bwr', 15) + title_txt = exptag + ' Tb O-F mean '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m') + if i == 1 and j == 0: + tile_data = OmF_stdv + crange =[0, 15] + colormap = plt.get_cmap ('jet',15) + title_txt = exptag + ' Tb O-F stdv '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m') + if i == 1 and j == 1: + tile_data = OmF_norm_stdv + crange =[0, 15] + colormap = plt.get_cmap ('jet',15) + title_txt = exptag + ' Tb normalized O-F stdv (approx!) '+ start_time.strftime('%Y%m%d')+'_'+end_time.strftime('%Y%m%d') + + colormap.set_bad(color='0.9') # light grey, 0-black, 1-white + + # Regrid 1d tile_data to 2d grid_data for map plots + if '_M09_' in domain: # special case + grid_data_M09 = np.zeros((1624, 3856)) + np.nan + grid_data_M09[tc['j_indg'],tc['i_indg']] = tile_data + + # Reshape the data into 4x4 blocks + reshaped = grid_data_M09.reshape(1624//4, 4, 3856//4, 4) + + # Combine each 4x4 M09 block into a M36 grid + #if i==0 and j==0: + # grid_data = np.sum(reshaped,axis=(1, 3)) + #else: + # grid_data = np.nanmean(reshaped,axis=(1, 3)) + + grid_data = grid_data_M09[1::4, 2::4] + + # NOT area weighted + wmean = np.nanmean(grid_data) + wabsmean = np.nanmean(np.abs(grid_data)) + if 'normalized' in title_txt: + wabsmean = np.nanmean(np.abs(grid_data-1.)) + + lat_M36, lon_M36 = EASEv2_ind2latlon(np.arange(406), np.arange(964),'M36') + lon_2d,lat_2d = np.meshgrid(lon_M36,lat_M36) + else: + grid_data, uy, ux = remap_1d_to_2d(tile_data, lat_1d = tc['com_lat'], lon_1d = tc['com_lon']) + lon_2d,lat_2d = np.meshgrid(ux, uy) + + # Aear weighted mean and mean(abs) + wmean = np.nansum( tile_data * tc['area'])/np.nansum(~np.isnan(tile_data)*tc['area']) + wabsmean = np.nansum(np.abs(tile_data) * tc['area'])/np.nansum(~np.isnan(tile_data)*tc['area']) + if 'normalized' in title_txt: + wabsmean = np.nansum(np.abs(tile_data-1.) * tc['area'])/np.nansum(~np.isnan(tile_data)*tc['area']) + + if 'normalized' in title_txt: + title_txt = title_txt + '\n' + "avg=%.3f, avg(abs(nstdv-1))=%.3f" % (wmean, wabsmean)+' '+units + elif 'mean' in title_txt: + title_txt = title_txt + '\n' + "avg=%.3f, avg(abs)=%.3f" % (wmean, wabsmean)+' '+units + else: + title_txt = title_txt + '\n' + "avg=%.2f" % (wmean) +' '+units + + if 'normalized' in title_txt: + grid_data = np.log10(grid_data) + crange = [-0.6, 0.45] + + mm, cs = plotMap(grid_data, ax =axes[i,j], lat=lat_2d, lon=lon_2d, cRange=crange, \ + title=title_txt, cmap=colormap, bounding=[-60, 80, -180,180]) + + plt.tight_layout() + # Save figure to file + fig.savefig(fig_path+'Map_OmF_'+ exptag +'_'+start_time.strftime('%Y%m')+'_'+\ + (end_time+timedelta(days=-1)).strftime('%Y%m')+'.png') + #plt.show() + plt.close(fig) + +# ----------------------------------------------------------------------------------------------- + +if __name__ == '__main__': + + from postproc_ObsFcstAna import postproc_ObsFcstAna + from user_config import get_config + + config = get_config() # edit user-defined inputs in user_config.py + + exp_list = config['exp_list'] + start_time = config['start_time'] + end_time = config['end_time'] + sum_path = config['sum_path'] + out_path = config['out_path'] + + # ------------------------------------------------------------------------------------ + # Initialize the postprocessing object + postproc = postproc_ObsFcstAna(exp_list, start_time, end_time, sum_path=sum_path) + + # Compute long-term temporal stats and plot maps + stats_file = out_path + 'temporal_stats_'+exp_list[0]['exptag']+'_'+ start_time.strftime('%Y%m%d')+'_'+ \ + (end_time+timedelta(days=-1)).strftime('%Y%m%d')+'.nc4' + + # temporal_stats is a dictionary that contains all mean/variances fields for computing long-term O-F/O-A stats + # each field has the dimension [N_tile, N_species] + + temporal_stats = postproc.calc_temporal_stats_from_sums(write_to_nc=True, fout_stats=stats_file) + + plot_OmF_maps(postproc, temporal_stats, fig_path=out_path ) + +# ====================== EOF ========================================================= diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_timeseries.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_timeseries.py new file mode 100644 index 00000000..db72e9f0 --- /dev/null +++ b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_timeseries.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 + +""" +Sample script for plotting monthly time series of spatially averaged data assimilation diagnostics. +Computes Nobs-weighted avg of each metric across all species. +Requires saved files with monthly sums (see Get_ObsFcstAna_stat.py). +""" + +import sys; sys.path.append('../../shared/python/') +import warnings; warnings.filterwarnings("ignore") +import os + +import numpy as np +import matplotlib.pyplot as plt + +from datetime import datetime, timedelta +from dateutil.relativedelta import relativedelta + +def Plot_monthly_OmF_bars(postproc_obj, fig_path='./'): + import pickle + exptag = postproc_obj.exptag + start_month = postproc_obj.start_time + end_month = postproc_obj.end_time + + stats_file = fig_path + 'spatial_stats_'+exptag+'_'+start_month.strftime('%Y%m')+ \ + '_'+(end_month+timedelta(days=-1)).strftime('%Y%m')+'.pkl' + + if os.path.isfile(stats_file): + + with open(stats_file,'rb') as file: + stats_dict = pickle.load(file) + date_vec = stats_dict['date_vec'] + + else: + + # Initialize monthly metrics as list + OmF_mean =[] + OmF_stdv =[] + OmA_mean =[] + OmA_stdv =[] + Ndata =[] + date_vec =[] + + # Time loop + current_month = start_month + while current_month < end_month: + print('compute monthly spatial mean stats for '+ current_month.strftime('%Y%m')) + + OmFm,OmFs,OmAm,OmAs,Nobsm = postproc_obj.calc_spatial_stats_from_sums(current_month) + + # Compute Nobs-weighted avg of each metric across all species. + # Typically used for SMAP Tb_h/h from asc and desc overpasses, + # or ASCAT soil moisture from Metop-A/B/C. + # DOES NOT MAKE SENSE IF, SAY, SPECIES HAVE DIFFERENT UNITS! + Nobsm = np.nansum( Nobsm) + OmFm = np.nansum(OmFm*Nobsm)/Nobsm + OmFs = np.nansum(OmFs*Nobsm)/Nobsm + OmAm = np.nansum(OmAm*Nobsm)/Nobsm + OmAs = np.nansum(OmAs*Nobsm)/Nobsm + + Ndata.append(Nobsm) + + OmF_mean.append(OmFm) + OmF_stdv.append(OmFs) + OmA_mean.append(OmAm) + OmA_stdv.append(OmAs) + + date_vec.append(current_month.strftime('%Y%m')) + current_month = current_month + relativedelta(months=1) + + # Store stats in a dictionary for easier saving and referencing + stats_dict = {"OmF_mean":OmF_mean, "OmF_stdv":OmF_stdv, + "OmA_mean":OmA_mean, "OmA_stdv":OmA_stdv, + "Ndata": Ndata, "date_vec":date_vec} + + if stats_file is not None: + with open(stats_file,'wb') as file: + pickle.dump(stats_dict,file) + + plot_var = 'OmF_mean' + + fig, ax = plt.subplots(figsize=(10,4)) + bars = ax.bar(date_vec, stats_dict[plot_var]) + + plt.grid(True, linestyle='--', alpha=0.5) + + plt.xticks(ticks=date_vec[::2], labels=date_vec[::2]) + plt.title(exptag+ 'monthly '+plot_var) + plt.xlim(-1, len(date_vec)+1) + plt.ylim(-.1, 2.) + + plt.tight_layout() + #plt.show() + plt.savefig(fig_path+'Bars_'+plot_var+exptag+date_vec[0]+'_'+\ + date_vec[-1]+'.png') + plt.close() + +if __name__ == "__main__": + + from postproc_ObsFcstAna import postproc_ObsFcstAna + from user_config import get_config + + config = get_config() # edit user-defined inputs in user_config.py + + exp_list = config['exp_list'] + start_time = config['start_time'] + end_time = config['end_time'] + sum_path = config['sum_path'] + out_path = config['out_path'] + + # ------------------------------------------------------------------------------------ + # + # Initialize the postprocessing object + postproc = postproc_ObsFcstAna(exp_list, start_time, end_time, sum_path=sum_path) + + # Compute spatial stats and plot monthly O-F stats bars + Plot_monthly_OmF_bars(postproc, fig_path=out_path) + + +# ====================== EOF ========================================================= diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/helper/write_nc4.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/helper/write_nc4.py new file mode 100644 index 00000000..8090cc89 --- /dev/null +++ b/GEOSldas_App/util/postproc/ObsFcstAna_stats/helper/write_nc4.py @@ -0,0 +1,62 @@ + +# collection of write routines for ObsFcstAna sums and stats + +from netCDF4 import Dataset, date2num + +import numpy as np + +# -------------------------------------------------------------------------------- + +def write_sums_nc4(file_path, N_data, data_sum, data2_sum, oxf_sum, oxa_sum, fxa_sum, obs_param): + + nc = Dataset(file_path,'w',formate='NETCDF4') + + tile = nc.createDimension('tile', N_data.shape[0]) + species = nc.createDimension('species', N_data.shape[1]) + + data = nc.createVariable('obs_param_assim', 'c', ('species', ), zlib=True, complevel=4) + for i in range(len(obs_param)): + data[i] = obs_param[i]['assim'] + + data = nc.createVariable('N_data', 'i4', ('tile','species', ), zlib=True, complevel=4) + data[:,:] = N_data + + data = nc.createVariable('obsxfcst_sum', 'f4', ('tile','species', ), zlib=True, complevel=4) + data[:,:] = oxf_sum + + data = nc.createVariable('obsxana_sum', 'f4', ('tile','species', ), zlib=True, complevel=4) + data[:,:] = oxa_sum + + data = nc.createVariable('fcstxana_sum', 'f4', ('tile','species', ), zlib=True, complevel=4) + data[:,:] = fxa_sum + + for key, value in data_sum.items(): + varname = key+'_sum' + data = nc.createVariable(varname, 'f4', ('tile','species', ), zlib=True, complevel=4) + data[:,:] = value + + for key, value in data2_sum.items(): + varname = key+'2_sum' + data = nc.createVariable(varname, 'f4', ('tile','species', ), zlib=True, complevel=4) + data[:,:] = value + + nc.close() + +# -------------------------------------------------------------------------------- + +def write_stats_nc4(file_path, stats): + + nc = Dataset(file_path,'w',formate='NETCDF4') + + tile = nc.createDimension('tile', stats['N_data'].shape[0]) + species = nc.createDimension('species', stats['N_data'].shape[1]) + + for key, value in stats.items(): + data = nc.createVariable(key,'f4',('tile','species', ), \ + fill_value=-9999.0, zlib=True, complevel=4) + data[:,:] = value + + nc.close() + + +# ======================= EOF ========================================================= diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py new file mode 100644 index 00000000..41b86df9 --- /dev/null +++ b/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py @@ -0,0 +1,423 @@ + +# Tool to preprocess GEOSldas ObsFcstAna output into monthly sums, sums of squares, and sums of cross-products +# +# qliu, amfox, rreichle - May 2025 + +import numpy as np +import os +import yaml + +import sys; sys.path.append('../../shared/python/') + +import warnings; warnings.filterwarnings("ignore") + +from datetime import datetime, timedelta +from dateutil.relativedelta import relativedelta +from netCDF4 import Dataset, date2num +from read_GEOSldas import read_ObsFcstAna, read_tilecoord, read_obs_param + +from helper.write_nc4 import write_sums_nc4, write_stats_nc4 + +class postproc_ObsFcstAna: + + def __init__(self, exp_list, start_time, end_time, sum_path='./'): + self.exp_list = exp_list + self.expdir_list = [item['expdir'] for item in exp_list] + self.expid_list = [item['expid'] for item in exp_list] + self.exptag = exp_list[0]['exptag'] + self.domain = exp_list[0]['domain'] + self.start_time = start_time + self.end_time = end_time + self.da_t0 = exp_list[0]['da_t0'] + self.da_dt = exp_list[0]['da_dt'] + self.var_list = ['obs_obs','obs_obsvar','obs_fcst','obs_fcstvar','obs_ana','obs_anavar'] + self.tilecoord = exp_list[0]['tilecoord'] + self.obsparam_list = [item['obsparam'] for item in exp_list] + self.sum_path = sum_path + + # Determine experiment that supplies obs data + self.obs_from = -1 + for exp_idx, exp in enumerate(exp_list): + if exp.get('use_obs', None): # found use_obs=True + if self.obs_from >= 0: + print("ERROR: use_obs=True in multiple experiments. Edit user_config.py to remove conflict." ) + sys.exit() + else: + self.obs_from = exp_idx + if self.obs_from < 0: self.obs_from = 0 # by default, obs data are from exp_list[0] + print(f"obs data are from {exp_list[self.obs_from]['expid']}") + + # Verify the configuration every time when current class is initialized + # to avoid saving sums with different configs in the same directory + + # Same configurations should have identical values for these fields + config_verify = ['expdir','expid','exptag','domain','use_obs','species_list'] + + # Construct config for each experiment + config_list = [] + for exp in exp_list: + config_list.append({var:exp[var] for var in config_verify if var in exp}) + + # File of configuration for verification + f_config = self.sum_path + '/' + self.exptag + '_config.yaml' + + # Save a new file or compare current configuration with previously saved + if not os.path.exists(f_config): + with open(f_config, 'w') as f: + yaml.dump(config_list, f, default_flow_style=False) + print(f'Configuration saved to {f_config}') + else: + with open(f_config,'r') as f: + saved_exp_list = yaml.safe_load(f) + print(f'Found saved configuration') + + if config_list != saved_exp_list: + print('user configuration is different from previously saved '+f_config) + sys.exit() + # ---------------------------------------------------------------------------------------------------------- + # + # Function to compute monthly sums from x-hourly ObsFcstAna data for a single month. + + def compute_monthly_sums(self, date_time): + + expdir_list = self.expdir_list + expid_list = self.expid_list + domain = self.domain + tc = self.tilecoord + obsparam_list = self.obsparam_list + var_list = self.var_list + obs_from = self.obs_from + da_dt = self.da_dt + + n_tile = tc['N_tile'] + n_spec = len(obsparam_list[0]) + + date_time = date_time.replace(hour=self.da_t0) + end_time = date_time + relativedelta(months=1) + + data_sum = {} + data2_sum = {} + + N_data = np.zeros((n_tile, n_spec)) + oxf_sum = np.zeros((n_tile, n_spec)) + oxa_sum = np.zeros((n_tile, n_spec)) + fxa_sum = np.zeros((n_tile, n_spec)) + + for var in var_list: + data_sum[ var] = np.zeros((n_tile, n_spec)) + data2_sum[var] = np.zeros((n_tile, n_spec)) + + while date_time < end_time: + + # read the list of experiments at each time step (OFA="ObsFcstAna") + OFA_list = [] + for i in range(len(expdir_list)): + fname = expdir_list[i]+expid_list[i]+'/output/'+domain+'/ana/ens_avg/Y'+ \ + date_time.strftime('%Y') + '/M' + \ + date_time.strftime('%m') + '/' + \ + expid_list[i]+'.ens_avg.ldas_ObsFcstAna.' + \ + date_time.strftime('%Y%m%d_%H%M') +'z.bin' + if os.path.isfile(fname): + print('read '+fname) + OFA_list.append(read_ObsFcstAna(fname)) + + data_all=[] + for OFA, obs_param in zip(OFA_list,obsparam_list): + + # Initialize full size variable for an experiment + data_tile={} + for var in var_list: + data_tile[var] = np.zeros((n_tile, n_spec)) + np.nan + + if len(OFA['obs_tilenum']) > 0: + for ispec in np.arange(n_spec): + # check species overall "assim" flag for masking + this_species = int(obs_param[ispec]['species']) + masked_data = {} + + # get mask + if obs_param[ispec]['assim'] == 'T': + is_valid = np.logical_and(OFA['obs_species'] == this_species, OFA['obs_assim']==1) + else: + is_valid = OFA['obs_species'] == this_species + + tile_idx = OFA['obs_tilenum'][is_valid]-1 + + for var in var_list: + masked_data[var] = OFA[var][is_valid] + + for var in var_list: + data_tile[var][tile_idx, ispec] = masked_data[var] + + data_all.append(data_tile) + + # cross-mask over all experiments + is_cross_valid = ~np.isnan(data_all[0]['obs_obs']) + for data in data_all[1:]: + mask = ~np.isnan(data['obs_obs']) + is_cross_valid = np.logical_and(is_cross_valid,mask) + + # reconstruct the output variable dictionary based on input options; + # obs_obs and obs_obsvar are from exp_list[obs_from], the rest are from exp_list[0] + data_tile = {} + for var in var_list: + if 'obs_obs' in var: + data_tile[var] = data_all[obs_from][var] + else: + data_tile[var] = data_all[0][var] + + is_valid = is_cross_valid + + N_data[ is_valid] += 1 + oxf_sum[is_valid] += data_tile['obs_obs' ][is_valid] * data_tile['obs_fcst'][is_valid] + oxa_sum[is_valid] += data_tile['obs_obs' ][is_valid] * data_tile['obs_ana' ][is_valid] + fxa_sum[is_valid] += data_tile['obs_fcst'][is_valid] * data_tile['obs_ana' ][is_valid] + + for var in var_list: + data_sum[ var][is_valid] += data_tile[var][is_valid] + data2_sum[var][is_valid] += data_tile[var][is_valid] **2 + + date_time = date_time + timedelta(seconds=da_dt) + + return N_data, data_sum, data2_sum, oxf_sum, oxa_sum, fxa_sum + + # ---------------------------------------------------------------------------------------------------------- + # + # Function to compute monthly sums and save results in nc4 files for all months in [start_time, end_time]. + # + # Skips computation/saving of monthly sums output if file already exists. + + def save_monthly_sums(self): + expdir_list = self.expdir_list + expid_list = self.expid_list + domain = self.domain + var_list = self.var_list + tc = self.tilecoord + obsparam_list = self.obsparam_list + + date_time = self.start_time + + while date_time < self.end_time: # loop through months + + # make monthly file output directory + mo_path = self.sum_path + '/Y'+ date_time.strftime('%Y') + '/M' + date_time.strftime('%m') + '/' + os.makedirs(mo_path, exist_ok=True) + + fout = self.exptag + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4' + + fout = mo_path + fout + + # skip if output file already exists + if not os.path.isfile(fout): + print('computing monthly sums...') + # compute monthly sums + mN_data, mdata_sum, mdata2_sum, moxf_sum, moxa_sum, mfxa_sum = \ + self.compute_monthly_sums(date_time) + + # save monthly sums in nc4 file + write_sums_nc4(fout, mN_data,mdata_sum, mdata2_sum, moxf_sum, moxa_sum, mfxa_sum, obsparam_list[0]) + else: + print('file exists, skipping '+fout) + + date_time = date_time + relativedelta(months=1) + + # ---------------------------------------------------------------------------------------------------------- + # + # Function to compute long-term temporal statistics of individual species based on monthly sums. + # + # Assumes that monthly sums files have been saved [see save_monthly_sums()]. + + def calc_temporal_stats_from_sums(self, write_to_nc=True, fout_stats='./stats.nc4'): + + if os.path.isfile(fout_stats): + + print('reading stats nc4 file '+ fout_stats) + stats = {} + with Dataset(fout_stats,'r') as nc: + for key, value in nc.variables.items(): + stats[key] = value[:].filled(np.nan) + + else: + # Variable list for computing sum and sum-of-squares + var_list = self.var_list + + # Read tilecoord and obsparam for tile and obs species information + n_tile = self.tilecoord['N_tile'] + n_spec = len(self.obsparam_list[0]) + + # Initialize statistical metrics + data_sum = {} + data2_sum = {} + N_data = np.zeros((n_tile, n_spec)) + oxf_sum = np.zeros((n_tile, n_spec)) + oxa_sum = np.zeros((n_tile, n_spec)) + fxa_sum = np.zeros((n_tile, n_spec)) + + for var in var_list: + data_sum[ var] = np.zeros((n_tile, n_spec)) + data2_sum[var] = np.zeros((n_tile, n_spec)) + + # Time loop: processing data at monthly time step + + date_time = self.start_time + + while date_time < self.end_time: # loop through months + + fpath = self.sum_path + '/Y'+ date_time.strftime('%Y') + '/M' + date_time.strftime('%m') + '/' + + fout = self.exptag + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4' + + fout = fpath + fout + + # Read monthly data + if os.path.isfile(fout): + print('read sums from monthly file: '+fout) + mdata_sum = {} + mdata2_sum = {} + with Dataset(fout,'r') as nc: + mN_data = nc.variables['N_data' ][:] + moxf_sum = nc.variables['obsxfcst_sum'][:] + moxa_sum = nc.variables['obsxana_sum' ][:] + mfxa_sum = nc.variables['fcstxana_sum'][:] + for var in var_list: + mdata_sum[ var] = nc.variables[var+'_sum' ][:] + mdata2_sum[var] = nc.variables[var+'2_sum'][:] + + # Aggregate monthly data + N_data += mN_data + oxf_sum += moxf_sum + oxa_sum += moxa_sum + fxa_sum += mfxa_sum + + for var in var_list: + data_sum[ var] += mdata_sum[ var] + data2_sum[var] += mdata2_sum[var] + else: + raise FileNotFoundError(f"File {fout} does not exist, run save_monthly_sums() first") + + date_time = date_time + relativedelta(months=1) + + # Compute stats (DA diagnostics) from accumulated sums data. + data_mean = {} + data2_mean = {} + data_var = {} + + # calculation + for var in var_list: + data_sum[var][ N_data == 0] = np.nan + data2_sum[var][N_data == 0] = np.nan + + data_mean[ var] = data_sum[ var] / N_data + data2_mean[var] = data2_sum[var] / N_data + # var(x) = E[x2] - (E[x])^2 + data_var[var] = data2_mean[var] - data_mean[var]**2 + + oxf_sum[N_data == 0] = np.nan + oxa_sum[N_data == 0] = np.nan + fxa_sum[N_data == 0] = np.nan + # E[xy] + oxf_mean = oxf_sum / N_data + oxa_mean = oxa_sum / N_data + fxa_mean = fxa_sum / N_data + + stats = {} + for var in var_list: + stats[var[4:]+'_mean'] = data_mean[var] + stats[var[4:]+'_variance'] = data_var[ var] + stats['oxf_mean'] = oxf_mean + stats['oxa_mean'] = oxa_mean + stats['fxa_mean'] = fxa_mean + stats['N_data'] = N_data + + if write_to_nc: + print('writing stats nc4 file: '+fout_stats) + write_stats_nc4(fout_stats, stats) + + return stats + + # ---------------------------------------------------------------------------------------------------------- + # + # Function to compute the O-F/O-A *spatial* statistics for a *single* month based on + # previously saved monthly sums. + # Individual temporal and grid cell DA diagnostic values within a month are aggregated first; + # the monthly Ndata/mean/stdv are derived from the aggregated sample. Consequently, in the + # final DA diagnostics, each obs value gets equal weight. Note that this differs from computing + # the straight spatial avg across a map of a given DA diagnostic. The latter approach gives + # the same weight to each location, regardless of how many obs are available at the location. + + def calc_spatial_stats_from_sums(self, date_time): + + var_list = ['obs_obs', 'obs_fcst','obs_ana'] + + mo_path = self.sum_path + '/Y'+ date_time.strftime('%Y') + '/M' + date_time.strftime('%m') + '/' + fnc4_sums = mo_path + self.exptag + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4' + + mdata_sum = {} + mdata2_sum = {} + + try: + with Dataset(fnc4_sums,'r') as nc: + mN_data = nc.variables['N_data' ][:] + moxf_sum = nc.variables['obsxfcst_sum'][:] + moxa_sum = nc.variables['obsxana_sum' ][:] + moxf_sum[mN_data == 0] = np.nan + moxa_sum[mN_data == 0] = np.nan + for var in var_list: + mdata_sum[ var] = nc.variables[var+'_sum' ][:] + mdata2_sum[var] = nc.variables[var+'2_sum'][:] + mdata_sum[ var][mN_data == 0] = np.nan + mdata2_sum[var][mN_data == 0] = np.nan + except FileNotFoundError: + print(f"Error: File '{filename}' not found.") + sys.exit(1) + except Exception as e: + print(f"An unexpected error occurred: {e}") + sys.exit(1) + + # Make sure only aggregate tiles with valid values for all variables + for var in var_list: + mN_data = mN_data.astype(float) + mN_data[np.isnan(mdata_sum[var])] = np.nan + mN_data[mN_data == 0] = np.nan + + # cross mask before aggregating tile values + for var in var_list: + mdata_sum[ var][np.isnan(mN_data)] = np.nan + mdata2_sum[var][np.isnan(mN_data)] = np.nan + moxf_sum[np.isnan(mN_data)] = np.nan + moxa_sum[np.isnan(mN_data)] = np.nan + + # Aggregate data of all tiles + N_data = np.nansum(mN_data, axis=0) + OxF_mean = np.nansum(moxf_sum,axis=0)/N_data + OxA_mean = np.nansum(moxa_sum,axis=0)/N_data + data_mean = {} + data2_mean = {} + data_var = {} + for var in var_list: + data_mean[ var] = np.nansum(mdata_sum[var ],axis=0)/N_data + data2_mean[var] = np.nansum(mdata2_sum[var],axis=0)/N_data + # var(x) = E[x2] - (E[x])^2 + data_var[var] = data2_mean[var] - data_mean[var]**2 + + # Compute metrics of O-F, O-A, etc. based on above stats + O_mean = data_mean['obs_obs'] + F_mean = data_mean['obs_fcst'] + A_mean = data_mean['obs_ana'] + + O_var = data_var[ 'obs_obs'] + F_var = data_var[ 'obs_fcst'] + A_var = data_var[ 'obs_ana'] + + # mean(x-y) = E[x] - E[y] + OmF_mean = O_mean - F_mean + OmA_mean = O_mean - A_mean + + # var(x-y) = var(x) + var(y) - 2cov(x,y) + # cov(x,y) = E[xy] - E[x]E[y] + OmF_stdv = np.sqrt(O_var + F_var - 2 * (OxF_mean - O_mean*F_mean)) + OmA_stdv = np.sqrt(O_var + A_var - 2 * (OxA_mean - O_mean*A_mean)) + + return OmF_mean, OmF_stdv, OmA_mean, OmA_stdv, N_data + +# ============== EOF ==================================================================== diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/user_config.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/user_config.py new file mode 100644 index 00000000..27e484fc --- /dev/null +++ b/GEOSldas_App/util/postproc/ObsFcstAna_stats/user_config.py @@ -0,0 +1,146 @@ + +import sys; sys.path.append('../../shared/python/') +import warnings; warnings.filterwarnings("ignore") +import os + +import numpy as np + +from datetime import datetime, timedelta +from read_GEOSldas import read_tilecoord, read_obs_param +from postproc_ObsFcstAna import postproc_ObsFcstAna + +def get_config(): + + # ========================================================================================= + # + # User-defined inputs for post-processing of ObsFcstAna output + + # Range of months to process: + + start_year = 2015 + start_month = 4 + last_year = 2016 + last_month = 4 + + # Sums or stats will be processed for exp_main: + + exp_main = { + 'expdir' : '/discover/nobackup/projects/gmao/merra/iau/merra_land/SMAP_runs/SMAP_Nature_v11/', + 'expid' : 'DAv8_SMOSSMAP', # GEOSldas exp ID of simulation + 'exptag' : 'DAMulti_SMAP', # string used in output file names for sums & stats, + # can be same as expid or different -- e.g., reflect info + # about subset of included species or about cross-masking + 'domain' : 'SMAP_EASEv2_M36_GLOBAL', + 'da_t0' : 3, # (fractional) UTC hour of first land analysis + 'da_dt' : 10800, # ObsFcstAna file interval in seconds + 'species_list' : [5,6,7,8], # indices of species to be processed + 'obsparam_time' : "20150401_0000" # time stamp of obsparam file (YYYYMMDD_HHMM) + } + + # Optional experiment(s) can be added for cross-masking or extracting obs from a different experiment. + # + # All optional experiments and the main experiment must have identical tilecoords. + # If the default "species" number/order do not match, set "species_list" accordingly to force a match. + # Output will be cross-masked between all specified experiments. + # + # Forecasts and analyses are always from the main experiment. + # Observations are from the experiment with 'use_obs' set to True (default is exp_main). The most + # likely use case for reading obs from a supplemental experiment is when computing OmF etc diagnostics + # for an open loop experiment that only has unscaled obs, and _scaled_ obs must be read from a + # coresponding DA experiment. + + exp_sup1 = { + 'expdir' : '/discover/nobackup/projects/gmao/merra/iau/merra_land/SMAP_runs/SMAP_Nature_v11/', + 'expid' : 'DAv8_M36', + 'domain' : 'SMAP_EASEv2_M36_GLOBAL', + 'use_obs' : True, # if True, use obs data from this exp + 'species_list' : [1,2,3,4], # indices of species to be processed, + # must identify same species as selected in main exp + 'obsparam_time' : "20150401_0000" # time stamp of obsparam file (YYYYMMDD_HHMM) + } + + # Convert experiments input to a list; first entry must be exp_main: + + #exp_list = [exp_main] # no cross-masking + exp_list = [exp_main, exp_sup1] # cross-mask exp_main with exp_sup1 + + # Top level directory for all output from this package: + + out_path = '/discover/nobackup/[USERNAME]/SMAP_test/' + + # Directory for monthly sum files: + # - Can use the experiment directory or a different path. + # - Automatically appends /Yyyyy/Mmm/ for individual months. + + sum_path = out_path+'/'+exp_main['expid']+'/output/'+exp_main['domain']+'/ana/ens_avg/' + + # + # ===================== end of user-defined inputs ================================================= + + # process time range info; end_time is first of month after (end_year, end_month) + + if last_month==12 : + end_month = 1 + end_year = last_year + 1 + else : + end_month = last_month + 1 + end_year = last_year + + start_time = datetime( start_year, start_month, 1) + end_time = datetime( end_year, end_month, 1) + + # Get tilecoord and obsparam information for each experiment + + for exp in exp_list: + expdir = exp['expdir'] + expid = exp['expid'] + domain = exp['domain'] + + YYYY = exp['obsparam_time'][0:4] + MM = exp['obsparam_time'][4:6] + + fop = expdir+expid+'/output/'+domain+'/rc_out/Y'+YYYY+'/M'+MM+'/'+expid+'.ldas_obsparam.'+exp['obsparam_time']+'z.txt' + obsparam_orig = read_obs_param(fop) + + # get the species list, default is all species + species_list = exp.get('species_list', [ int(obsparam_orig[i]['species']) for i in range(len(obsparam_orig)) ]) + + # subset obsparam_orig based on species_list; keep order of obsparam_orig (independent of order of species_list) + obsparam = [] + for i in range(len(obsparam_orig)): + if int(obsparam_orig[i]['species']) in species_list: + obsparam.append(obsparam_orig[i]) + + ftc = expdir+expid+'/output/'+ domain+'/rc_out/'+ expid+'.ldas_tilecoord.bin' + tc = read_tilecoord(ftc) + + # add tilecoord and obsparam into to exp + exp.update({'tilecoord':tc, 'obsparam':obsparam}) + + # verify that obs species match across experiments + for exp_idx, exp in enumerate(exp_list) : + obsparam = exp.get('obsparam') + if exp_idx==0 : + obsparam_main = obsparam + else: + if len(obsparam) != len(obsparam_main) : + print("ERROR: 'obsparam' mismatch (length)! Check 'select_species' input in user_config.py." ) + sys.exit() + else: + for a, b in zip(obsparam, obsparam_main) : + if a['descr'] != b['descr'] : + print("ERROR: 'obsparam' mismatch (descr)! Check 'select_species' input in user_config.py." ) + sys.exit() + + # wrap up config + config ={ + 'exp_list' : exp_list, + 'start_time' : start_time, + 'end_time' : end_time, + 'sum_path' : sum_path, + 'out_path' : out_path, + } + + return config + +# ====================== EOF ========================================================= diff --git a/GEOSldas_App/util/shared/matlab/EASEv2_ind2latlon.m b/GEOSldas_App/util/shared/matlab/EASEv2_ind2latlon.m index c2d61140..dacbb4ec 100644 --- a/GEOSldas_App/util/shared/matlab/EASEv2_ind2latlon.m +++ b/GEOSldas_App/util/shared/matlab/EASEv2_ind2latlon.m @@ -1,202 +1,202 @@ -% -% SMAPEASE2INVERSE The principal function is to perform inverse transformation -% from (row,col)'s to (lat,lon)'s for a set of nested EASE -% grids defined at 1, 3, 9, and 36km grid resolutions. These -% grids are all based on the EASE-Grid 2.0 specification (WGS84 -% ellipsoid). -% -% SYNTAX [lat,lon] = smapease2inverse(row,col,gridid) -% -% where gridid is a 3-character string enclosed in single -% quotes, in the form of {M|N|S}{01,03,09,36}. This subroutine -% accepts vector inputs and produce vector outputs. -% -% HISTORY This subroutine was adapted from the offical EASE-Grid-2.0 -% conversion utilities (written in IDL) developed by the -% NSIDC. -% -% Note that in NSIDC's original implementation, (row,col) are -% zero-based. In other words, the first cell is (0,0) and the -% last cell is (N-1,M-1), where N and M are the row and column -% dimensions of the array. In this MATLAB implementation, the -% same convention is used. In other words, the end point of -% the first cell is located at (r,c) = (-0.5,-0.5) whereas the -% end point of the last cell is located at (r,c) = (14615.5, -% 34703.5). Thus, -% -% [lat,lon] = smapease2inverse(-0.5,-0.5,'M01') returns: -% lat = 85.044566407398861 -% lon = 1.799999999999994e+02 -% -% [lat,lon] = smapease2inverse(14615.5,34703.5,'M01') returns: -% lat = -85.044566407398861 -% lon = -1.799999999999994e+02 -% -% The polar grids, on the other hand, are more complete in -% terms of latitude coverage: -% -% [lat,lon] = smapease2inverse(8999,8999,'N01') -% lat = 89.993669248945238 -% lon = -135 -% [lat,lon] = smapease2inverse(9000,9000,'N01') -% lat = 89.993669248945238 -% lon = 45 -% -% [lat,lon] = smapease2inverse(8999,8999,'S01') -% lat = -89.993669248945238 -% lon = -45 -% [lat,lon] = smapease2inverse(9000,9000,'S01') -% lat = -89.993669248945238 -% lon = 135 -% -% UPDATE North/south polar projections were added. (03/2012) -% -% REFERENCE Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. -% Savoie (2012): EASE-Grid 2.0: Incremental but Significant -% Improvements for Earth-Gridded Data Sets. ISPRS International -% Journal of Geo-Information, vol. 1, no. 1, pp. 32-45, -% http://www.mdpi.com/2220-9964/1/1/32/ -% -% Steven Chan, 11/2011 -% Email: steven.k.chan@jpl.nasa.gov - -function [lat,lon] = EASEv2_ind2latlon(row,col,gridid) - -% Constants returned by EASE2_GRID_INFO.PRO -projection = gridid(1); -switch gridid - case 'M36' - map_scale_m = 36032.220840584; - cols = 964; - rows = 406; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'M09' - map_scale_m = 9008.055210146; - cols = 3856; - rows = 1624; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'M03' - map_scale_m = 3002.6850700487; - cols = 11568; - rows = 4872; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'M01' - map_scale_m = 1000.89502334956; - cols = 34704; - rows = 14616; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N36' - map_scale_m = 36000.0; - cols = 500; - rows = 500; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N09' - map_scale_m = 9000.0; - cols = 2000; - rows = 2000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N03' - map_scale_m = 3000.0; - cols = 6000; - rows = 6000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N01' - map_scale_m = 1000.0; - cols = 18000; - rows = 18000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S36' - map_scale_m = 36000.0; - cols = 500; - rows = 500; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S09' - map_scale_m = 9000.0; - cols = 2000; - rows = 2000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S03' - map_scale_m = 3000.0; - cols = 6000; - rows = 6000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S01' - map_scale_m = 1000.0; - cols = 18000; - rows = 18000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - otherwise - disp(['ERROR: Incompatible grid specification.']); -end - -% Constants returned by EASE2_MAP_INFO.PRO -map_equatorial_radius_m = 6378137.0; -map_eccentricity = 0.081819190843; -e2 = map_eccentricity^2; -switch projection - case 'M' - map_reference_latitude = 0.0; - map_reference_longitude = 0.0; - map_second_reference_latitude = 30.0; - sin_phi1 = sin(map_second_reference_latitude*pi/180); - cos_phi1 = cos(map_second_reference_latitude*pi/180); - kz = cos_phi1/sqrt(1.0-e2*sin_phi1*sin_phi1); - case 'N' - map_reference_latitude = 90.0; - map_reference_longitude = 0.0; - case 'S' - map_reference_latitude = -90.0; - map_reference_longitude = 0.0; -end - -% Selected calculations inside WGS84_INVERSE.PRO -x = (col-r0)*map_scale_m; -y = (s0-row)*map_scale_m; - -% Selected calculations inside WGS84_INVERSE_XY.PRO -e4 = map_eccentricity^4; -e6 = map_eccentricity^6; -qp = (1.0-e2)*( (1.0/(1.0-e2))-(1.0/(2.0*map_eccentricity))*log((1.0-map_eccentricity)/(1.0+map_eccentricity))); -switch projection - case 'M' - beta = asin(2.0*y*kz/(map_equatorial_radius_m*qp)); - lam = x/(map_equatorial_radius_m*kz); - case 'N' - rho = sqrt(x.^2+y.^2); - beta = asin(1.0-(rho.^2)/(qp*map_equatorial_radius_m^2)); - lam = atan2(x,-y); - case 'S' - rho = sqrt(x.^2+y.^2); - beta = -asin(1.0-(rho.^2)/(qp*map_equatorial_radius_m^2)); - lam = atan2(x,y); -end -phi = beta+(((e2/3.0)+((31.0/180.0)*e4)+((517.0/5040.0)*e6))*sin(2.0*beta))+...\ - ((((23.0/360.0)*e4)+((251.0/3780.0)*e6))*sin(4.0*beta))+(((761.0/45360.0)*e6)*sin(6.0*beta)); -lat = phi*180.0/pi; -lon = map_reference_longitude+(lam*180.0/pi); -msk1 = lon < -180.0; lon(msk1) = lon(msk1) + 360.0; -msk2 = lon > +180.0; lon(msk2) = lon(msk2) - 360.0; -switch projection - case 'N' - idx = lat < 0.0; - lat(idx) = NaN; - lon(idx) = NaN; - case 'S' - idx = lat > 0.0; - lat(idx) = NaN; - lon(idx) = NaN; -end - -% ========= EOF ========================================================= +% +% SMAPEASE2INVERSE The principal function is to perform inverse transformation +% from (row,col)'s to (lat,lon)'s for a set of nested EASE +% grids defined at 1, 3, 9, and 36km grid resolutions. These +% grids are all based on the EASE-Grid 2.0 specification (WGS84 +% ellipsoid). +% +% SYNTAX [lat,lon] = smapease2inverse(row,col,gridid) +% +% where gridid is a 3-character string enclosed in single +% quotes, in the form of {M|N|S}{01,03,09,36}. This subroutine +% accepts vector inputs and produce vector outputs. +% +% HISTORY This subroutine was adapted from the offical EASE-Grid-2.0 +% conversion utilities (written in IDL) developed by the +% NSIDC. +% +% Note that in NSIDC's original implementation, (row,col) are +% zero-based. In other words, the first cell is (0,0) and the +% last cell is (N-1,M-1), where N and M are the row and column +% dimensions of the array. In this MATLAB implementation, the +% same convention is used. In other words, the end point of +% the first cell is located at (r,c) = (-0.5,-0.5) whereas the +% end point of the last cell is located at (r,c) = (14615.5, +% 34703.5). Thus, +% +% [lat,lon] = smapease2inverse(-0.5,-0.5,'M01') returns: +% lat = 85.044566407398861 +% lon = 1.799999999999994e+02 +% +% [lat,lon] = smapease2inverse(14615.5,34703.5,'M01') returns: +% lat = -85.044566407398861 +% lon = -1.799999999999994e+02 +% +% The polar grids, on the other hand, are more complete in +% terms of latitude coverage: +% +% [lat,lon] = smapease2inverse(8999,8999,'N01') +% lat = 89.993669248945238 +% lon = -135 +% [lat,lon] = smapease2inverse(9000,9000,'N01') +% lat = 89.993669248945238 +% lon = 45 +% +% [lat,lon] = smapease2inverse(8999,8999,'S01') +% lat = -89.993669248945238 +% lon = -45 +% [lat,lon] = smapease2inverse(9000,9000,'S01') +% lat = -89.993669248945238 +% lon = 135 +% +% UPDATE North/south polar projections were added. (03/2012) +% +% REFERENCE Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. +% Savoie (2012): EASE-Grid 2.0: Incremental but Significant +% Improvements for Earth-Gridded Data Sets. ISPRS International +% Journal of Geo-Information, vol. 1, no. 1, pp. 32-45, +% http://www.mdpi.com/2220-9964/1/1/32/ +% +% Steven Chan, 11/2011 +% Email: steven.k.chan@jpl.nasa.gov + +function [lat,lon] = EASEv2_ind2latlon(row,col,gridid) + +% Constants returned by EASE2_GRID_INFO.PRO +projection = gridid(1); +switch gridid + case 'M36' + map_scale_m = 36032.220840584; + cols = 964; + rows = 406; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'M09' + map_scale_m = 9008.055210146; + cols = 3856; + rows = 1624; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'M03' + map_scale_m = 3002.6850700487; + cols = 11568; + rows = 4872; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'M01' + map_scale_m = 1000.89502334956; + cols = 34704; + rows = 14616; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N36' + map_scale_m = 36000.0; + cols = 500; + rows = 500; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N09' + map_scale_m = 9000.0; + cols = 2000; + rows = 2000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N03' + map_scale_m = 3000.0; + cols = 6000; + rows = 6000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N01' + map_scale_m = 1000.0; + cols = 18000; + rows = 18000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S36' + map_scale_m = 36000.0; + cols = 500; + rows = 500; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S09' + map_scale_m = 9000.0; + cols = 2000; + rows = 2000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S03' + map_scale_m = 3000.0; + cols = 6000; + rows = 6000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S01' + map_scale_m = 1000.0; + cols = 18000; + rows = 18000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + otherwise + disp(['ERROR: Incompatible grid specification.']); +end + +% Constants returned by EASE2_MAP_INFO.PRO +map_equatorial_radius_m = 6378137.0; +map_eccentricity = 0.081819190843; +e2 = map_eccentricity^2; +switch projection + case 'M' + map_reference_latitude = 0.0; + map_reference_longitude = 0.0; + map_second_reference_latitude = 30.0; + sin_phi1 = sin(map_second_reference_latitude*pi/180); + cos_phi1 = cos(map_second_reference_latitude*pi/180); + kz = cos_phi1/sqrt(1.0-e2*sin_phi1*sin_phi1); + case 'N' + map_reference_latitude = 90.0; + map_reference_longitude = 0.0; + case 'S' + map_reference_latitude = -90.0; + map_reference_longitude = 0.0; +end + +% Selected calculations inside WGS84_INVERSE.PRO +x = (col-r0)*map_scale_m; +y = (s0-row)*map_scale_m; + +% Selected calculations inside WGS84_INVERSE_XY.PRO +e4 = map_eccentricity^4; +e6 = map_eccentricity^6; +qp = (1.0-e2)*( (1.0/(1.0-e2))-(1.0/(2.0*map_eccentricity))*log((1.0-map_eccentricity)/(1.0+map_eccentricity))); +switch projection + case 'M' + beta = asin(2.0*y*kz/(map_equatorial_radius_m*qp)); + lam = x/(map_equatorial_radius_m*kz); + case 'N' + rho = sqrt(x.^2+y.^2); + beta = asin(1.0-(rho.^2)/(qp*map_equatorial_radius_m^2)); + lam = atan2(x,-y); + case 'S' + rho = sqrt(x.^2+y.^2); + beta = -asin(1.0-(rho.^2)/(qp*map_equatorial_radius_m^2)); + lam = atan2(x,y); +end +phi = beta+(((e2/3.0)+((31.0/180.0)*e4)+((517.0/5040.0)*e6))*sin(2.0*beta))+...\ + ((((23.0/360.0)*e4)+((251.0/3780.0)*e6))*sin(4.0*beta))+(((761.0/45360.0)*e6)*sin(6.0*beta)); +lat = phi*180.0/pi; +lon = map_reference_longitude+(lam*180.0/pi); +msk1 = lon < -180.0; lon(msk1) = lon(msk1) + 360.0; +msk2 = lon > +180.0; lon(msk2) = lon(msk2) - 360.0; +switch projection + case 'N' + idx = lat < 0.0; + lat(idx) = NaN; + lon(idx) = NaN; + case 'S' + idx = lat > 0.0; + lat(idx) = NaN; + lon(idx) = NaN; +end + +% ========= EOF ========================================================= diff --git a/GEOSldas_App/util/shared/matlab/EASEv2_latlon2ind.m b/GEOSldas_App/util/shared/matlab/EASEv2_latlon2ind.m index e8afdd7b..a4080d4b 100644 --- a/GEOSldas_App/util/shared/matlab/EASEv2_latlon2ind.m +++ b/GEOSldas_App/util/shared/matlab/EASEv2_latlon2ind.m @@ -1,218 +1,218 @@ -% -% SMAPEASE2FORWARD The principal function is to perform forward transformation -% from (lat,lon)'s to (row,col)'s for a set of nested EASE -% grids defined at 1, 3, 9, and 36km grid resolutions. These -% grids are all based on the EASE-Grid 2.0 specification (WGS84 -% ellipsoid). -% -% SYNTAX [row,col] = smapease2forward(lat,lon,gridid) -% -% where gridid is a 3-character string enclosed in single -% quotes, in the form of {M|N|S}{01,03,09,36}. This subroutine -% accepts vector inputs and produce vector outputs. -% -% HISTORY This subroutine was adapted from the offical EASE-Grid-2.0 -% conversion utilities (written in IDL) developed by the -% NSIDC. -% -% Note that in NSIDC's original implementation, (row,col) are -% zero-based. In other words, the first cell is (0,0) and the -% last cell is (N-1,M-1), where N and M are the row and column -% dimensions of the array. In this MATLAB implementation, the -% same convention is used. In other words, the end point of -% the first cell is located at (r,c) = (-0.5,-0.5) whereas the -% end point of the last cell is located at (r,c) = (14615.5, -% 34703.5). Thus, -% -% [lat,lon] = smapease2inverse(-0.5,-0.5,'M01') returns: -% lat = 85.044566407398861 -% lon = 1.799999999999994e+02 -% -% [lat,lon] = smapease2inverse(14615.5,34703.5,'M01') returns: -% lat = -85.044566407398861 -% lon = -1.799999999999994e+02 -% -% The polar grids, on the other hand, are more complete in -% terms of latitude coverage: -% -% [lat,lon] = smapease2inverse(8999,8999,'N01') -% lat = 89.993669248945238 -% lon = -135 -% [lat,lon] = smapease2inverse(9000,9000,'N01') -% lat = 89.993669248945238 -% lon = 45 -% -% [lat,lon] = smapease2inverse(8999,8999,'S01') -% lat = -89.993669248945238 -% lon = -45 -% [lat,lon] = smapease2inverse(9000,9000,'S01') -% lat = -89.993669248945238 -% lon = 135 -% -% UPDATE North/south polar projections were added. (03/2012) -% -% REFERENCE Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. -% Savoie (2012): EASE-Grid 2.0: Incremental but Significant -% Improvements for Earth-Gridded Data Sets. ISPRS International -% Journal of Geo-Information, vol. 1, no. 1, pp. 32-45, -% http://www.mdpi.com/2220-9964/1/1/32/ -% -% Steven Chan, 11/2011 -% Email: steven.k.chan@jpl.nasa.gov - -function [row,col] = EASEv2_latlon2ind(lat,lon,gridid,return_rounded) - -% By design, [row, col] are real numbers, with the fractional portion indicating -% the position of the specified [lat, lon] coordinates between adjacent grid -% cell centers. E.g., col=0.5 indicates that the input longitude is on the -% boundary between grid cells associated with col=0 and col=1. -% If the [optional] input argument 'return_rounded' is present and ~=0, then -% [row, col] are rounded to the nearest integer. - -% Constants returned by EASE2_GRID_INFO.PRO -projection = gridid(1); -switch gridid - case 'M36' - map_scale_m = 36032.220840584; - cols = 964; - rows = 406; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'M09' - map_scale_m = 9008.055210146; - cols = 3856; - rows = 1624; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'M03' - map_scale_m = 3002.6850700487; - cols = 11568; - rows = 4872; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'M01' - map_scale_m = 1000.89502334956; - cols = 34704; - rows = 14616; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N36' - map_scale_m = 36000.0; - cols = 500; - rows = 500; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N09' - map_scale_m = 9000.0; - cols = 2000; - rows = 2000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N03' - map_scale_m = 3000.0; - cols = 6000; - rows = 6000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'N01' - map_scale_m = 1000.0; - cols = 18000; - rows = 18000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S36' - map_scale_m = 36000.0; - cols = 500; - rows = 500; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S09' - map_scale_m = 9000.0; - cols = 2000; - rows = 2000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S03' - map_scale_m = 3000.0; - cols = 6000; - rows = 6000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - case 'S01' - map_scale_m = 1000.0; - cols = 18000; - rows = 18000; - r0 = (cols-1)/2; - s0 = (rows-1)/2; - otherwise - disp(['ERROR: Incompatible grid specification.']); -end - -% Constants returned by EASE2_MAP_INFO.PRO -epsilon = 1.0e-6; -map_equatorial_radius_m = 6378137.0; -map_eccentricity = 0.081819190843; -e2 = map_eccentricity^2; -switch projection - case 'M' - map_reference_latitude = 0.0; - map_reference_longitude = 0.0; - map_second_reference_latitude = 30.0; - sin_phi1 = sin(map_second_reference_latitude*pi/180); - cos_phi1 = cos(map_second_reference_latitude*pi/180); - kz = cos_phi1/sqrt(1.0-e2*sin_phi1*sin_phi1); - case 'N' - map_reference_latitude = 90.0; - map_reference_longitude = 0.0; - case 'S' - map_reference_latitude = -90.0; - map_reference_longitude = 0.0; -end - -% Selected calculations inside WGS84_CONVERT.PRO and WGS84_CONVERT_XY.PRO -dlon = lon - map_reference_longitude; -msk1 = dlon < -180.0; dlon(msk1) = dlon(msk1) + 360.0; -msk2 = dlon > +180.0; dlon(msk2) = dlon(msk2) - 360.0; -phi = lat*pi/180.0; -lam = dlon*pi/180.0; -sin_phi = sin(phi); -q = (1.0-e2)*((sin_phi./(1.0-e2*sin_phi.*sin_phi))-(1.0/(2.0*map_eccentricity))*log((1.0-map_eccentricity*sin_phi)./(1.0+map_eccentricity*sin_phi))); -qp = 1.0-((1.0-e2)/(2.0*map_eccentricity)*log((1.0-map_eccentricity)/(1.0+map_eccentricity))); -switch projection - case 'M' - x = map_equatorial_radius_m*kz*lam; - y = (map_equatorial_radius_m*q)/(2.0*kz); - case 'N' - tmp = qp - q; - tmp(abs(tmp) < epsilon) = 0.0; - rho = map_equatorial_radius_m*sqrt(tmp); - x = rho.*sin(lam); - y = -rho.*cos(lam); - case 'S' - tmp = qp + q; - tmp(abs(tmp) < epsilon) = 0.0; - rho = map_equatorial_radius_m*sqrt(tmp); - x = rho.*sin(lam); - y = rho.*cos(lam); -end -row = s0-(y/map_scale_m); -col = r0+(x/map_scale_m); -switch projection - case 'N' - idx = lat < 0.0; - row(idx) = NaN; - col(idx) = NaN; - case 'S' - idx = lat > 0.0; - row(idx) = NaN; - col(idx) = NaN; -end - -if exist('return_rounded','var') - if return_rounded - col=round(col); - row=round(row); - end -end - -% ========= EOF ========================================================= +% +% SMAPEASE2FORWARD The principal function is to perform forward transformation +% from (lat,lon)'s to (row,col)'s for a set of nested EASE +% grids defined at 1, 3, 9, and 36km grid resolutions. These +% grids are all based on the EASE-Grid 2.0 specification (WGS84 +% ellipsoid). +% +% SYNTAX [row,col] = smapease2forward(lat,lon,gridid) +% +% where gridid is a 3-character string enclosed in single +% quotes, in the form of {M|N|S}{01,03,09,36}. This subroutine +% accepts vector inputs and produce vector outputs. +% +% HISTORY This subroutine was adapted from the offical EASE-Grid-2.0 +% conversion utilities (written in IDL) developed by the +% NSIDC. +% +% Note that in NSIDC's original implementation, (row,col) are +% zero-based. In other words, the first cell is (0,0) and the +% last cell is (N-1,M-1), where N and M are the row and column +% dimensions of the array. In this MATLAB implementation, the +% same convention is used. In other words, the end point of +% the first cell is located at (r,c) = (-0.5,-0.5) whereas the +% end point of the last cell is located at (r,c) = (14615.5, +% 34703.5). Thus, +% +% [lat,lon] = smapease2inverse(-0.5,-0.5,'M01') returns: +% lat = 85.044566407398861 +% lon = 1.799999999999994e+02 +% +% [lat,lon] = smapease2inverse(14615.5,34703.5,'M01') returns: +% lat = -85.044566407398861 +% lon = -1.799999999999994e+02 +% +% The polar grids, on the other hand, are more complete in +% terms of latitude coverage: +% +% [lat,lon] = smapease2inverse(8999,8999,'N01') +% lat = 89.993669248945238 +% lon = -135 +% [lat,lon] = smapease2inverse(9000,9000,'N01') +% lat = 89.993669248945238 +% lon = 45 +% +% [lat,lon] = smapease2inverse(8999,8999,'S01') +% lat = -89.993669248945238 +% lon = -45 +% [lat,lon] = smapease2inverse(9000,9000,'S01') +% lat = -89.993669248945238 +% lon = 135 +% +% UPDATE North/south polar projections were added. (03/2012) +% +% REFERENCE Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. +% Savoie (2012): EASE-Grid 2.0: Incremental but Significant +% Improvements for Earth-Gridded Data Sets. ISPRS International +% Journal of Geo-Information, vol. 1, no. 1, pp. 32-45, +% http://www.mdpi.com/2220-9964/1/1/32/ +% +% Steven Chan, 11/2011 +% Email: steven.k.chan@jpl.nasa.gov + +function [row,col] = EASEv2_latlon2ind(lat,lon,gridid,return_rounded) + +% By design, [row, col] are real numbers, with the fractional portion indicating +% the position of the specified [lat, lon] coordinates between adjacent grid +% cell centers. E.g., col=0.5 indicates that the input longitude is on the +% boundary between grid cells associated with col=0 and col=1. +% If the [optional] input argument 'return_rounded' is present and ~=0, then +% [row, col] are rounded to the nearest integer. + +% Constants returned by EASE2_GRID_INFO.PRO +projection = gridid(1); +switch gridid + case 'M36' + map_scale_m = 36032.220840584; + cols = 964; + rows = 406; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'M09' + map_scale_m = 9008.055210146; + cols = 3856; + rows = 1624; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'M03' + map_scale_m = 3002.6850700487; + cols = 11568; + rows = 4872; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'M01' + map_scale_m = 1000.89502334956; + cols = 34704; + rows = 14616; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N36' + map_scale_m = 36000.0; + cols = 500; + rows = 500; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N09' + map_scale_m = 9000.0; + cols = 2000; + rows = 2000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N03' + map_scale_m = 3000.0; + cols = 6000; + rows = 6000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'N01' + map_scale_m = 1000.0; + cols = 18000; + rows = 18000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S36' + map_scale_m = 36000.0; + cols = 500; + rows = 500; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S09' + map_scale_m = 9000.0; + cols = 2000; + rows = 2000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S03' + map_scale_m = 3000.0; + cols = 6000; + rows = 6000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + case 'S01' + map_scale_m = 1000.0; + cols = 18000; + rows = 18000; + r0 = (cols-1)/2; + s0 = (rows-1)/2; + otherwise + disp(['ERROR: Incompatible grid specification.']); +end + +% Constants returned by EASE2_MAP_INFO.PRO +epsilon = 1.0e-6; +map_equatorial_radius_m = 6378137.0; +map_eccentricity = 0.081819190843; +e2 = map_eccentricity^2; +switch projection + case 'M' + map_reference_latitude = 0.0; + map_reference_longitude = 0.0; + map_second_reference_latitude = 30.0; + sin_phi1 = sin(map_second_reference_latitude*pi/180); + cos_phi1 = cos(map_second_reference_latitude*pi/180); + kz = cos_phi1/sqrt(1.0-e2*sin_phi1*sin_phi1); + case 'N' + map_reference_latitude = 90.0; + map_reference_longitude = 0.0; + case 'S' + map_reference_latitude = -90.0; + map_reference_longitude = 0.0; +end + +% Selected calculations inside WGS84_CONVERT.PRO and WGS84_CONVERT_XY.PRO +dlon = lon - map_reference_longitude; +msk1 = dlon < -180.0; dlon(msk1) = dlon(msk1) + 360.0; +msk2 = dlon > +180.0; dlon(msk2) = dlon(msk2) - 360.0; +phi = lat*pi/180.0; +lam = dlon*pi/180.0; +sin_phi = sin(phi); +q = (1.0-e2)*((sin_phi./(1.0-e2*sin_phi.*sin_phi))-(1.0/(2.0*map_eccentricity))*log((1.0-map_eccentricity*sin_phi)./(1.0+map_eccentricity*sin_phi))); +qp = 1.0-((1.0-e2)/(2.0*map_eccentricity)*log((1.0-map_eccentricity)/(1.0+map_eccentricity))); +switch projection + case 'M' + x = map_equatorial_radius_m*kz*lam; + y = (map_equatorial_radius_m*q)/(2.0*kz); + case 'N' + tmp = qp - q; + tmp(abs(tmp) < epsilon) = 0.0; + rho = map_equatorial_radius_m*sqrt(tmp); + x = rho.*sin(lam); + y = -rho.*cos(lam); + case 'S' + tmp = qp + q; + tmp(abs(tmp) < epsilon) = 0.0; + rho = map_equatorial_radius_m*sqrt(tmp); + x = rho.*sin(lam); + y = rho.*cos(lam); +end +row = s0-(y/map_scale_m); +col = r0+(x/map_scale_m); +switch projection + case 'N' + idx = lat < 0.0; + row(idx) = NaN; + col(idx) = NaN; + case 'S' + idx = lat > 0.0; + row(idx) = NaN; + col(idx) = NaN; +end + +if exist('return_rounded','var') + if return_rounded + col=round(col); + row=round(row); + end +end + +% ========= EOF ========================================================= diff --git a/GEOSldas_App/util/shared/python/EASEv2.py b/GEOSldas_App/util/shared/python/EASEv2.py new file mode 100644 index 00000000..31cd4392 --- /dev/null +++ b/GEOSldas_App/util/shared/python/EASEv2.py @@ -0,0 +1,199 @@ +import numpy as np + +# SMAPEASE2INVERSE The principal function is to perform inverse transformation +# from (row,col)'s to (lat,lon)'s for a set of nested EASE +# grids defined at 1, 3, 9, and 36km grid resolutions. These +# grids are all based on the EASE-Grid 2.0 specification (WGS84 +# ellipsoid). + +# SYNTAX [lat,lon] = smapease2inverse(row,col,gridid) + +# where gridid is a 3-character string enclosed in single +# quotes, in the form of {M|N|S}{01,03,09,36}. This subroutine +# accepts vector inputs and produce vector outputs. + +# HISTORY This subroutine was adapted from the offical EASE-Grid-2.0 +# conversion utilities (written in IDL) developed by the +# NSIDC. + +# Note that in NSIDC's original implementation, (row,col) are +# zero-based. In other words, the first cell is (0,0) and the +# last cell is (N-1,M-1), where N and M are the row and column +# dimensions of the array. In this MATLAB implementation, the +# same convention is used. In other words, the end point of +# the first cell is located at (r,c) = (-0.5,-0.5) whereas the +# end point of the last cell is located at (r,c) = (14615.5, +# 34703.5). Thus, + +# [lat,lon] = smapease2inverse(-0.5,-0.5,'M01') returns: +# lat = 85.044566407398861 +# lon = 1.799999999999994e+02 + +# [lat,lon] = smapease2inverse(14615.5,34703.5,'M01') returns: +# lat = -85.044566407398861 +# lon = -1.799999999999994e+02 + +# The polar grids, on the other hand, are more complete in +# terms of latitude coverage: + +# [lat,lon] = smapease2inverse(8999,8999,'N01') +# lat = 89.993669248945238 +# lon = -135 +# [lat,lon] = smapease2inverse(9000,9000,'N01') +# lat = 89.993669248945238 +# lon = 45 + +# [lat,lon] = smapease2inverse(8999,8999,'S01') +# lat = -89.993669248945238 +# lon = -45 +# [lat,lon] = smapease2inverse(9000,9000,'S01') +# lat = -89.993669248945238 +# lon = 135 + +# UPDATE North/south polar projections were added. (03/2012) + +# REFERENCE Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. +# Savoie (2012): EASE-Grid 2.0: Incremental but Significant +# Improvements for Earth-Gridded Data Sets. ISPRS International +# Journal of Geo-Information, vol. 1, no. 1, pp. 32-45, +# http://www.mdpi.com/2220-9964/1/1/32/ +# +# Steven Chan, 11/2011 +# Email: steven.k.chan@jpl.nasa.gov + +def EASEv2_ind2latlon(row=None,col=None,gridid=None,*args,**kwargs): + + # Constants returned by EASE2_GRID_INFO.PRO + projection=gridid[0] + if 'M36' == gridid: + map_scale_m=36032.220840584 + cols=964 + rows=406 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'M09' == gridid: + map_scale_m=9008.055210146 + cols=3856 + rows=1624 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'M03' == gridid: + map_scale_m=3002.6850700487 + cols=11568 + rows=4872 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'M01' == gridid: + map_scale_m=1000.89502334956 + cols=34704 + rows=14616 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'N36' == gridid: + map_scale_m=36000.0 + cols=500 + rows=500 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'N09' == gridid: + map_scale_m=9000.0 + cols=2000 + rows=2000 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'N03' == gridid: + map_scale_m=3000.0 + cols=6000 + rows=6000 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'N01' == gridid: + map_scale_m=1000.0 + cols=18000 + rows=18000 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'S36' == gridid: + map_scale_m=36000.0 + cols=500 + rows=500 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'S09' == gridid: + map_scale_m=9000.0 + cols=2000 + rows=2000 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'S03' == gridid: + map_scale_m=3000.0 + cols=6000 + rows=6000 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + elif 'S01' == gridid: + map_scale_m=1000.0 + cols=18000 + rows=18000 + r0=(cols - 1) / 2 + s0=(rows - 1) / 2 + else: + print('ERROR: Incompatible grid specification.') + + # Constants returned by EASE2_MAP_INFO.PRO + map_equatorial_radius_m=6378137.0 + map_eccentricity=0.081819190843 + e2=map_eccentricity ** 2 + if 'M' == projection: + map_reference_latitude=0.0 + map_reference_longitude=0.0 + map_second_reference_latitude=30.0 + sin_phi1=np.sin(map_second_reference_latitude*np.pi / 180) + cos_phi1=np.cos(map_second_reference_latitude*np.pi / 180) + kz=cos_phi1 / np.sqrt(1.0 - e2*sin_phi1*sin_phi1) + elif 'N' == projection: + map_reference_latitude=90.0 + map_reference_longitude=0.0 + elif 'S' == projection: + map_reference_latitude=- 90.0 + map_reference_longitude=0.0 + + # Selected calculations inside WGS84_INVERSE.PRO + x= (col - r0) * map_scale_m + y= (s0 - row) * map_scale_m + # Selected calculations inside WGS84_INVERSE_XY.PRO + e4=map_eccentricity ** 4 + e6=map_eccentricity ** 6 + qp= (1.0 - e2)*((1.0 / (1.0 - e2)) - (1.0 / (2.0*map_eccentricity))*np.log((1.0 - map_eccentricity) / (1.0 + map_eccentricity))) + if 'M' == projection: + beta= np.arcsin(2.0 * y * kz / (map_equatorial_radius_m * qp)) + lam= x / (map_equatorial_radius_m*kz) + elif 'N' == projection: + rho=np.sqrt(x ** 2 + y ** 2) + beta=np.arcsin(1.0 - (rho ** 2) / (qp * map_equatorial_radius_m ** 2)) + lam=np.arctan2(x,- y) + elif 'S' == projection: + rho=np.sqrt(x ** 2 + y ** 2) + beta=- np.arcsin(1.0 - (rho ** 2) / (qp * map_equatorial_radius_m ** 2)) + lam=np.arctan2(x,y) + + phi= beta + (((e2 / 3.0) + ((31.0 / 180.0) * e4) + ((517.0 / 5040.0) * e6)) * np.sin(2.0 * beta)) + \ + ((((23.0 / 360.0) * e4) + ((251.0 / 3780.0) * e6)) * np.sin(4.0*beta)) + (((761.0 / 45360.0) * e6) * np.sin(6.0 * beta)) + lat= phi * 180.0 / np.pi + lon= map_reference_longitude + (lam * 180.0 / np.pi) + msk1= np.where(lon < - 180.0) + lon[msk1]= lon[msk1] + 360.0 + msk2= np.where(lon > 180.0) + lon[msk2]=lon[msk2] - 360.0 + if 'N' == projection: + idx=np.where(lat < 0.0) + lat[idx]= float('nan') + lon[idx]= float('nan') + elif 'S' == projection: + idx= np.where(lat > 0.0) + lat[idx]= float('nan') + lon[idx]= float('nan') + + return lat, lon + +# ============= EOF ============================================ diff --git a/GEOSldas_App/util/shared/python/plot.py b/GEOSldas_App/util/shared/python/plot.py new file mode 100644 index 00000000..5dff0892 --- /dev/null +++ b/GEOSldas_App/util/shared/python/plot.py @@ -0,0 +1,80 @@ +# collection of py functions for plotting + +import numpy as np +import matplotlib.pyplot as plt + +from mpl_toolkits import basemap + +def plotMap( + data, + *, + ax=None, + lat=None, + lon=None, + title=None, + cRange=None, + figsize=(8, 4), + clbar=True, + cRangeint=False, + cmap=plt.cm.jet, + bounding=None, + prj="cyl", +): + + # color range + if cRange is not None: + vmin = cRange[0] + vmax = cRange[1] + else: + temp = flatData(data) + vmin = np.percentile(temp, 5) + vmax = np.percentile(temp, 95) + if cRangeint is True: + vmin = int(round(vmin)) + vmax = int(round(vmax)) + if ax is None: + fig = plt.figure(figsize=figsize) + ax = fig.subplots() + + # map boundary + if bounding is None: + bounding = [ + np.min(lat) - 0.5, + np.max(lat) + 0.5, + np.min(lon) - 0.5, + np.max(lon) + 0.5, + ] + + # add basemap + mm = basemap.Basemap( + llcrnrlat=bounding[0], + urcrnrlat=bounding[1], + llcrnrlon=bounding[2], + urcrnrlon=bounding[3], + projection=prj, + resolution="c", + ax=ax, + ) + mm.drawcoastlines() + #mm.drawstates(linestyle="dashed") + #mm.drawcountries(linewidth=1.0, linestyle="-.") + # plot data on basemap + cs = mm.pcolormesh(lon, lat, data, cmap=cmap, vmin=vmin, vmax=vmax) + + # colorbar + if clbar is True: + cb = mm.colorbar(cs, pad="5%", location="bottom") + if 'normalized' in title: + cb.set_ticks(np.linspace(vmin,vmax,6)) + cb.set_ticklabels([f'{10**x:.2f}' for x in np.linspace(vmin, vmax, 6)]) + + # plot title, return objects in case plot needs adjustment after function call + if title is not None: + ax.set_title(title) + if ax is None: + return fig, ax, mm + else: + return mm, cs + + +# ================ EOF ================================================= diff --git a/GEOSldas_App/util/shared/python/read_GEOSldas.py b/GEOSldas_App/util/shared/python/read_GEOSldas.py new file mode 100644 index 00000000..96ba4173 --- /dev/null +++ b/GEOSldas_App/util/shared/python/read_GEOSldas.py @@ -0,0 +1,239 @@ + +# collection of readers for GEOSldas output files + +import numpy as np +import struct +import os +import numpy as np + +# ---------------------------------------------------------------------------- +# +# reader for GEOSldas obsparam file (ASCII) + +def read_obs_param(fname): + print(f"Reading {fname}") + + with open(fname, 'r') as fid: + N_obs_param = int(fid.readline().strip()) + + obs_param = [] + for _ in range(N_obs_param): + param = {} + param['descr'] = fid.readline().strip().strip('"') + param['species'] = float(fid.readline().strip()) + param['orbit'] = float(fid.readline().strip()) + param['pol'] = float(fid.readline().strip()) + + param['N_ang'] = int(float(fid.readline().strip())) + + param['ang'] = np.array([float(x) for x in fid.readline().split()]) + + param['freq'] = float(fid.readline().strip()) + param['FOV'] = float(fid.readline().strip()) + param['FOV_units'] = fid.readline().strip().strip('"') + param['assim'] = fid.readline().strip() + param['scale'] = fid.readline().strip() + param['getinnov'] = fid.readline().strip() + param['RTM_ID'] = float(fid.readline().strip()) + param['bias_Npar'] = float(fid.readline().strip()) + param['bias_trel'] = float(fid.readline().strip()) + param['bias_tcut'] = float(fid.readline().strip()) + param['nodata'] = float(fid.readline().strip()) + param['varname'] = fid.readline().strip().strip('"') + param['units'] = fid.readline().strip().strip('"') + param['path'] = fid.readline().strip().strip('"') + param['name'] = fid.readline().strip().strip('"') + param['maskpath'] = fid.readline().strip().strip('"') + param['maskname'] = fid.readline().strip().strip('"') + param['scalepath'] = fid.readline().strip().strip('"') + param['scalename'] = fid.readline().strip().strip('"') + param['flistpath'] = fid.readline().strip().strip('"') + param['flistname'] = fid.readline().strip().strip('"') + param['errstd'] = float(fid.readline().strip()) + param['std_normal_max'] = float(fid.readline().strip()) + param['zeromean'] = fid.readline().strip() + param['coarsen_pert'] = fid.readline().strip() + param['xcorr'] = float(fid.readline().strip()) + param['ycorr'] = float(fid.readline().strip()) + param['adapt'] = float(fid.readline().strip()) + + obs_param.append(param) + + print(f"Done reading obs_param for {N_obs_param} species") + + return obs_param + +# ---------------------------------------------------------------------------- +# +# reader for GEOSldas tilecoord file (binary) + +def read_tilecoord(fname): + int_precision = 'i' + float_precision = 'f' + + # SPECIFY endianness + machfmt = '<' # '>' for big-endian, '<' for little-endian + + print(f"reading from {fname}") + + tile_coord = {} + + with open(fname, 'rb') as ifp: + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + tile_coord['N_tile'] = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + Nt = tile_coord['N_tile'] + + fields = ['tile_id', 'typ', 'pfaf', 'com_lon', 'com_lat', 'min_lon', 'max_lon', + 'min_lat', 'max_lat', 'i_indg', 'j_indg', 'frac_cell', 'frac_pfaf', + 'area', 'elev'] + + for field in fields: + this_dtype = int_precision if field in ['tile_id', 'typ', 'pfaf', 'i_indg', 'j_indg'] else float_precision + + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + tile_coord[field] = np.frombuffer(ifp.read(Nt * 4), dtype=f'{machfmt}{this_dtype}') + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + print("done reading file") + return tile_coord + +# ---------------------------------------------------------------------------- +# +# reader for GEOSldas ObsFcstAna file (binary) + +def read_ObsFcstAna(fname, isLDASsa=False): + + # Initialize outputs + nodata = -9999 + + date_time = { + 'year' : nodata, + 'month' : nodata, + 'day' : nodata, + 'hour' : nodata, + 'min' : nodata, + 'sec' : nodata, + 'dofyr' : nodata, + 'pentad': nodata + } + + obs_assim = [] + obs_species = [] + obs_tilenum = [] + obs_lon = [] + obs_lat = [] + obs_obs = [] + obs_obsvar = [] + obs_fcst = [] + obs_fcstvar = [] + obs_ana = [] + obs_anavar = [] + + # SPECIFY endianness + machfmt = '<' # '>' for big-endian, '<' for little-endian + + if os.path.exists(fname): + print(f"reading from {fname}") + + with open(fname, 'rb') as ifp: + # Read N_obs and time stamp entry + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + N_obs = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + year, month, day, hour, minute, second, dofyr, pentad = struct.unpack(f'{machfmt}8i', ifp.read(32)) + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + date_time = { + 'year' : year, + 'month' : month, + 'day' : day, + 'hour' : hour, + 'min' : minute, + 'sec' : second, + 'dofyr' : dofyr, + 'pentad': pentad + } + + # Read observation assim flag + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + tmp_data = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}i').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + obs_assim = np.zeros(N_obs, dtype=int) + obs_assim[tmp_data != 0] = 1 + + # Read species information + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_species = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}i').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read tile number information + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_tilenum = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}i').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read longitude + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_lon = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read latitude + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_lat = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read observation value + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_obs = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read observation variance + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_obsvar = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read observation-space model forecast value + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_fcst = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read observation-space model forecast variance + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_fcstvar = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read observation-space analysis value + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_ana = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # Read observation-space analysis variance + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + obs_anavar = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy() + fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0] + + # No-data check + obs_obsvar[ obs_obsvar == nodata] = np.nan + obs_fcst[ obs_fcst == nodata] = np.nan + obs_fcstvar[obs_fcstvar == nodata] = np.nan + obs_ana[ obs_ana == nodata] = np.nan + obs_anavar[ obs_anavar == nodata] = np.nan + + else: + print(f"file does not exist: {fname}") + + return {'date_time' : date_time, + 'obs_assim' : obs_assim, + 'obs_species': obs_species, + 'obs_tilenum': obs_tilenum, + 'obs_lon' : obs_lon, + 'obs_lat' : obs_lat, + 'obs_obs' : obs_obs, + 'obs_obsvar' : obs_obsvar, + 'obs_fcst' : obs_fcst, + 'obs_fcstvar': obs_fcstvar, + 'obs_ana' : obs_ana, + 'obs_anavar' : obs_anavar} + +# ================ EOF ================================================= diff --git a/GEOSldas_App/util/shared/python/remap_1d_to_2d.py b/GEOSldas_App/util/shared/python/remap_1d_to_2d.py new file mode 100644 index 00000000..d0a1d9fc --- /dev/null +++ b/GEOSldas_App/util/shared/python/remap_1d_to_2d.py @@ -0,0 +1,29 @@ +import os +import numpy as np + +# convert data from 1-d array ("compressed") format to 2-d array (grid) format based on lat/lon coordinates +# +# 1-d input data can have one additional dimension (e.g., time) +# +# 2-d output data is lat-by-lon[-by-time] + +def remap_1d_to_2d( data_1d, *, lat_1d, lon_1d): + + # Extract unique coordinates + unique_lats, indY = np.unique(lat_1d, return_inverse=True) + unique_lons, indX = np.unique(lon_1d, return_inverse=True) + + ny = len(unique_lats) + nx = len(unique_lons) + + if data_1d.ndim == 2: #[n_1d, ntime] + ntime = data_1d.shape[1] + data_2d = np.full([ny, nx, ntime], np.nan) + data_2d[indY, indX, :] = data_1d + elif data_1d.ndim == 1: #[n_1d] + data_2d = np.full([ny, nx], np.nan) + data_2d[indY, indX ] = data_1d + + return data_2d, unique_lats, unique_lons + +