diff --git a/CHANGELOG.md b/CHANGELOG.md index 969914d..6a9d43e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added Matlab and python readers for binary Tb scaling parameters files. - Added python reader for binary catparam files. - Added QC of SMAP L1C_TB using max value for Tb_error. +- Added SMOS Tb preprocessing scripts. ### Changed diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/.gitignore b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/.gitignore new file mode 100644 index 0000000..c201240 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/.gitignore @@ -0,0 +1,4 @@ +__pycache__/ +*.pyc +out.log +*.swp diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/README.md b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/README.md new file mode 100644 index 0000000..a47b975 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/README.md @@ -0,0 +1,65 @@ +# SMOS_preproc + +## 1. Purpose +Preprocesses SMOS L1C brightness temperature data for near-real-time (NRT) +ingestion into the MERRA21C-land pipeline: + +``` +EE .zip --(smos-ee-to-nc.sh)--> NetCDF --(preprocess_nc)--> EASEv2 M36 "REG" binaries --(SCLF1C_reg2fit)--> Tb40 "FIT" binaries +``` + +Normal operation is a **daily cron job** that processes the previous day's +data. It also supports on-demand backfill of a single day or a date range. + +## 2. Requirements +- Runs on Discover/NCCS. +- Environment: `module load python/GEOSpyD` (provides numpy, scipy, netCDF4, + PyYAML — no other Python environment should be needed). +- External script dependency: `smos-ee-to-nc.sh` (path set in `config.yaml`, + currently maintained under `/discover/nobackup/projects/gmao/smap/...` — + **not part of this repo**, +- Static input: `data/GEOSIT_to_EASEv2_M36.mat` — pre-generated EASEv2 + regridding weights, link to the actual location on Discover +- Static Aux files: `data/SM_OPER_AUX_GAL_SM_20050101T000000_20500101T000000_001_003_3` + link to the actual location on Discover. + +## 3. Configuration (`config.yaml`) +| Key | Meaning | +|-----------------------------|-------------------------------------------------------------------------------------------------| +| `ee_to_nc_script` | Path to the external `smos-ee-to-nc.sh` converter. | +| `smos_base_path` | Where incoming SMOS `SM_*_MIR_SCLF1C_*.zip` (EE) files land, organized `/Y/M/`. | +| `tmp_nc_path` | Scratch directory for converted NetCDF files. A `to_delete/` subfolder is created here. Files therein are safe to delete after completion. | +| `out_reg_path` | Output root for REG binaries (`SMOS_reg_Tb_*.bin`, organized by `//`). FIT binaries are written to a sibling directory: `_reg_` in this path is replaced by `_fit_` and further nested under `SMOS_fit_poly2//`. | +| `GEOSIT_path` | Path to GEOS-IT data. | +| `GEOSIT_to_EASEv2_M36_file` | Path and name of file mapping from GEOS-IT output grid to EASEv2_M36 grid. | +| `SM_OPER_AUX_GAL_SM_path` | Path to SMOS galaxy correction files. | + + +Current values are set for the production Discover paths under +`/discover/nobackup/projects/gmao/smap/SMAP_Nature/SMOS/...` and +`/discover/nobackup/dao_ops/SMOS/AOSMOS.4662/SCLF1C/`. + +## 4. Usage +``` +module load python/GEOSpyD + +python SMOSproc_main.py # NRT: process yesterday's data (normal cron mode) +python SMOSproc_main.py --date 20260801 # Backfill a single day +python SMOSproc_main.py --start 20260801 --end 20260803 # Backfill an inclusive date range +``` +`--date` and `--start`/`--end` are mutually exclusive; `--end` requires +`--start`. + +## 5. Layout +``` +SMOSproc_main.py top-level driver +config.yaml paths +data/ static regridding weights (GEOSIT_to_EASEv2_M36.mat) and AUX +src/preprocess_nc.py NetCDF -> EASEv2 M36 REG binaries +src/SCLF1C_reg2fit.py REG -> Tb40 FIT binaries +src/readwrite/ I/O helpers (SMOS NetCDF, REG binary, GEOS-IT, aux Gal SM) +src/helper/ grid/geometry/time utilities (EASEv2 indexing, h/v tile + conversion, celestial angle calc, galactic+atmospheric correction) +``` + + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/SMOSproc_main.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/SMOSproc_main.py new file mode 100644 index 0000000..19bb729 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/SMOSproc_main.py @@ -0,0 +1,164 @@ +import os +import shutil +import logging +import multiprocessing +import yaml + +import sys; sys.path.append('../../../shared/python/') + +from src import preprocess_nc, SCLF1C_reg2fit +from src.helper.util import * + +# Configure logging +logging.basicConfig( + level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + datefmt='%Y-%m-%d %H:%M:%S' +) + +# --------------------------------------------------------- +# Configuration +# --------------------------------------------------------- +CONFIG_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)),'config.yaml') +with open(CONFIG_PATH) as f: + config = yaml.safe_load(f)['paths'] + +EE_TO_NC_SCRIPT = config['ee_to_nc_script'] +SMOS_BASE_PATH = config['smos_base_path'] +TMP_NC_PATH = config['tmp_nc_path'] +OUT_REG_PATH = config['out_reg_path'] + +def run_in_isolated_process(func, *args): + """ + Runs a function in an isolated child process and waits for it to finish. + This prevents memory leaks or state changes in external libraries from + affecting the main script. + """ + p = multiprocessing.Process(target=func, args=args) + p.start() + p.join() + + if p.exitcode != 0: + logging.error(f"Process running {func.__name__} failed with exit code {p.exitcode}") + raise RuntimeError(f"{func.__name__} failed during execution.") + + logging.info(f"Successfully completed {func.__name__}.") + +def process_ee_to_nc(date_time: datetime) -> list: + """ + Finds .zip ee files for the given date, converts them to .nc if needed, + and returns a list of resulting netcdf files. + """ + # Build search pattern for EE files + year, month, day = date_time.strftime('Y%Y'), date_time.strftime('M%m'), date_time.strftime('%d') + date_str = date_time.strftime('%Y%m%d') + + search_pattern = os.path.join( + SMOS_BASE_PATH, year, month, + f'SM_*_MIR_SCLF1C_{date_str}*_724_*zip' + ) + + eeflist = sorted(glob.glob(search_pattern)) + + logging.info(f"[{date_str}] Found {len(eeflist)} zip files to process.") + + if not eeflist: + raise RuntimeError(f"No SMOS EE zip files found for {date_str}") + + # Convert EE to NC + for fee in eeflist: + base_name = os.path.basename(fee)[:-3] + 'nc' + target_nc_file = os.path.join(TMP_NC_PATH, base_name) + + if not os.path.isfile(target_nc_file): + logging.info(f"Running ee_to_nc conversion on: {fee}") + result = subprocess.run([ + EE_TO_NC_SCRIPT, + '--target-directory', TMP_NC_PATH, + fee + ], capture_output=True, text=True) + + if result.returncode != 0: + logging.error(f"Conversion script failed for {fee}:\n{result.stderr}") + raise RuntimeError("ee_to_nc conversion corrupted") + + logging.info(f"Done ee_to_nc conversion: {fee}") + + # Return list of resulting NC files + nc_search_pattern = os.path.join(TMP_NC_PATH, f'SM_*_MIR_SCLF1C_{date_str}*724*.nc') + ncflist = sorted(glob.glob(nc_search_pattern)) + + if len(eeflist) != len(ncflist): + logging.warning(f"File count mismatch: {len(eeflist)} zip files vs {len(ncflist)} nc files.") + + return ncflist + +def main(): + args = parse_args() + start_time, end_time = get_time_range(args) + logging.info( + f"Processing {start_time:%Y-%m-%d} through " + f"{(end_time - timedelta(days=1)):%Y-%m-%d}" + ) + + # Ensure necessary directories exist + os.makedirs(os.path.join(TMP_NC_PATH, 'to_delete'), exist_ok=True) + os.makedirs(OUT_REG_PATH, exist_ok=True) + + current_date = start_time + + while current_date < end_time: + date_str = current_date.strftime('%Y%m%d') + logging.info(f"=== Starting processing for {date_str} ===") + + # Step 1: Convert files to NetCDF + try: + ncflist = process_ee_to_nc(current_date) + except RuntimeError as e: + logging.error(f"Halting processing for {date_str} due to error: {e}") + current_date += timedelta(days=1) + continue + + # Step 2: Preprocess NetCDF files into REG + for fnc in ncflist: + logging.info(f"Preprocessing NC file: {fnc}") + try: + run_in_isolated_process(preprocess_nc, fnc, config) + except RuntimeError as e: + logging.error(f"Skipping {fnc} due to error: {e}") + continue + + # Move processed file to 'to_delete' + dest = os.path.join(TMP_NC_PATH, 'to_delete', os.path.basename(fnc)) + shutil.move(fnc, dest) + logging.info(f"Moved {fnc} to cleanup directory.") + + # Step 3: Run REG to FIT processing for Ascending and Descending + next_date = current_date + timedelta(days=1) + + logging.info("Running SCLF1C_reg2fit for Ascending (_A)") + try: + run_in_isolated_process(SCLF1C_reg2fit, OUT_REG_PATH, current_date, next_date, '_A') + except RuntimeError as e: + logging.error(f"SCLF1C_reg2fit failed for Ascending (_A) on {date_str}: {e}") + + logging.info("Running SCLF1C_reg2fit for Descending (_D)") + try: + run_in_isolated_process(SCLF1C_reg2fit, OUT_REG_PATH, current_date, next_date, '_D') + except RuntimeError as e: + logging.error(f"SCLF1C_reg2fit failed for Descending (_D) on {date_str}: {e}") + + logging.info(f"=== Completed processing to Tb40 for {date_str} ===") + + # Advance to the next day + current_date += timedelta(days=1) + + +if __name__ == '__main__': + # Set the multiprocessing start method exactly once here + try: + multiprocessing.set_start_method('spawn') + except RuntimeError: + pass # Context was already set + + main() diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/config.yaml b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/config.yaml new file mode 100644 index 0000000..4b2be1b --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/config.yaml @@ -0,0 +1,8 @@ +paths: + ee_to_nc_script: /discover/nobackup/projects/gmao/smap/SMAP_Nature/SMOS/smos-ee-to-netcdf/smos-ee-to-nc.sh + smos_base_path: /discover/nobackup/dao_ops/SMOS/AOSMOS.4662/SCLF1C/ + tmp_nc_path: /discover/nobackup/projects/gmao/smap/SMAP_Nature/SMOS/netcdf/ + out_reg_path: /discover/nobackup/projects/gmao/smap/SMAP_Nature/SMOS/EASEv2/ESA/SMOS_M36_SCLF1C_reg_nosky_noatm_v724_ESA_v102/ + GEOSIT_path: /discover/nobackup/projects/gmao/geos-it/dao_ops/archive/ + GEOSIT_to_EASEv2_M36_file: /discover/nobackup/projects/gmao/smap/SMAP_Nature/SMOS/Regrid_data/GEOSIT_to_EASEv2_M36.mat + SM_OPER_AUX_GAL_SM_path: /discover/nobackup/projects/gmao/smap/SMAP_Nature/SMOS/SMOS_ESA/AUX_DATA/SM_OPER_AUX_GAL_SM_20050101T000000_20500101T000000_001_003_3/ diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/SCLF1C_reg2fit.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/SCLF1C_reg2fit.py new file mode 100644 index 0000000..1192db0 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/SCLF1C_reg2fit.py @@ -0,0 +1,193 @@ +from datetime import datetime,timedelta +import numpy as np +import os +import sys +from src.readwrite import read_bin_SMOS_reg, write_bin_SMOS_reg + +def SCLF1C_reg2fit(reg_path,start_time,end_time,Asc): + #Read preprocessed 'regular' SMOS Tb data - all angles + #Interpolation / fit -> extract 40^o angle, and all other angles + #======================================================== + + in_path=reg_path + out_path=reg_path.replace('_reg_','_fit_') + #file_in_prefix = 'SMOS_reg_nosky_noatm_Tb_'; + #file_out_prefix = 'SMOS_fit_nosky_noatm_Tb_'; + file_in_prefix='SMOS_reg_Tb_' + file_out_prefix='SMOS_fit_Tb_' + + # half a hour + dtstep=1800 + #---------------------------------- + TbHVconstrained='' + #TbHVconstrained = '_TbHVconstrained'; #constrained fitting of multiple Tb-variables, or empty if individual + + #EXPONENTIAL fit: y = b0 + b1.exp(-b3.x) + #tunable coefficients appear non-linearly -> non-linear fitting takes iteration time + #---------------------------------- + #fit_type = 'exp1'; + + #QUADRATIC fit: y =ax^2 + bx + c + #tunable coefficients appear linearly -> quick matrix inversion + #---------------------------------- + fit_type='poly2' + + read_ind_latlon='latlon' + write_ind_latlon='latlon' + + overwrite=1 + + #------------------------------------------------------ + # QC + #------------------------------------------------------ + + N_gt30_lt50=10 #number of data at angles>30 and <50^o + + N_points_tot=15 #total number of angles used in the fit + + min_N_ang=5 #request at least 5 angles available at either side + #of each interpolated value + + #remove data that indicate of RFI: + cutoff_weight=2 / (7 ** 2) #eliminate data with excessive Tb variability + + cutoff_K=320 + + #------------------------------------------------------ + + N_out_fields=16 + field_tags=['Tb_H','Tb_V','std_h','std_v','c_h','c_v','RAh','RAv',\ + 'Tb3','Tb4','std_3','std_4','c_3','c_4','RA3','RA4'] + + #Variables to be fitted + #------------------------- + #Tb_ind = [0 1 8 9]; + #We really do not know the shape of T3, T4 - do not try to fit for now + + Tb_ind=[0,1] #Tbh, Tbv + + # std_ind = Tb_ind + 2; + # c_ind = Tb_ind + 4; + # Ra_ind = Tb_ind + 6; + + data_product='SCLF1C' + + date_time_new=start_time + + no_data_value=-999.0 + no_data_tol=0.0001 + + # Go through loop of orbits + while date_time_new < end_time: + # augment date_time and t_ind + date_time_old=date_time_new + date_time_new=date_time_old + timedelta(seconds=dtstep) + + fname = in_path+date_time_new.strftime('%Y%m')+'/'+ \ + file_in_prefix+date_time_new.strftime('%Y%m%d_%H%M')+\ + Asc+'.bin' + + if not os.path.exists(out_path+'/SMOS_fit_'+fit_type+TbHVconstrained+'/'+ \ + date_time_new.strftime('%Y%m')): + os.makedirs(out_path+'/SMOS_fit_'+fit_type+TbHVconstrained+'/'+ \ + date_time_new.strftime('%Y%m')) + + if os.path.isfile(fname): + out_filename=out_path+'/SMOS_fit_'+fit_type+TbHVconstrained+'/'+ \ + date_time_new.strftime('%Y%m')+'/'+file_out_prefix+ \ + date_time_new.strftime('%Y%m%d_%H%M')+Asc+'.bin' + + time=date_time_new.strftime('%Y%m%d%H%M') + + data,inc_ang,col_ind,row_ind,asc_flag,version,prep_version,o_start_time,o_end_time,N_grid=\ + read_bin_SMOS_reg(fname,N_out_fields,read_ind_latlon,data_product,nargout=10) + + #inc_ang=inc_ang.T + data_out=data.copy() + + if len(data_out.shape) < 3: + continue + + for v in range(len(Tb_ind)): + data_out[Tb_ind[v],:,:]=np.nan + data_out[Tb_ind[v] + 2,:,:]=np.nan + data_out[Tb_ind[v] + 4,:,:]=np.nan + data_out[Tb_ind[v] + 6,:,:]=np.nan + #-------------------------------------------------------------------- + #weighted interpolation + #weight = number of points in average / stdv + #-------------------------------------------------------------------- + data[np.abs(data - no_data_value) < no_data_tol]=np.nan + + #for all grid cells, overwrite the original 'data' + #with new 'data' containing fitted information + for i in range(len(col_ind)): + if (~np.isnan(data[:,i,:])).any(): + if TbHVconstrained == '': + #for each of the 4 or 2 Tb-variables individually + for v in range(len(Tb_ind)): + #1) actual Tb + Tb_data=data[Tb_ind[v],i,:] + Tb_data[Tb_data < 0]=np.nan + + #stdv+RA + #total_prior_uncert = sqrt(squeeze(data(Tb_ind(v)+2,i,:)).^2 + squeeze(data(Tb_ind(v)+6,i,:)).^2); + total_prior_uncert=data[Tb_ind[v]+2,i,:]**2 + \ + data[Tb_ind[v]+6,i,:]**2 + #Each angle has a different number of contributing DGG cells, + #because it is collecting data over a [x-0.5 x+0.5]-window and + #from different instants (seconds) in time (fwd, backward looking) + weights=np.sqrt(data[Tb_ind[v]+4,i,:] / total_prior_uncert) + good=np.nonzero(~np.isnan(Tb_data) * (Tb_data < cutoff_K) * \ + (weights>cutoff_weight) * ~np.isnan(weights) )[0] + + Tb_data=Tb_data[good] + weights=weights[good] + x=inc_ang[good] + + if ((x[(x >= 30) * (x<= 50)].size >= N_gt30_lt50 )* + (x.size >= N_points_tot) * \ + ((fit_type == 'exp1') + \ + ((fit_type == 'poly2') * \ + (x[x <= 40].size >= min_N_ang) * \ + (x[x >= 40].size >= min_N_ang)))): + #quadratic fit + if fit_type == 'poly2': + #Quadratic fit, no constraints. + p_w=np.polyfit(x,Tb_data,2, w=weights) + Tb_data_fit=np.polyval(p_w,inc_ang) + #available at either side of the interpolation + count_angle=np.cumsum(1 + 0.*good) + good_4=np.nonzero(np.logical_and(count_angle > min_N_ang, (count_angle + min_N_ang) <= count_angle[-1]))[0] + only_good=good[good_4] + else: + sys.exit('fitfunction not available') + #==> replace Radiometric Accuracy w/ sqrt(average(Ra)^2) + RA_data=data[Tb_ind[v]+6,i,:] + RA_data[RA_data < 0]=np.nan + RA_data=RA_data[only_good] + RA_data_new=np.sqrt(np.mean(RA_data ** 2)) + stdv_data_new=np.sqrt(np.mean((Tb_data[good_4] - Tb_data_fit[only_good]) ** 2)) + #counts and stdv (equal for all angles) + data_out[Tb_ind[v],i,:]=np.nan + data_out[Tb_ind[v],i,only_good]=Tb_data_fit[only_good] + data_out[Tb_ind[v] + 2,i,:]=stdv_data_new + data_out[Tb_ind[v] + 4,i,:]=only_good.size + data_out[Tb_ind[v] + 6,i,:]=RA_data_new + else: + data_out[Tb_ind[v],i,:]=np.nan + data_out[Tb_ind[v] + 2,i,:]=np.nan + data_out[Tb_ind[v] + 4,i,:]=np.nan + data_out[Tb_ind[v] + 6,i,:]=np.nan + else: + sys.exit('Revise fitting options') + + data=data_out + if (~np.isnan(data)).any(): + data[np.isnan(data)]=no_data_value + write_bin_SMOS_reg(out_filename,col_ind,row_ind,inc_ang,data,\ + asc_flag,version,prep_version,o_start_time,\ + o_end_time,overwrite,N_out_fields,\ + write_ind_latlon,'SCLF1C') + + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/__init__.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/__init__.py new file mode 100644 index 0000000..8338639 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/__init__.py @@ -0,0 +1,2 @@ +from .preprocess_nc import preprocess_nc +from .SCLF1C_reg2fit import SCLF1C_reg2fit diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/__init__.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/__init__.py new file mode 100644 index 0000000..bdeb3c5 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/__init__.py @@ -0,0 +1,3 @@ +from .xy2hv import xy2hv +from .latlontime_to_celestial import latlontime_to_celestial +from .gal_atm_correction import gal_atm_correction diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/gal_atm_correction.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/gal_atm_correction.py new file mode 100644 index 0000000..fca4c10 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/gal_atm_correction.py @@ -0,0 +1,232 @@ +import os +import numpy as np +import xml.etree.ElementTree as ET +from src.readwrite import read_aux_gal_sm +from src.helper import latlontime_to_celestial + + +#------------------------------------------------------------------------- +# Remove galactic and atmospheric contributions. +# Correct SMOS Tb in order to mimic SMAP. +#------------------------------------------------------------------------- + +# INPUT: Tb_H, Tb_V Tb [K] +# lat, lon latitude, longitude [degree] +# angle, phi incidence, azimuth angle [degree] +# time time [Y M D HH MM SS] +# T_air surface air temperature [oC] +# T_surf surface temperature [oC] +# P_surf surface pressure [mbar] +# V_surf surface water vapour density [g/m^3] +# C efficiency [-] + +#========================================================================= +def sub2ind(array_shape, rows, cols): + ind = rows*array_shape[1] + cols + ind[ind < 0] = -1 + ind[ind >= array_shape[0]*array_shape[1]] = -1 + return ind + +def ind2sub(array_shape, ind): + ind[ind < 0] = -1 + ind[ind >= array_shape[0]*array_shape[1]] = -1 + rows = (ind.astype('int') / array_shape[1]) + cols = ind % array_shape[1] + return (rows, cols) + +def gal_atm_correction(Tb_H=None,Tb_V=None,lat=None,lon=None,angle=None, \ + phi=None,time=None,SMOS_AUX_path=None,T_air=None,T_surf=None, \ + P_surf=None,V_surf=None,C=None,*args,**kwargs): + + d2r=np.pi / 180 + C_2_K=273.15 + +# 0) reduce to good data: +#---------------------------- + + length_org=Tb_H.size + limited_ind = np.nonzero( \ + np.logical_not(np.isnan(lat)) & \ + np.logical_not(np.isnan(lon)) & \ + np.logical_not(np.isnan(angle)) & \ + np.logical_not(np.isnan(phi)) & \ + np.logical_not(np.isnan(T_surf)) & \ + np.logical_not(np.isnan(P_surf)) & \ + np.logical_not(np.isnan(V_surf)))[0] + + #print('galactic and atmospheric correction on ', str(len(limited_ind)),' of ',length_org) + Tb_ap_H_out = np.full([length_org],np.nan) + Tb_ap_V_out = np.full([length_org],np.nan) + Tb_BOA_H_out = np.full([length_org],np.nan) + Tb_BOA_V_out = np.full([length_org],np.nan) + + Tb_H=Tb_H[limited_ind] + Tb_V=Tb_V[limited_ind] + lat=lat[limited_ind] + lon=lon[limited_ind] + angle=angle[limited_ind] + phi=phi[limited_ind] + T_air=T_air[limited_ind] + T_surf=T_surf[limited_ind] + P_surf=P_surf[limited_ind] + V_surf=V_surf[limited_ind] + + if Tb_H.size > 0: + # 1) Remove galactic terms: + + # Get Tc+Tgal from SMOS AUX data + + # v724 + current_dir = os.path.dirname(os.path.abspath(__file__)) + data_file=SMOS_AUX_path + + sidx = data_file.index('SM_OPER_AUX_GAL_SM_') + p=data_file+'/'+data_file[sidx:sidx+60]+'.HDR' + doc = ET.parse(p) + root = doc.getroot() + + for tmp in root.iter("{http://213.170.46.150/smos/schemas}Min_RA"): + Min_ra = float(tmp.text) + del tmp + + for tmp in root.iter("{http://213.170.46.150/smos/schemas}Max_RA"): + Max_ra = float(tmp.text) + del tmp + + for tmp in root.iter("{http://213.170.46.150/smos/schemas}Min_DEC"): + Min_dec = float(tmp.text) + del tmp + + for tmp in root.iter("{http://213.170.46.150/smos/schemas}Max_DEC"): + Max_dec = float(tmp.text) + del tmp + + for tmp in root.iter("{http://213.170.46.150/smos/schemas}DELTA_RA"): + D_ra = float(tmp.text) + del tmp + + for tmp in root.iter("{http://213.170.46.150/smos/schemas}DELTA_DEC"): + D_dec = float(tmp.text) + del tmp + + TB_Sky_H,TB_Sky_V=read_aux_gal_sm(data_file+'/'+data_file[sidx:sidx+60]+'.DBL') + +#First row corresponds to largest declination. +#We need to reverse the row order from the smallest to the largest +#declination. This was *not* done in the orginal code used for the GRSL +#paper [bug in paper]. + + TB_Sky_H=TB_Sky_H[::-1,:] + TB_Sky_V=TB_Sky_V[::-1,:] + + alpha,delta,Th_G0,UTC=latlontime_to_celestial(lat,lon,angle,phi,time) + + idx = np.ravel_multi_index( \ + [np.round((delta-Min_dec)/D_dec).astype('int'), \ + np.round((alpha-Min_ra)/D_ra).astype('int')], \ + TB_Sky_H.shape, mode='raise', order='F') + + TcTgal_H=TB_Sky_H.flatten('F')[idx] + TcTgal_V=TB_Sky_V.flatten('F')[idx] + + #print('min-max Tsky (H) = ' + str(TcTgal_H.min()) + ' ' + str(TcTgal_H.max())) + #print('min-max Tsky (V) = ' + str(TcTgal_V.min()) + ' ' + str(TcTgal_V.max())) + + # Estimate epsilon + eps_H = C * Tb_H / (T_surf + C_2_K) + eps_V = C * Tb_V / (T_surf + C_2_K) + + #print('reset eps (H) = ' + str(np.nansum(eps_H > 1)) + ' out of ' + str(eps_H.size)) + #print('reset eps (V) = ' + str(np.nansum(eps_V > 1)) + ' out of ' + str(eps_V.size)) + + eps_H[eps_H > 1]=1 + eps_V[eps_V > 1]=1 + + #print('min-max eps (H) = ' + str(eps_H.min()) + ' ' + str(eps_H.max())) + #print('min-max eps (V) = ' + str(eps_V.min()) + ' ' + str(eps_V.max())) + + # Estimate atmospheric loss factor: + # coefficients according to L1B ATBD (slightly different from Peng et al., 2013) + L=1.00938 - 2.9626*10**(-5)*T_air + 1.6521*10**(-5)*(P_surf-900) + 1.0712*10**(-5)*V_surf + + #print('min-max T_air = ' + str(T_air.min()) + ' ' + str(T_air.max())) + #print('min-max T_surf = ' + str(T_surf.min())+ ' ' + str(T_surf.max())) + #print('min-max V_surf = ' + str(V_surf.min())+ ' ' + str(V_surf.max())) + + ang_corr=np.cos(40.0*d2r) * 1/np.cos(angle*d2r) + L = L**ang_corr + + #print('min-max angle = ' +str(angle.min()) + ' ' + str(angle.max())) + #print('min-max atm loss L = ' + str(L.min()) + ' ' +str(L.max())) + + # Remove galactic and cosmic radiation + H_gal_corr = TcTgal_H * (1 - eps_H) / (L*L) + V_gal_corr = TcTgal_V * (1 - eps_V) / (L*L) + + Tb_ap_H = Tb_H - H_gal_corr + Tb_ap_V = Tb_V - V_gal_corr + + #print('min-max gal corr (H) = -(' + str(H_gal_corr.min()) + ' ' + str(H_gal_corr.max())+')') + #print('min-max gal corr (V) = -(' + str(V_gal_corr.min()) + ' ' + str(V_gal_corr.max())+')') + + # 2) Remove atmospheric effects: + + # Estimate upwelling Tb: + # T_up, temporary, coefficients according to L1B ATBD (slightly different from Peng et al., 2013) + T_up=2.3058 - 3.2735*10**(-3)*T_air + 4.233*10**(-3)*(P_surf - 900) + 1.4472*10**(-3)*V_surf + + # angular function to correct T_up, coefficients according to Peng et al., 2013 + f_ang=np.nan*angle + + ind = np.nonzero(angle < 20)[0] + + if ind.size: + f_ang[ind] = 1.2855*10**(-4)*angle[ind]**2 - 1.3361*10**(-4)*angle[ind] + 0.7625 + + ind=np.nonzero(np.logical_and( angle >= 20, angle <= 60))[0] + + if ind.size: + f_ang[ind] = 8.2724*10**(-6)*angle[ind]**3 - 5.7129*10**(-4)*angle[ind]**2 + 2.0411*10**(-2)*angle[ind] + 0.5655 + + ind=np.nonzero(np.logical_and(angle > 60,angle <= 70))[0] + + if ind.size: + f_ang[ind]=2.4189*10**(-3)*angle[ind]**2 - 0.2458*angle[ind]+ 7.5624 + + # get final T_up + T_up = T_up*f_ang + + #print('min-max T_up = ' + str(T_up.min()) + ' ' + str(T_up.max())) + + # Remove atmospheric correction, using Eq. below Eq. 5.55 in L1B ATBD + Tb_BOA_H = (T_surf + C_2_K)*(Tb_ap_H*L - (1 + L)*T_up) / (T_surf + C_2_K - T_up) + Tb_BOA_V = (T_surf + C_2_K)*(Tb_ap_V*L - (1 + L)*T_up) / (T_surf + C_2_K - T_up) + + H_atm_corr=Tb_BOA_H - Tb_ap_H + V_atm_corr=Tb_BOA_V - Tb_ap_V + + # prevent an increase in Tb + # an increase would happen occasionally, + # especically for V-pol and higher incidence angles, + # when T_up becomes larger and Tb_BOA_V > Tb_ap_V + #print('reset atm corr (H) = ' + str(np.nansum(H_atm_corr > 0)) + ' out of ' + str(H_atm_corr.size)) + #print('reset atm corr (V) = ' + str(np.nansum(V_atm_corr > 0)) + ' out of ' + str(V_atm_corr.size)) + + H_atm_corr[H_atm_corr > 0]=0 + V_atm_corr[V_atm_corr > 0]=0 + + Tb_BOA_H=Tb_ap_H + H_atm_corr + Tb_BOA_V=Tb_ap_V + V_atm_corr + + #print('min-max atm corr (H) = ' + str(H_atm_corr.min()) + ' ' + str(H_atm_corr.max())) + #print('min-max atm corr (V) = ' + str(V_atm_corr.min()) + ' ' + str(V_atm_corr.max())) + # 3) Expand to original length: +#-------------------------- + Tb_ap_H_out[limited_ind]=Tb_ap_H + Tb_ap_V_out[limited_ind]=Tb_ap_V + + Tb_BOA_H_out[limited_ind]=Tb_BOA_H + Tb_BOA_V_out[limited_ind]=Tb_BOA_V + + return Tb_ap_H_out, Tb_ap_V_out, Tb_BOA_H_out, Tb_BOA_V_out + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/latlontime_to_celestial.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/latlontime_to_celestial.py new file mode 100644 index 0000000..fe304f6 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/latlontime_to_celestial.py @@ -0,0 +1,92 @@ +import numpy as np +from datetime import datetime + +def latlontime_to_celestial(lat=None,lon=None,theta=None,phi=None,time=None,*args,**kwargs): + +# INPUT: +# lat = latitude [degree] +# lon = longitude [degree] +# theta = incidence angle [degree] +# phi = azimuth [degree] +# time = vector: date + time [Y M D HH MM SS] + +# OUTPUT: +# alpha = right ascension [degree] +# delta = declination [degree] +#--------------------------------------------------------------------- + + rad2deg = 180/np.pi + theta = theta/rad2deg + phi = phi/rad2deg + lat = lat/rad2deg + lon = lon/rad2deg + +# find elapsed minutes since Jan 1, 2000 (reference epoch) +#------------------------------------------- + + t1 = datetime(2000,1,1,0,0,0) + t2 = time + +# time diff between t1 and t2 + + UTC = t2-t1 + UTC = UTC.total_seconds()/86400 #in day + +# convert (lat,lon,az,inc,time) to celestial coordinates +#------------------------------------------- +#(theta, pi+phi) = specular direction points to the location in the celestial sky +#standard astronomical reference local geographic coordinates, +#elevation el, and azimuth Az + + el = np.pi/2 - theta #elevation or: 90-theta_in_degrees + Az = phi #astronom azimuth(0 towards south?) + +#celestial coordinates = declination (delta), right ascension (alpha) + +#--- alpha (RA) +# The right ascension is the angle between an object and the location +# of the vernal equinox (First Point in Aries) measured eastward along +# the celestial equator in hours, minutes, and seconds of sidereal time. +# Since the location of the vernal equinox changes due to the precession +# of the Earth's axis of rotation, coordinates must be given with +# reference to a date or epoch. + + JD2000_DAY = 2451544.5 #reference epoch, 1 Jan 2000 + JD_CENT = 36525 #Julian century + Omega_E = 0.2506846*60.0*24.0 #deg/day; Earth rotation rate; (15 degree/h) + + Y = time.year + M = time.month + D = time.day + HH = time.hour + MM = time.minute + SS = time.second + + C0 = 1721013.5 + D + (HH + MM/60 + SS/3600)/24 + JD = 367*Y - np.floor(1.75*(Y + np.floor((M + 9)/12))) + np.floor(275*M/9) + C0 #julian date + + if np.abs(JD - (UTC + JD2000_DAY)) > 1: + print('is the JD calculated correctly?') + + U0 = (JD - JD2000_DAY) / JD_CENT + Th_G0 = 100.46062 + 36000.77*U0 + 0.000388*U0**2 - 2.6*10**(-8)*U0**3 # degree + Th_G0 = np.mod(Th_G0,360) + + Th_L = Th_G0 + Omega_E*np.mod(UTC,1) + lon*rad2deg # degree; local sideral time + Th_L = np.mod(Th_L,360) + + tmp_H = np.sin(Az) / (np.tan(el)*np.cos(lat) + np.cos(Az)*np.sin(lat)) + H = np.arctan(tmp_H)*rad2deg + + alpha = np.mod(Th_L-H, 360) + +#--- delta (DEC) +# The declination of an object is its angle in degrees, minutes, and seconds of arc +# above or below the celestial equator. + + tmp_delta = np.sin(lat)*np.sin(el) - np.cos(lat)*np.cos(el)*np.cos(Az) + delta = np.arcsin(tmp_delta) * rad2deg + + return alpha, delta, Th_G0, UTC + + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/util.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/util.py new file mode 100644 index 0000000..944f2c8 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/util.py @@ -0,0 +1,56 @@ +import os +import glob +import argparse +import subprocess +import multiprocessing +from datetime import datetime, timedelta + +def parse_args(): + """ + No arguments: process yesterday only (the normal NRT cron behavior). + --date: backfill a single specific day. + --start [--end]: backfill an inclusive range of days. + """ + parser = argparse.ArgumentParser( + description="Process SMOS L1C data. With no arguments, processes " + "yesterday's data (NRT mode). Use --date or --start/--end " + "to backfill missing days." + ) + group = parser.add_mutually_exclusive_group() + group.add_argument('--date', metavar='YYYYMMDD', + help="Backfill a single specific day.") + group.add_argument('--start', metavar='YYYYMMDD', + help="First day of a backfill range (inclusive).") + parser.add_argument('--end', metavar='YYYYMMDD', + help="Last day of a backfill range (inclusive). " + "Only valid together with --start; defaults to --start.") + args = parser.parse_args() + + if args.end and not args.start: + parser.error("--end requires --start") + if args.date and args.end: + parser.error("--date cannot be combined with --end") + + return args + + +def get_time_range(args): + """ + Returns (start, end) datetimes, where `end` is exclusive (the day after + the last day to process). + """ + yesterday = datetime.today().replace(hour=0, minute=0, second=0, microsecond=0) - timedelta(days=1) + + if args.date: + start = datetime.strptime(args.date, '%Y%m%d') + end_day = start + elif args.start: + start = datetime.strptime(args.start, '%Y%m%d') + end_day = datetime.strptime(args.end, '%Y%m%d') if args.end else start + else: + start = yesterday + end_day = yesterday + + return start, end_day + timedelta(days=1) + + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/xy2hv.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/xy2hv.py new file mode 100644 index 0000000..6cd0aeb --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/helper/xy2hv.py @@ -0,0 +1,512 @@ +import numpy as np +from numpy import matlib + +def invMR4L(alpha = None,*args,**kwargs): + + cosdalpha = np.cos(np.deg2rad(alpha)) + sindalpha = np.sin(np.deg2rad(alpha)) + cosd2alpha = np.cos(np.deg2rad(2*alpha)) + sind2alpha = np.sin(np.deg2rad(2*alpha)) + + MR4 = np.full([len(alpha),4,4],np.nan) + + MR4[:,0,0] = cosdalpha ** 2 + MR4[:,0,1] = sindalpha ** 2 + MR4[:,0,2] = -cosdalpha * sindalpha + MR4[:,0,3] = 0 + MR4[:,1,0] = sindalpha ** 2 + MR4[:,1,1] = cosdalpha ** 2 + MR4[:,1,2] = cosdalpha*sindalpha + MR4[:,1,3] = 0 + MR4[:,2,0] = sind2alpha + MR4[:,2,1] = -sind2alpha + MR4[:,2,2] = cosd2alpha + MR4[:,2,3] = 0 + MR4[:,3,0] = 0 + MR4[:,3,1] = 0 + MR4[:,3,2] = 0 + MR4[:,3,3] = 1 + + IMR4 = np.full(MR4.shape,np.nan) + + IMR4[:,0,0]=MR4[:,1,2]*MR4[:,2,3]*MR4[:,3,1] - MR4[:,1,3]*MR4[:,2,2]*MR4[:,3,1] + MR4[:,1,3]*MR4[:,2,1]*MR4[:,3,2] - MR4[:,1,1]*MR4[:,2,3]*MR4[:,3,2] - MR4[:,1,2]*MR4[:,2,1]*MR4[:,3,3] + MR4[:,1,1]*MR4[:,2,2]*MR4[:,3,3] + IMR4[:,0,1]=MR4[:,0,3]*MR4[:,2,2]*MR4[:,3,1] - MR4[:,0,2]*MR4[:,2,3]*MR4[:,3,1] - MR4[:,0,3]*MR4[:,2,1]*MR4[:,3,2] + MR4[:,0,1]*MR4[:,2,3]*MR4[:,3,2] + MR4[:,0,2]*MR4[:,2,1]*MR4[:,3,3] - MR4[:,0,1]*MR4[:,2,2]*MR4[:,3,3] + IMR4[:,0,2]=MR4[:,0,2]*MR4[:,1,3]*MR4[:,3,1] - MR4[:,0,3]*MR4[:,1,2]*MR4[:,3,1] + MR4[:,0,3]*MR4[:,1,1]*MR4[:,3,2] - MR4[:,0,1]*MR4[:,1,3]*MR4[:,3,2] - MR4[:,0,2]*MR4[:,1,1]*MR4[:,3,3] + MR4[:,0,1]*MR4[:,1,2]*MR4[:,3,3] + IMR4[:,0,3]=MR4[:,0,3]*MR4[:,1,2]*MR4[:,2,1] - MR4[:,0,2]*MR4[:,1,3]*MR4[:,2,1] - MR4[:,0,3]*MR4[:,1,1]*MR4[:,2,2] + MR4[:,0,1]*MR4[:,1,3]*MR4[:,2,2] + MR4[:,0,2]*MR4[:,1,1]*MR4[:,2,3] - MR4[:,0,1]*MR4[:,1,2]*MR4[:,2,3] + IMR4[:,1,0]=MR4[:,1,3]*MR4[:,2,2]*MR4[:,3,0] - MR4[:,1,2]*MR4[:,2,3]*MR4[:,3,0] - MR4[:,1,3]*MR4[:,2,0]*MR4[:,3,2] + MR4[:,1,0]*MR4[:,2,3]*MR4[:,3,2] + MR4[:,1,2]*MR4[:,2,0]*MR4[:,3,3] - MR4[:,1,0]*MR4[:,2,2]*MR4[:,3,3] + IMR4[:,1,1]=MR4[:,0,2]*MR4[:,2,3]*MR4[:,3,0] - MR4[:,0,3]*MR4[:,2,2]*MR4[:,3,0] + MR4[:,0,3]*MR4[:,2,0]*MR4[:,3,2] - MR4[:,0,0]*MR4[:,2,3]*MR4[:,3,2] - MR4[:,0,2]*MR4[:,2,0]*MR4[:,3,3] + MR4[:,0,0]*MR4[:,2,2]*MR4[:,3,3] + IMR4[:,1,2]=MR4[:,0,3]*MR4[:,1,2]*MR4[:,3,0] - MR4[:,0,2]*MR4[:,1,3]*MR4[:,3,0] - MR4[:,0,3]*MR4[:,1,0]*MR4[:,3,2] + MR4[:,0,0]*MR4[:,1,3]*MR4[:,3,2] + MR4[:,0,2]*MR4[:,1,0]*MR4[:,3,3] - MR4[:,0,0]*MR4[:,1,2]*MR4[:,3,3] + IMR4[:,1,3]=MR4[:,0,2]*MR4[:,1,3]*MR4[:,2,0] - MR4[:,0,3]*MR4[:,1,2]*MR4[:,2,0] + MR4[:,0,3]*MR4[:,1,0]*MR4[:,2,2] - MR4[:,0,0]*MR4[:,1,3]*MR4[:,2,2] - MR4[:,0,2]*MR4[:,1,0]*MR4[:,2,3] + MR4[:,0,0]*MR4[:,1,2]*MR4[:,2,3] + IMR4[:,2,0]=MR4[:,1,1]*MR4[:,2,3]*MR4[:,3,0] - MR4[:,1,3]*MR4[:,2,1]*MR4[:,3,0] + MR4[:,1,3]*MR4[:,2,0]*MR4[:,3,1] - MR4[:,1,0]*MR4[:,2,3]*MR4[:,3,1] - MR4[:,1,1]*MR4[:,2,0]*MR4[:,3,3] + MR4[:,1,0]*MR4[:,2,1]*MR4[:,3,3] + IMR4[:,2,1]=MR4[:,0,3]*MR4[:,2,1]*MR4[:,3,0] - MR4[:,0,1]*MR4[:,2,3]*MR4[:,3,0] - MR4[:,0,3]*MR4[:,2,0]*MR4[:,3,1] + MR4[:,0,0]*MR4[:,2,3]*MR4[:,3,1] + MR4[:,0,1]*MR4[:,2,0]*MR4[:,3,3] - MR4[:,0,0]*MR4[:,2,1]*MR4[:,3,3] + IMR4[:,2,2]=MR4[:,0,1]*MR4[:,1,3]*MR4[:,3,0] - MR4[:,0,3]*MR4[:,1,1]*MR4[:,3,0] + MR4[:,0,3]*MR4[:,1,0]*MR4[:,3,1] - MR4[:,0,0]*MR4[:,1,3]*MR4[:,3,1] - MR4[:,0,1]*MR4[:,1,0]*MR4[:,3,3] + MR4[:,0,0]*MR4[:,1,1]*MR4[:,3,3] + IMR4[:,2,3]=MR4[:,0,3]*MR4[:,1,1]*MR4[:,2,0] - MR4[:,0,1]*MR4[:,1,3]*MR4[:,2,0] - MR4[:,0,3]*MR4[:,1,0]*MR4[:,2,1] + MR4[:,0,0]*MR4[:,1,3]*MR4[:,2,1] + MR4[:,0,1]*MR4[:,1,0]*MR4[:,2,3] - MR4[:,0,0]*MR4[:,1,1]*MR4[:,2,3] + IMR4[:,3,0]=MR4[:,1,2]*MR4[:,2,1]*MR4[:,3,0] - MR4[:,1,1]*MR4[:,2,2]*MR4[:,3,0] - MR4[:,1,2]*MR4[:,2,0]*MR4[:,3,1] + MR4[:,1,0]*MR4[:,2,2]*MR4[:,3,1] + MR4[:,1,1]*MR4[:,2,0]*MR4[:,3,2] - MR4[:,1,0]*MR4[:,2,1]*MR4[:,3,2] + IMR4[:,3,1]=MR4[:,0,1]*MR4[:,2,2]*MR4[:,3,0] - MR4[:,0,2]*MR4[:,2,1]*MR4[:,3,0] + MR4[:,0,2]*MR4[:,2,0]*MR4[:,3,1] - MR4[:,0,0]*MR4[:,2,2]*MR4[:,3,1] - MR4[:,0,1]*MR4[:,2,0]*MR4[:,3,2] + MR4[:,0,0]*MR4[:,2,1]*MR4[:,3,2] + IMR4[:,3,2]=MR4[:,0,2]*MR4[:,1,1]*MR4[:,3,0] - MR4[:,0,1]*MR4[:,1,2]*MR4[:,3,0] - MR4[:,0,2]*MR4[:,1,0]*MR4[:,3,1] + MR4[:,0,0]*MR4[:,1,2]*MR4[:,3,1] + MR4[:,0,1]*MR4[:,1,0]*MR4[:,3,2] - MR4[:,0,0]*MR4[:,1,1]*MR4[:,3,2] + IMR4[:,3,3]=MR4[:,0,1]*MR4[:,1,2]*MR4[:,2,0] - MR4[:,0,2]*MR4[:,1,1]*MR4[:,2,0] + MR4[:,0,2]*MR4[:,1,0]*MR4[:,2,1] - MR4[:,0,0]*MR4[:,1,2]*MR4[:,2,1] - MR4[:,0,1]*MR4[:,1,0]*MR4[:,2,2] + MR4[:,0,0]*MR4[:,1,1]*MR4[:,2,2] + + return IMR4 + +def xy2hv(a_BT_Real=None,a_BT_Imag=None,a_RA=None,a_Theta=None,\ + a_Az=None,a_alpha=None, a_Snap_ID=None,a_t_smos_sec=None,\ + lat=None,lon=None,a_flag_15=None,a_flag_16=None,\ + mask_ok=None,dgg_list=None,BT_count=None,\ + in_dgg=None, *args, **kwargs): + ''' + # This function do the rotation of SMOS data from the antenna reference + # (XY) to Earth reference (HV) - Top Of Atmosphere. + # + # INPUTS : + # = output from call to TB nc file + + # OUTPUT : + # = specific variables in a long 1D array + + # Interpolation and rotation can only be done if the number of available + # angles is greater than 6 for full pol. If the roation + # cannot be performed, the strucutre will be filled with NaN and a message + # will be displayed on screen. + # + # Based on: + # ========== + # Authors : Delphine Leroux (delphine.leroux@cesbio.cnes.fr) + # Creation : 29/10/2010 (v1.1) + # Version : 1.8 + # Comments : - this fnunction has been tested on 64-bit Linux with RWAPI + # v1.3 and Matlab R2010a + # - in dual polarization, the user should be aware that the + # rotation matrix MR2 is not invertible around +-45° so the + # inversion might not succeed or will give very high values of + # brightness temperature. Therefore, the user can set a + # threshold or can filter the results. + # v1.1 : no need to put the polarization as an input (directly read in the + # filename) - 29/10/2010 + # v1.2 : in the full pol part, the beginning of the interpolation has been + # changed (only the time field gives the order to sort all the other + # fields) & lat, lon, lat and mask are stored - 02/11/2010 + # v1.3 : a reciprocal condition number threshold has been added such that + # if the rcond number of MR2 is less than this threshold, the MR2 matrix is + # not inverted - only for DUAL pol - 02/11/2010 + # v1.4 : The flag for each measurement is now stored and can be used to + # look at the RFI flag & a test has been added to check if the DGG of + # interest are located in the input product - 04/11/2010 + # v1.5 : Bug fixed (lat/long were not stored correctly) & problem with the + # alpha angle used to do the rotation fixed - 08/11/2010 + # v1.6 : minor changes to gain some speed - 10/11/2010 + # V1.7 : convention of the rotation angle changing depending on the version + # of L1OP. For v344 or later a=Fa+Ge and before v344 a=Fa-Ge - 15/11/2010 + # v1.8 : The azimuth angles are now stored and the flag are stored + # differently (now it is a matrix) - 17/11/2010 + # CESBIO - 2010 + + # UPDATED: + #========== + # Gabrielle De Lannoy, NASA/GMAO, 17dec10: pre-allocated memory, + # removed loops, + # updated for F pol only + # Gabrielle De Lannoy, NASA/GMAO, 08mar11: take mex-output, + # insert optimization (FC) + # Gabrielle De Lannoy, NASA/GMAO, 24mar11: reduce data coming in from reader + # by a priori selection on flags, + # limit sort/find/loop - commands, + # correct RA (to match v1.10) + # Gabrielle De Lannoy, NASA/GMAO, 02dec14: added azimuth output + # Gabrielle De Lannoy, NASA/GMAO, 24aug15: added check on interpolation + #===================================================== + ''' + #nargin = XY2HV_F_mex_opt_2dec14.nargin + + # 10 K = 2*4K +2K + max_dev=10 + + # don't smooth when there is a gap of more than 1 degree + theta_threshold=1 + + #if nargin < 19 : + in_dgg = dgg_list + ia = np.arange(len(in_dgg)) + + if in_dgg.size >= 1: + ####################### + ## FULL POLARIZATION ## + ####################### + # Initializing the output structure + + # number of *valid* (no RFI, ALIAS-free) snapshots/inc angles per grid cell + max_el = 250 + max_tot_el = len(in_dgg) * max_el + + out_lon= np.full([max_tot_el], np.nan) + out_lat= np.full([max_tot_el], np.nan) + out_inc= np.full([max_tot_el], np.nan) + out_Az= np.full([max_tot_el], np.nan) + out_Tbh= np.full([max_tot_el], np.nan) + out_Tbv= np.full([max_tot_el], np.nan) + out_T3= np.full([max_tot_el], np.nan) + out_T4= np.full([max_tot_el], np.nan) + out_RAh= np.full([max_tot_el], np.nan) + out_RAv= np.full([max_tot_el], np.nan) + out_RA3= np.full([max_tot_el], np.nan) + out_RA4= np.full([max_tot_el], np.nan) + + out_Xorg= np.full([max_tot_el], 0.) + out_Yorg= np.full([max_tot_el], 0.) + out_RXYorg= np.full([max_tot_el], 0.) + out_IXYorg= np.full([max_tot_el], 0.) + + end_id_list = np.full([len(in_dgg)], np.nan) + + start_id = 0 + + for idx_in_dgg in np.arange(len(in_dgg)): + + idx_dgg = ia[idx_in_dgg] + TB_count = BT_count[idx_dgg] + + # only select land points, + # and limit data to those points at a distance more than 40km away + # from 'coastline' (i.e. water) + # => taken together in mask_ok (alias and RFI flags are + # accounted for in reader) + + if (TB_count > 6 and mask_ok[idx_dgg] and lat[idx_dgg] > -60): + + flag_15 = a_flag_15[idx_dgg, :TB_count] + flag_16 = a_flag_16[idx_dgg, :TB_count] + BT_Real = a_BT_Real[idx_dgg, :TB_count] + BT_Imag = a_BT_Imag[idx_dgg, :TB_count] + + RA = a_RA[idx_dgg,:TB_count] + Theta = a_Theta[idx_dgg,:TB_count] + Az = a_Az[idx_dgg,:TB_count] + alpha = a_alpha[idx_dgg,:TB_count] + Snap_ID = a_Snap_ID[idx_dgg,:TB_count] + t_smos_sec = a_t_smos_sec[idx_dgg,:TB_count] + + # Only consider unique incidence angle, unique acquisition time and unique + # snapshot ID even for mixed snapshot (XX and XY or YY and XY) + # Unique and in the same order as found (sort(index)) + dum,b_tsmos = np.unique(t_smos_sec,return_index=True) + + b_tsmos = b_tsmos[np.logical_not(np.isnan(dum))] + ind = np.sort(b_tsmos) + t_smos_sec_uniq = t_smos_sec[ind] + Theta_uniq = Theta[ind] + Az_uniq = Az[ind] + alpha_uniq = alpha[ind] + Snap_ID_uniq = Snap_ID[ind] + + if Snap_ID_uniq.size == 0: + + end_id = start_id + 1 + end_id_list[idx_in_dgg] = end_id + out_lon[start_id:end_id] = np.nan + out_lat[start_id:end_id] = np.nan + out_Tbh[start_id:end_id] = np.nan + out_Tbv[start_id:end_id] = np.nan + out_T3[start_id:end_id] = np.nan + out_T4[start_id:end_id] = np.nan + out_RAh[start_id:end_id] = np.nan + out_RAv[start_id:end_id] = np.nan + out_RA3[start_id:end_id] = np.nan + out_RA4[start_id:end_id] = np.nan + + out_Xorg[start_id:end_id] = np.nan + out_Yorg[start_id:end_id] = np.nan + out_RXYorg[start_id:end_id] = np.nan + out_IXYorg[start_id:end_id] = np.nan + out_inc[start_id:end_id] = np.nan + out_Az[start_id:end_id] = np.nan + start_id = end_id + + else: + + # np.nonzero return a tuple, thus need to specify [0] for 1-D idx + idx_TBxx = np.nonzero(np.logical_and(flag_15==0, flag_16 == 0))[0] + idx_TByy = np.nonzero(np.logical_and(flag_15==0, flag_16 == 1))[0] + idx_TBxy = np.nonzero(np.logical_or(np.logical_and(flag_15==1,flag_16==0), \ + np.logical_and(flag_15==1,flag_16==1)))[0] + + # TBxy and Flags_final = output; since we filter out for + # the flags in the step above, there is no need to have flags in + # output + + # TBxy : 1-TBxx 2-TByy 3-Re(TBxy) 4-Im(TBxy) 5-Theta 6-SnapID 7-t_smos_sec + # 8-RA_TBxx 9-RA_TByy 10-RA_TBxy + + TBxy = np.full([len(Snap_ID_uniq),10],np.nan) + + TBxy[:,4] = Theta_uniq + TBxy[:,5] = Snap_ID_uniq + TBxy[:,6] = t_smos_sec_uniq + + for i in np.arange(len(idx_TBxx)): + idx_snap = np.nonzero(Snap_ID_uniq == Snap_ID[idx_TBxx[i]])[0] + TBxy[idx_snap,0] = BT_Real[idx_TBxx[i]] + TBxy[idx_snap,7] = RA[idx_TBxx[i]] + for i in np.arange(len(idx_TByy)): + idx_snap = np.nonzero(Snap_ID_uniq == Snap_ID[idx_TByy[i]])[0] + TBxy[idx_snap,1] = BT_Real[idx_TByy[i]] + TBxy[idx_snap,8] = RA[idx_TByy[i]] + for i in np.arange(len(idx_TBxy)): + idx_snap = np.nonzero(Snap_ID_uniq == Snap_ID[idx_TBxy[i]])[0] + TBxy[idx_snap,2] = BT_Real[idx_TBxy[i]] + TBxy[idx_snap,3] = BT_Imag[idx_TBxy[i]] + TBxy[idx_snap,9] = RA[idx_TBxy[i]] + + #------------------------------------------------------------------ + # Eliminate outlier data before interpolation + #------------------------------------------------------------------ + # For that we need to do moving average per angle, because the TB + # will not necessarily follow a smooth path over time, and there + # may be unforeseen angle-gaps when smoothing over snapshots... + # It would be harder to select sub-windows of 'similar' snapshots; + # now we smooth only over continuous sets of angles + # ==> this could be easily reverted by limiting the + # snapshot-windows to periods with minimal time differences... + + TBxy_thetasorted = TBxy[np.argsort(TBxy[:,4]),:] + theta_1 = TBxy_thetasorted[0:-1,4] + theta_2 = TBxy_thetasorted[1:,4] + smooth_end = np.nonzero(np.abs(theta_1-theta_2) > theta_threshold)[0] + + smooth_end = np.concatenate(([-1],smooth_end,[len(theta_1)+1]),axis=None) + + #Per continuous set of angles, do the smoothing and throw out outliers + N = 5 + + for i in np.arange(len(smooth_end)-1): + + # Nothing is done on outliers + # if there is no more than N consecutive angles + # if np.abs((smooth_end[i]+1) - smooth_end[i+1]) >= 5: + if np.abs((smooth_end[i]+1) - smooth_end[i+1]) > 5: + + tmp_TBxy = TBxy_thetasorted[smooth_end[i]+1:smooth_end[i+1],:4] + tmp = np.nancumsum(tmp_TBxy,0) + count_tmp = np.cumsum(np.logical_not(np.isnan(tmp_TBxy)),0) + + #Initial rows of tmp may contain 0, + #if there happened to be leading nan-values for some fields. + #You would hope to have at least 1 good value in the first + #window of N angles - if not, look for the first non-zero + #element and replace zero cumsum by this first value + #assign window-mean to the middle element + #if count_tmp(N+1:end,:) - count_tmp(1:end-N,:) happens to be + #zero, then NaN results in the window mean, which is fine, + #because there won't be any data to check in this window anyway + #In case there was no single good value in the first N angles: + #Replace leading zeros (~nan) with first non-zero element + + for jj in np.arange(4): + if (count_tmp[0,jj] == 0): + # index of first non-zero only + ind=np.nonzero(tmp[:,jj] != 0)[0] + if ind.size > 0: + tmp[np.arange(ind[0]),jj]=tmp[ind[0],jj] + count_tmp[np.arange(ind[0]),jj]=1 + + aux = tmp[N,:] / count_tmp[N,:] #for first window middle + tmp[int(np.ceil(N/2)):-int(np.floor(N/2)),:] = \ + (tmp[N:,:]-tmp[0:-N,:]) / (count_tmp[N:,:]-count_tmp[0:-N,:]) + tmp[int(np.ceil(N/2))-1,:]=aux + + # Assign constant (mean) values to tails of the moving window + tmp[0:int(np.ceil(N/2))-1,:]=np.matlib.repmat(tmp[int(np.ceil(N/2))-1,:], int(np.ceil(N/2))-1,1) + tmp[-int(np.floor(N/2)):,:]=np.matlib.repmat(tmp[-int(np.floor(N/2))-1,:],int(np.floor(N/2)),1) + + # Remove outliers + tmp_TBxy[np.abs(tmp_TBxy-tmp) > max_dev] = np.nan + TBxy_thetasorted[smooth_end[i]+1:smooth_end[i+1],0:4]=tmp_TBxy + + del tmp, count_tmp + + # sort back on time + TBxy = TBxy_thetasorted[np.argsort(TBxy_thetasorted[:,6]),:] + + #-------------------------------------------------------------------- + + # Locate the X, Y and XY polarizations + bool_Flags_X = np.logical_not(np.isnan(TBxy[:,0])) + bool_Flags_Y = np.logical_not(np.isnan(TBxy[:,1])) + bool_Flags_XY = np.logical_and(np.logical_not(np.isnan(TBxy[:,2])),\ + np.logical_not(np.isnan(TBxy[:,3]))) + + # To keep track of where original data were before interpolation + bool_Flags_RXY = np.logical_not(np.isnan(TBxy[:,2])) + bool_Flags_IXY = np.logical_not(np.isnan(TBxy[:,3])) + + # Identify if the possibleX, Y and XY polarizations are for the + # interpolation + bool_Flags_pre2_X = bool_Flags_X[0:-4] + bool_Flags_pre1_X = bool_Flags_X[1:-3] + bool_Flags_fol1_X = bool_Flags_X[3:-1] + bool_Flags_fol2_X = bool_Flags_X[4:] + + bool_Flags_pre2_Y = bool_Flags_Y[0:-4] + bool_Flags_pre1_Y = bool_Flags_Y[1:-3] + bool_Flags_fol1_Y = bool_Flags_Y[3:-1] + bool_Flags_fol2_Y = bool_Flags_Y[4:] + + bool_Flags_pre1_XY = bool_Flags_XY[1:-3] + bool_Flags_fol1_XY =bool_Flags_XY[3:-1] + + # dentify if the possible polarizations for interpolation are not too far + # considering time acquisition + bool_t_smos_pre2 = (t_smos_sec_uniq[2:-2]-t_smos_sec_uniq[0:-4]) <= 2.5 + bool_t_smos_pre1 = (t_smos_sec_uniq[2:-2]-t_smos_sec_uniq[1:-3]) <= 1.3 + bool_t_smos_fol1 = (t_smos_sec_uniq[3:-1]-t_smos_sec_uniq[2:-2]) <= 1.3 + bool_t_smos_fol2 = (t_smos_sec_uniq[4:] -t_smos_sec_uniq[2:-2]) <= 2.5 + + # Interpolation possibilities for each polarization X, Y and XY + bool_interp_X = np.concatenate(([1<0],[1<0],\ + np.logical_and(np.logical_or(np.logical_and(bool_Flags_pre2_X,bool_t_smos_pre2), \ + np.logical_and(bool_Flags_pre1_X,bool_t_smos_pre1)),\ + np.logical_or(np.logical_and(bool_Flags_fol1_X,bool_t_smos_fol1), \ + np.logical_and(bool_Flags_fol2_X,bool_t_smos_fol2))), \ + [1<0],[1<0]),axis=None) + bool_interp_Y = np.concatenate(([1<0],[1<0],\ + np.logical_and(np.logical_or(np.logical_and(bool_Flags_pre2_Y,bool_t_smos_pre2), \ + np.logical_and(bool_Flags_pre1_Y,bool_t_smos_pre1)), \ + np.logical_or(np.logical_and(bool_Flags_fol1_Y,bool_t_smos_fol1), \ + np.logical_and(bool_Flags_fol2_Y,bool_t_smos_fol2))), \ + [1<0],[1<0]),axis=None) + bool_interp_XY= np.concatenate(([1<0],[1<0], \ + np.logical_and(np.logical_and(bool_Flags_pre1_XY,bool_t_smos_pre1), \ + np.logical_and(bool_Flags_fol1_XY,bool_t_smos_fol1)), \ + [1<0],[1<0]),axis=None) + + if (np.sum(bool_Flags_X) >= 2 and np.sum(bool_Flags_Y) >= 2 and \ + np.sum(bool_Flags_XY) >= 2): + + # Interpolation of X, Y and XY where needed + I_TBx = np.interp(t_smos_sec_uniq[bool_interp_X], t_smos_sec_uniq[bool_Flags_X],TBxy[bool_Flags_X,0]) + I_TBy = np.interp(t_smos_sec_uniq[bool_interp_Y],t_smos_sec_uniq[bool_Flags_Y],TBxy[bool_Flags_Y,1]) + I_TBre = np.interp(t_smos_sec_uniq[bool_interp_XY],t_smos_sec_uniq[bool_Flags_XY],TBxy[bool_Flags_XY,2]) + I_TBim = np.interp(t_smos_sec_uniq[bool_interp_XY],t_smos_sec_uniq[bool_Flags_XY],TBxy[bool_Flags_XY,3]) + + # Interpolation of the radiometric accuracies + I_RATBx = np.interp(t_smos_sec_uniq[bool_interp_X],t_smos_sec_uniq[bool_Flags_X],TBxy[bool_Flags_X,7]) + I_RATBy = np.interp(t_smos_sec_uniq[bool_interp_Y],t_smos_sec_uniq[bool_Flags_Y],TBxy[bool_Flags_Y,8]) + I_RATBxy = np.interp(t_smos_sec_uniq[bool_interp_XY],t_smos_sec_uniq[bool_Flags_XY],TBxy[bool_Flags_XY,9]) + + # Saving the interpolated values + TBxy[bool_interp_X,0]=I_TBx + TBxy[bool_interp_Y,1]=I_TBy + TBxy[bool_interp_XY,2]=I_TBre + TBxy[bool_interp_XY,3]=I_TBim + TBxy[bool_interp_X,7]=I_RATBx + TBxy[bool_interp_Y,8]=I_RATBy + TBxy[bool_interp_XY,9]=I_RATBxy + + # ROTATION & ERROR PROPAGATION + # Transformation + TBhv = np.full([len(Snap_ID_uniq),4],np.nan) + RATBhv = np.full([len(Snap_ID_uniq),4],np.nan) + + # Parallelizing the inversion of MR4 for all angles/Snap_ID - FC 01/11 + IMR4 = invMR4L(alpha_uniq) + TBhv[:,0] = IMR4[:,0,0]*TBxy[:,0] + IMR4[:,0,1]*TBxy[:,1] + IMR4[:,0,2]*TBxy[:,2]*2 - IMR4[:,0,3]*TBxy[:,3]*2 + TBhv[:,1] = IMR4[:,1,0]*TBxy[:,0] + IMR4[:,1,1]*TBxy[:,1] + IMR4[:,1,2]*TBxy[:,2]*2 - IMR4[:,1,3]*TBxy[:,3]*2 + TBhv[:,2] = IMR4[:,2,0]*TBxy[:,0] + IMR4[:,2,1]*TBxy[:,1] + IMR4[:,2,2]*TBxy[:,2]*2 - IMR4[:,2,3]*TBxy[:,3]*2 + TBhv[:,3] = IMR4[:,3,0]*TBxy[:,0] + IMR4[:,3,1]*TBxy[:,1] + IMR4[:,3,2]*TBxy[:,2]*2 - IMR4[:,3,3]*TBxy[:,3]*2 + + RATBhv[:,0] = (IMR4[:,0,0]**2*TBxy[:,7]**2 + IMR4[:,0,1]**2*TBxy[:,8]**2 + 4*(IMR4[:,0,2]**2 + IMR4[:,0,3]**2)*TBxy[:,9]**2)**0.5 + RATBhv[:,1] = (IMR4[:,1,0]**2*TBxy[:,7]**2 + IMR4[:,1,1]**2*TBxy[:,8]**2 + 4*(IMR4[:,1,2]**2 + IMR4[:,1,3]**2)*TBxy[:,9]**2)**0.5 + RATBhv[:,2] = (IMR4[:,2,0]**2*TBxy[:,7]**2 + IMR4[:,2,1]**2*TBxy[:,8]**2 + 4*(IMR4[:,2,2]**2 + IMR4[:,2,3]**2)*TBxy[:,9]**2)**0.5 + RATBhv[:,3] = (IMR4[:,3,0]**2*TBxy[:,7]**2 + IMR4[:,3,1]**2*TBxy[:,8]**2 + 4*(IMR4[:,3,2]**2 + IMR4[:,3,3]**2)*TBxy[:,9]**2)**0.5 + + # OUTPUT STRUCTURE + end_id = start_id + len(Snap_ID_uniq) + end_id_list[idx_in_dgg] = end_id + out_lon[start_id:end_id]=np.matlib.repmat(lon[idx_dgg],1,len(Snap_ID_uniq)) + out_lat[start_id:end_id]=np.matlib.repmat(lat[idx_dgg],1,len(Snap_ID_uniq)) + out_Tbh[start_id:end_id]=TBhv[:,0] + out_Tbv[start_id:end_id]=TBhv[:,1] + out_T3[start_id:end_id]=TBhv[:,2] + out_T4[start_id:end_id]=TBhv[:,3] + out_RAh[start_id:end_id]=RATBhv[:,0] + out_RAv[start_id:end_id]=RATBhv[:,1] + out_RA3[start_id:end_id]=RATBhv[:,2] + out_RA4[start_id:end_id]=RATBhv[:,3] + + out_Xorg[start_id + np.nonzero(bool_Flags_X == 1)[0]] = 1 + out_Yorg[start_id + np.nonzero(bool_Flags_Y == 1)[0]] = 1 + out_RXYorg[start_id + np.nonzero(bool_Flags_RXY == 1)[0]] = 1 + out_IXYorg[start_id + np.nonzero(bool_Flags_IXY == 1)[0]] = 1 + + out_inc[start_id:end_id] = Theta_uniq + out_Az[start_id:end_id] = Az_uniq + start_id=end_id + + else: + + end_id = start_id + 1 + end_id_list[idx_in_dgg] = end_id + out_lon[start_id:end_id] = np.nan + out_lat[start_id:end_id] = np.nan + out_Tbh[start_id:end_id] = np.nan + out_Tbv[start_id:end_id] = np.nan + out_T3[start_id:end_id] = np.nan + out_T4[start_id:end_id] = np.nan + out_RAh[start_id:end_id] = np.nan + out_RAv[start_id:end_id] = np.nan + out_RA3[start_id:end_id] = np.nan + out_RA4[start_id:end_id] = np.nan + + out_Xorg[start_id:end_id] = np.nan + out_Yorg[start_id:end_id] = np.nan + out_RXYorg[start_id:end_id] = np.nan + out_IXYorg[start_id:end_id] = np.nan + out_inc[start_id:end_id] = np.nan + out_Az[start_id:end_id] = np.nan + start_id = end_id + + else: + + end_id = start_id +1 + end_id_list[idx_in_dgg] = end_id + out_lon[start_id:end_id] = np.nan + out_lat[start_id:end_id] = np.nan + out_Tbh[start_id:end_id] = np.nan + out_Tbv[start_id:end_id] = np.nan + out_T3[start_id:end_id] = np.nan + out_T4[start_id:end_id] = np.nan + out_RAh[start_id:end_id] = np.nan + out_RAv[start_id:end_id] = np.nan + out_RA3[start_id:end_id] = np.nan + out_RA4[start_id:end_id] = np.nan + + out_Xorg[start_id:end_id] = np.nan + out_Yorg[start_id:end_id] = np.nan + out_RXYorg[start_id:end_id] = np.nan + out_IXYorg[start_id:end_id] = np.nan + out_inc[start_id:end_id] = np.nan + out_Az[start_id:end_id] = np.nan + start_id = end_id + + #only output the useful data + out_lon = out_lon[0:end_id] + out_lat = out_lat[0:end_id] + out_Tbh = out_Tbh[0:end_id] + out_Tbv = out_Tbv[0:end_id] + out_T3 = out_T3[0:end_id] + out_T4 = out_T4[0:end_id] + out_RAh = out_RAh[0:end_id] + out_RAv = out_RAv[0:end_id] + out_RA3 = out_RA3[0:end_id] + out_RA4 = out_RA4[0:end_id] + + out_Xorg = out_Xorg[0:end_id] + out_Yorg = out_Yorg[0:end_id] + out_RXYorg = out_RXYorg[0:end_id] + out_IXYorg = out_IXYorg[0:end_id] + + out_inc = out_inc[0:end_id] + out_Az = out_Az[0:end_id] + + del dgg_list + return out_lon, out_lat, out_inc, out_Az, end_id_list, out_Tbh, out_Tbv, \ + out_T3, out_T4, out_RAh, out_RAv, out_RA3, out_RA4, \ + out_Xorg, out_Yorg, out_RXYorg, out_IXYorg diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/preprocess_nc.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/preprocess_nc.py new file mode 100644 index 0000000..0c11043 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/preprocess_nc.py @@ -0,0 +1,323 @@ +import os +import numpy as np +import numpy.matlib as matlib +from scipy.io import loadmat +from netCDF4 import Dataset +from datetime import timedelta, datetime +from src.readwrite import read_GEOSIT, read_SCLF1C_nc, write_bin_SMOS_reg +from src.helper import gal_atm_correction, xy2hv +from EASEv2 import EASEv2_latlon2ind, EASEv2_ind2latlon +import warnings; warnings.filterwarnings("ignore") + +def preprocess_nc(smos_nc, config): + + # Code to preprocess SMOS data to the M36 EASEv2 grid. + + # Input: 1 SMOS-files: *.nc + + # Output: 1) M36, angular binned Tb w/ Faraday and geometric rotation + # 2) M36, angular binned Tb w/ Faraday and geometric rotation, + # and removal of galactic, sun, moon,... + # 3) M36, angular binned Tb w/ Faraday and geometric rotation, + # and removal of galactic, sun, moon,... (aux data) + # and removal of atmospheric correction (Peng et al., 2013) + + # *** SCLF1C *** + + # Qing Liu NASA/GSFC 2021 Convert matlab code to python + # Gabrielle De Lannoy - NASA/GSFC, based on code 21sept10 + + #======================================================== + # + Apply Tx/Ty conversion to Tbh/Tbv + # + Loop through SMOS pixels, assign EASE-indices, + # add (average) them per EASE-grid + # AND average them per degree bin + # + Sequential binary output (Fortran) + + #======================================================== + # GDL, 24Nov14: - add additional output files (2, 3) + # GDL, 21Aug15: - fully remove any pointing to RWAPI-library + # - instead include xml-sparsing based on matlab routines + # - add preprocessor version in output + # - bugfix: + # + corrected spatial stdv of DGG values within EASE pixel + # + possibly corrected Tb_Sky radiation + + outpath = config['out_reg_path'] + GEOSIT_path = config['GEOSIT_path'] + GEOSIT_to_EASEv2_M36_file = config['GEOSIT_to_EASEv2_M36_file'] + SMOS_AUX_path = config['SM_OPER_AUX_GAL_SM_path'] + + overwrite=1 + + # 1 = do overwrite to get the latest processed file at the same time stamp!! + + prep_version='v102' + prep_version=int(prep_version[1:]) + + # OUTPUT: + #-------------------------------------------------------- + N_out_fields=16 + + # Auxiliary info + grid_name='M36' + + if grid_name == 'M36': + N_ease_lat=406 + N_ease_lon=964 + else: + raise RuntimeError('not ready for any other resolution than EASEv2 M36') + + write_ind_latlon='latlon' + + #print('WRITE OUT...'+write_ind_latlon) + + # INPUT SPECS: + N_sec_half_orbit=3239 + + one_time_tag=14 + all_time_tag=2*one_time_tag + 1 + 1 + K_err=75 + + max_K=375 + + max_el=150000 + N_angle=200 + + max_tot_el=max_el*N_angle + out_min_res=30.0 + + angle_step=1 + + start_angle=19.5 + + end_angle=60.5 + + inc_range=np.arange(start_angle,end_angle+0.1,angle_step) + N_angle=np.size(inc_range) - 1 + inc_angle=np.arange(start_angle + angle_step / 2.0,end_angle - angle_step / 2.0 + 0.01,angle_step) + if (N_angle != np.size(inc_angle)): + raise RuntimeError('angle bins badly defined') + + correct_galaxy_only=0 + + Asc='N' + int_Asc=np.nan + d_lon=np.full(max_tot_el,np.nan) + d_lat=np.full(max_tot_el,np.nan) + d_Tbh=np.full(max_tot_el,np.nan) + d_Tbv=np.full(max_tot_el,np.nan) + d_inc=np.full(max_tot_el,np.nan) + d_Az=np.full(max_tot_el,np.nan) + + N_dat = N_ease_lat*N_ease_lon + + data_all=np.full((N_out_fields,N_dat,N_angle),0.) + + ind_i_tmp=matlib.repmat(np.array([np.arange(N_ease_lat)]).T,1,N_ease_lon) + ind_j_tmp=matlib.repmat(np.array([np.arange(N_ease_lon)]),N_ease_lat,1) + + ind_i_all=np.reshape(ind_i_tmp,N_dat,order='F') + ind_j_all=np.reshape(ind_j_tmp,N_dat,order='F') + + N_f=8 + data_gridded_bin =np.full((N_f,N_angle,N_ease_lat,N_ease_lon),0.) + data_gridded_bin_sq=np.full((N_f,N_angle,N_ease_lat,N_ease_lon),0.) + N_data_bin =np.full((N_f,N_angle,N_ease_lat,N_ease_lon),0.) + + # 2 datasets in the SCL product: dealt with in XY2HV + #print('READING') + + # read SMOS SCLF1C netcdf file + TB_real,TB_imag,RA,Theta,Az,Fa,Ge,Snap_ID,t_smos_sec,lat,lon,flag_15,flag_16,mask_ok,Grid_ID,BT_count=read_SCLF1C_nc(smos_nc) + + fin = Dataset(smos_nc,'r') + Creator_Version = fin.getncattr('Fixed_Header:Source:Creator_Version') + Asc = fin.getncattr('Variable_Header:Specific_Product_Header:Main_Info:Time_Info:Ascending_Flag') + fname = fin.getncattr('Fixed_Header:File_Name') + + idx = fname.index('SCLF1C_') + + start_time = datetime(int(fname[idx+7:idx+11]), + int(fname[idx+11:idx+13]), + int(fname[idx+13:idx+15]), + int(fname[idx+16:idx+18]), + int(fname[idx+18:idx+20]), + int(fname[idx+20:idx+22])) + + end_time = datetime(int(fname[idx+23:idx+27]), + int(fname[idx+27:idx+29]), + int(fname[idx+29:idx+31]), + int(fname[idx+32:idx+34]), + int(fname[idx+34:idx+36]), + int(fname[idx+36:idx+38])) + + N_sec_date_diff = end_time - start_time; + + if (N_sec_date_diff.total_seconds() < 0.9*N_sec_half_orbit): + print('--This is an incomplete half orbit with '+ + str(N_sec_date_diff.total_seconds())+ + 'sec instead of '+str(N_sec_half_orbit)) + + t_mid = start_time + timedelta(seconds=np.round(N_sec_date_diff.total_seconds()/2)) + t_out = t_mid.replace(second=0) + + if np.round(t_mid.minute/out_min_res) < 1: + t_out = t_out.replace(minute = 0) + elif np.round(t_mid.minute/out_min_res) == 1: + t_out = t_out.replace(minute = 30) + else: + t_out = t_out.replace(minute = 0) + t_out = t_out + timedelta(seconds=3600); + + if int(Creator_Version) >= 344: + alpha=Fa + Ge + else: + alpha=Ge - Fa + + del Fa,Ge + + d_lon,d_lat,d_inc,d_Az,end_id,d_Tbh,\ + d_Tbv,d_T3,d_T4,d_RAh,d_RAv,d_RA3,d_RA4,\ + d_Xorg,d_Yorg,d_RXYorg,d_IXYorg = xy2hv \ + (TB_real,TB_imag,RA,Theta,Az,alpha,Snap_ID,t_smos_sec,\ + lat,lon,flag_15,flag_16,mask_ok,Grid_ID,BT_count,nargout=17) + + d_Tbh[d_Xorg == 0]=np.nan + d_Tbv[d_Yorg == 0]=np.nan + d_T3[d_RXYorg == 0]=np.nan + d_T4[d_IXYorg == 0]=np.nan + + dat_l=max(end_id) + + if Asc =='A': + int_Asc=1 + else: + int_Asc=0 + + N_points=len(d_Tbh) + check=len(d_Tbv) + + if ((N_points != check) or (N_points != dat_l)): + raise RuntimeError('error in variable dimensions coming out of rotation') + + #print('Working on '+str(N_points)+' data points ; max='+str(max_tot_el)) + C=1 + + # obtain GEOSIT latlon -> EASEv2 grid remapping info: + # the mat file is pre-generated usign the "aux_grids_SMOS_prep" matlab script + + m2d = loadmat(GEOSIT_to_EASEv2_M36_file) + vegcls_grid = m2d['vegcls_grid'] + NN_grid = m2d['NN_grid'] + + T_air,T_surf,V_surf,P_surf = read_GEOSIT(t_out,vegcls_grid,NN_grid,d_lat, d_lon,GEOSIT_path) + + Tb_ap_H,Tb_ap_V,Tb_BOA_H,Tb_BOA_V = \ + gal_atm_correction(d_Tbh,d_Tbv,d_lat,d_lon,d_inc,d_Az,\ + t_out,SMOS_AUX_path,T_air,T_surf,P_surf,V_surf,C,nargout=4) + + #if correct_galaxy_only: + # print('Galact corr (H): maxdiff '+str(np.nanmax(Tb_ap_H - d_Tbh))+' mindiff '+str(np.nanmin(Tb_ap_H - d_Tbh))) + #else: + # print('Finished galact+atm correction') + #print('Galact+atm corr (H): maxdiff '+str(np.nanmax(Tb_BOA_H - d_Tbh))+' mindiff '+str(np.nanmin(Tb_BOA_H - d_Tbh))) + + #print('Project the SMOS data on a M36 grid...') + ind_row,ind_col=EASEv2_latlon2ind(d_lat,d_lon,grid_name,nargout=2) + if correct_galaxy_only: + Tb_data=np.array([Tb_ap_H[np.arange(N_points)], + Tb_ap_V[np.arange(N_points)], + d_RAh[ np.arange(N_points)], + d_RAv[ np.arange(N_points)], + np.nan*d_T3[ np.arange(N_points)], + np.nan*d_T4[ np.arange(N_points)], + np.nan*d_RA3[np.arange(N_points)], + np.nan*d_RA4[np.arange(N_points)]]) + else: + Tb_data=np.array([Tb_BOA_H[np.arange(N_points)], + Tb_BOA_V[np.arange(N_points)], + d_RAh[ np.arange(N_points)], + d_RAv[ np.arange(N_points)], + np.nan*d_T3[ np.arange(N_points)], + np.nan*d_T4[ np.arange(N_points)], + np.nan*d_RA3[np.arange(N_points)], + np.nan*d_RA4[np.arange(N_points)]]) + + coord=np.array([ind_col[np.arange(N_points)], + ind_row[np.arange(N_points)], + d_inc[np.arange(N_points)]]) + + for i in np.arange(N_points): + ind_j=coord[0,i].astype('int') + ind_i=coord[1,i].astype('int') + ind_a=np.ceil((coord[2,i] - start_angle) / angle_step) + ind_a=ind_a.astype('int') + + if np.any(np.logical_not(np.isnan(Tb_data[:,i]))): + if ind_a > 0 and ind_a <= N_angle and \ + np.nanmax(Tb_data[:,i]) < max_K: + ind_f=np.nonzero(np.logical_not(np.isnan(Tb_data[:,i])))[0] + data_gridded_bin[ind_f,ind_a-1,ind_i,ind_j]=data_gridded_bin[ind_f,ind_a-1,ind_i,ind_j] + Tb_data[ind_f,i] + data_gridded_bin_sq[ind_f,ind_a-1,ind_i,ind_j]=data_gridded_bin_sq[ind_f,ind_a-1,ind_i,ind_j] + Tb_data[ind_f,i]*Tb_data[ind_f,i] + N_data_bin[ind_f,ind_a-1,ind_i,ind_j]=N_data_bin[ind_f,ind_a-1,ind_i,ind_j] + 1 + + # Could refine the angle bins, filter out more RFI, i.e. + # using K_err as maximum diff from mean, + # and then reshape the data afterwards to get less bins + N_data_bin[N_data_bin==0] = np.nan + data_gridded_bin=data_gridded_bin / N_data_bin + data_gridded_bin_sq=data_gridded_bin_sq / N_data_bin + data_gridded_bin_sq=np.sqrt((data_gridded_bin_sq - data_gridded_bin ** 2)*N_data_bin/ (N_data_bin - 1)) + + # Set #points to zero, so that RFI-affected Tb + # won't show up in output later. + N_data_bin[data_gridded_bin_sq > K_err]=0 + + for a in np.arange(N_angle): + data_all[ 0,:,a]=np.reshape(data_gridded_bin[0,a,:,:],[N_dat],order='F') + data_all[ 1,:,a]=np.reshape(data_gridded_bin[1,a,:,:],[N_dat],order='F') + data_all[ 2,:,a]=np.reshape(data_gridded_bin_sq[0,a,:,:],[N_dat],order='F') + data_all[ 3,:,a]=np.reshape(data_gridded_bin_sq[1,a,:,:],[N_dat],order='F') + data_all[ 4,:,a]=np.reshape(N_data_bin[0,a,:,:],[N_dat],order='F') + data_all[ 5,:,a]=np.reshape(N_data_bin[1,a,:,:],[N_dat],order='F') + data_all[ 6,:,a]=np.reshape(data_gridded_bin[2,a,:,:],[N_dat],order='F') + data_all[ 7,:,a]=np.reshape(data_gridded_bin[3,a,:,:],[N_dat],order='F') + data_all[ 8,:,a]=np.reshape(data_gridded_bin[4,a,:,:],[N_dat],order='F') + data_all[ 9,:,a]=np.reshape(data_gridded_bin[5,a,:,:],[N_dat],order='F') + data_all[10,:,a]=np.reshape(data_gridded_bin_sq[4,a,:,:],[N_dat],order='F') + data_all[11,:,a]=np.reshape(data_gridded_bin_sq[5,a,:,:],[N_dat],order='F') + data_all[12,:,a]=np.reshape(N_data_bin[4,a,:,:],[N_dat],order='F') + data_all[13,:,a]=np.reshape(N_data_bin[5,a,:,:],[N_dat],order='F') + data_all[14,:,a]=np.reshape(data_gridded_bin[6,a,:,:],[N_dat],order='F') + data_all[15,:,a]=np.reshape(data_gridded_bin[7,a,:,:],[N_dat],order='F') + + ind=np.nonzero(np.any(data_all[4:6,:,:] > 0, axis=(0,2)))[0] + + data_all[np.isnan(data_all)]= -999.0 + stamp1=t_out.strftime("%Y") + stamp2=t_out.strftime("%m") + stamp3=t_out.strftime("%d") + stamp4=t_out.strftime("%H") + stamp5=t_out.strftime("%M") + + if not os.path.exists(outpath+'/'+stamp1+stamp2): + os.makedirs(outpath+'/'+stamp1+stamp2) + + if correct_galaxy_only: + out_filename=outpath+'/'+stamp1+stamp2+'/SMOS_reg_nosky_Tb_'+stamp1+stamp2+stamp3+'_'+stamp4+stamp5+'_'+Asc+'.bin' + else: + out_filename=outpath+'/'+stamp1+stamp2+'/SMOS_reg_Tb_'+stamp1+stamp2+stamp3+'_'+stamp4+stamp5+'_'+Asc+'.bin' + + if ind.size > 0: + if write_ind_latlon == 'ind': + write_bin_SMOS_reg(out_filename,ind_j_all[ind],ind_i_all[ind],inc_angle,data_all[:,ind,:],int_Asc,Creator_Version,prep_version,start_time,end_time,overwrite,N_out_fields,write_ind_latlon,'SCLF1C') + else: + if write_ind_latlon == 'latlon': + lat_out,lon_out=EASEv2_ind2latlon(ind_i_all[ind],ind_j_all[ind],grid_name,nargout=2) + write_bin_SMOS_reg(out_filename,lon_out,lat_out,inc_angle,data_all[:,ind,:],int_Asc,int(Creator_Version),prep_version,start_time,end_time,overwrite,N_out_fields,write_ind_latlon,'SCLF1C') + else: + raise RuntimeError('Unknown format of indices / latlon output') + #else: + # print('Zero valid data, NO output file written') + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/__init__.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/__init__.py new file mode 100644 index 0000000..4d7f3ee --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/__init__.py @@ -0,0 +1,5 @@ +from .read_aux_gal_sm import read_aux_gal_sm +from .read_GEOSIT import read_GEOSIT +from .write_bin_SMOS_reg import write_bin_SMOS_reg +from .read_bin_SMOS_reg import read_bin_SMOS_reg +from .read_SCLF1C_nc import read_SCLF1C_nc diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_GEOSIT.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_GEOSIT.py new file mode 100644 index 0000000..a7981e3 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_GEOSIT.py @@ -0,0 +1,114 @@ +import numpy as np +from netCDF4 import Dataset +from EASEv2 import EASEv2_latlon2ind + +def read_GEOSIT(date_time=None,vegcls_grid=None,EASEv2_M36_to_M2_grid=None,lat=None,lon=None,GEOSIT_path=None): + + # Read fields from a GEOSIT + + # INPUT: + # date_time = #structure with year, month, day, hour,... + # vegcls_grid = #2-d vegcls on EASEv2_M36 output grid + # EASEv2_M36_to_M2_grid = #2-d EASEv2_M36 grid where each index + # refers to the nearest GEOSIT LL grid cell + # lat = [ : ]; #EASEv2 M36 indices for each element in swath + # lon = [ : ]; + + # OUTPUT: + # Tair, Ts, Vs, Ps + # - in units needed for atm and gal correction + # - for locations within SMOS swath only + #----------------------------------------------------------------- + + C_2_K=273.15 + + if (date_time.year < 2008): + print('did not look into the stream sequence for older years than 2010') + else: + if (date_time.year < 2018): + stream='d5294_geosit_jan08' + else: + stream='d5294_geosit_jan18' + + #GEOS-IT: + + path=GEOSIT_path+stream+'/diag/Y'+date_time.strftime("%Y")+'/M'+date_time.strftime("%m")+'/' + + date_string=date_time.strftime("%Y-%m-%dT%H") + + filename= Dataset(path+'/'+stream+'.lnd_tavg_1hr_glo_L576x361_slv.'+date_string+'30Z.nc4','r') + + Tp = filename.variables['TSURF'] + + Ts = filename.variables['TSOIL1'] + + filename = Dataset(path+'/'+stream+'.slv_tavg_1hr_glo_L576x361_slv.'+date_string+'30Z.nc4','r') + + Ps = filename.variables['PS'] + + SH = filename.variables['QV2M'] + + Tair = filename.variables['T2M'] + + #print('Read GEOSIT files - '+stream+'.lnd_tavg_1hr_glo_L576x361_slv.'+date_string+'30Z.nc4') + # Interpolate, limit to swath + #----------------------------------------------------------------- + + row_ind,col_ind = EASEv2_latlon2ind(lat,lon,'M36') + idx_M36 = np.ravel_multi_index([row_ind,col_ind],EASEv2_M36_to_M2_grid.shape, mode='raise', order='F') + #idx = sub2ind(EASEv2_M36_to_M2_grid.shape,row_ind,col_ind) + + idx = EASEv2_M36_to_M2_grid.ravel(order='F')[idx_M36] + idx = idx -1 # change to 0-based index + + #--- 0) Surface air temperature in [^o C]------------------------- + + Tair = Tair[0,:,:]- C_2_K + Tair = Tair.flatten('F')[idx] + #--- 1) Surface temperature in [^o C]----------------------------- + + Ts = Ts[0,:,:] - C_2_K + Tp = Tp[0,:,:] - C_2_K + + Ts=Ts.flatten('F')[idx] + Tp=Tp.flatten('F')[idx] + Ts[vegcls_grid.flatten('F')[idx_M36] != 1] = \ + (Tp[vegcls_grid.flatten('F')[idx_M36] != 1] + \ + Ts[vegcls_grid.flatten('F')[idx_M36] != 1]) / 2.0 + + if any(Ts > 70) or any(Ts < - 80): + print('WARNING: Surface temperature out of range') + + #print('Surface temperature min '+str(Ts.min())+' max '+str(Ts.max())+' [^o C]') + #--- 2) Surface water vapour density [g/m3]------------------------ + + # SH = specific humidity, ratio of the water vapor content of the + # mixture to the total air content on a mass basis. + # Assuming water density at 1 kg/m^3: + + Vs = SH[0,:,:] * 10**3 + + Vs=Vs.flatten('F')[idx] + + if any(Vs > 80): + print('WARNING: Surface vapour pressure density out of range') + + if any(Vs < 0): + print('WARNING: min vapour pressure density will be reset. min: '+min(Vs)) + Vs[Vs < 0] = 0 + + #print('Surface water vapour density min '+str(Vs.min())+' max '+str(Vs.max())+' [g/m3]') + #--- 3) Surface pressure [mbar]------------------------------------ + + Ps = Ps[0,:,:] * 10**(-5) * 10**3 # [Pa] --> [mbar] + + + Ps = Ps.flatten('F')[idx] + + if any(Ps > 1200) or any(Ps < 400): + print('WARNING: Surface pressure out of range') + + #print('Surface pressure min '+str(Ps.min())+' max '+str(Ps.max())+' [mbar]') + #----------------------------------------------------------------- + return Tair, Ts, Vs, Ps + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_SCLF1C_nc.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_SCLF1C_nc.py new file mode 100644 index 0000000..e38b54c --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_SCLF1C_nc.py @@ -0,0 +1,83 @@ +from netCDF4 import Dataset +import numpy as np + +def isKthBitSet(number, kth): + if number & (1 << (kth - 1)): + return True + else: + return False + +def bitget(number, kth): + return (number >> kth-1) & 1 + +def ismember(a, b): + bind = {} + for i, elt in enumerate(b): + if elt not in bind: + bind[elt] = i + # None can be replaced by any other "not in b" value + return [bind.get(itm, None) for itm in a] + +def read_SCLF1C_nc(fname=None,*args,**kwargs): + + fin = Dataset(fname,'r') + + lat = fin.variables['Grid_Point_Latitude'][:] + lon = fin.variables['Grid_Point_Longitude'][:] + + mask_grid = fin.variables['Grid_Point_Mask'][:] + + mask_ok = np.logical_and(bitget(mask_grid,2), np.logical_not(bitget(mask_grid,6))) + + Grid_ID = fin.variables['Grid_Point_ID'][:] + flag_grid = fin.variables['Flags'][:] + + RA = fin.variables['Radiometric_Accuracy_of_Pixel'][:] + + Theta = 1.*fin.variables['Incidence_Angle'][:] + Az = fin.variables['Azimuth_Angle'][:] + Fa = fin.variables['Faraday_Rotation_Angle'][:] + Ge = fin.variables['Geometric_Rotation_Angle'][:] + + Seconds = fin.variables['Seconds'][:] + Microseconds = fin.variables['Microseconds'][:] + + BT_Value_Real = fin.variables['BT_Value_Real'][:] + BT_Value_Imag = fin.variables['BT_Value_Imag'][:] + + # v620 + #mask_nonRFI = np.all((np.logical_not(bitget(flag_grid,12)),np.logical_not(bitget(flag_grid,16)), \ + # bitget(flag_grid,11), BT_Value_Real < 325., BT_Value_Imag < 325.), axis=0) + # v724 + mask_nonRFI = np.all((np.logical_not(bitget(flag_grid,12)),np.logical_not(bitget(flag_grid,7)), \ + bitget(flag_grid,11), BT_Value_Real < 325., BT_Value_Imag < 325.), axis=0) + + TB_real = BT_Value_Real + TB_imag = BT_Value_Imag + + TB_real[np.logical_not(mask_nonRFI)] = np.nan + TB_imag[np.logical_not(mask_nonRFI)] = np.nan + + RA[np.logical_not(mask_nonRFI)] = np.nan + Theta[np.logical_not(mask_nonRFI)] = np.nan + Az[np.logical_not(mask_nonRFI)] = np.nan + Fa[np.logical_not(mask_nonRFI)] = np.nan + Ge[np.logical_not(mask_nonRFI)] = np.nan + + Snap_ID = fin.variables['Snapshot_ID'][:] + Snap_ID_of_Pixel = fin.variables['Snapshot_ID_of_Pixel'][:] + + idx_nonRFI = np.nonzero(mask_nonRFI) + t_smos_sec = np.full(TB_real.shape,np.nan) + loc_snap = ismember(Snap_ID_of_Pixel[idx_nonRFI], Snap_ID) + t_smos_sec[idx_nonRFI] = Seconds[loc_snap] + Microseconds[loc_snap] / 1.e6 + + flag_16 = bitget(flag_grid,1) + flag_15 = bitget(flag_grid,2) + + BT_count = np.sum(np.logical_not(np.isnan(TB_real)).astype(int),axis=1) + + return TB_real, TB_imag, RA, Theta, Az, Fa, Ge, Snap_ID_of_Pixel, t_smos_sec, \ + lat, lon, flag_15, flag_16, mask_ok, Grid_ID, BT_count + + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_aux_gal_sm.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_aux_gal_sm.py new file mode 100644 index 0000000..a063a1b --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/read_aux_gal_sm.py @@ -0,0 +1,11 @@ +import numpy as np + +def read_aux_gal_sm(filename=None): + N_row=721 + N_col=1441 + tmp = np.fromfile(filename, dtype=' repeated for T3 and T4 (9-16) + +#SM-files: +#------------------------------------------------------------------ +#N_out_fields # 1 - SM [N float] + # 2 - ST [N float] + # 3 - tau [N float] + # 4 - Tbh 42.5^o simulated on antenna reference frame [N float] + # 5 - Tbv 42.5^o simulated on antenna reference frame [N float] + # 6 - RSTDSM [N float] + # 7 - RSTDST [N float] + # 8 - RSTDtau [N float] + # 9 - stdv in SM EASE pixel (heterogeneity index) [N float] + # 10 - number of SMOS SM pixels per EASE grid cell [N int] + # 11 - science flag [N int] + +# ------------------------------------------------------------------ + + print('reading from '+fname) + + if path.isfile(fname): + with open(fname,'rb') as fin: + din=fin.read() + fin.close() + else: + sys.exit('file does not exist') + + # fortran tag before and after each record + byte_beg = 0; byte_end = 4 + + # byte size for record + byte_beg = byte_end + byte_end = byte_beg + 4*3 + asc_flag,version,prep_version = struct.unpack('>'+'i'*3,din[byte_beg:byte_end]) + + # 2 fortran tags between records + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + byte_end = byte_beg + 4*5 + + tmp = struct.unpack('>'+'i'*5, din[byte_beg:byte_end]) + start_time = datetime(tmp[0],tmp[1],tmp[2],tmp[3],tmp[4],0) + + # 2 fortran tags + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + byte_end = byte_beg + 4*5 + + tmp = struct.unpack('>'+'i'*5, din[byte_beg:byte_end]) + end_time = datetime(tmp[0],tmp[1],tmp[2],tmp[3],tmp[4],0) + + # 2 fortran tags + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + if data_product == 'scaling': + byte_end = byte_beg + 4*3 + N_grid,N_angle,N_tile = struct.uppack('>'+'i'*3, din[byte_beg:byte_end]) + else: + byte_end = byte_beg + 4*2 + N_grid,N_angle = struct.unpack('>'+'i'*2, din[byte_beg:byte_end]) + N_tile=1 + + print('N_grid and N_angle : '+ str(N_grid) +', '+ str(N_angle)) + + # read all records + + # 2 fortran tags + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + if (N_grid > 1): + byte_end = byte_beg + 4*N_angle + inc_angle = struct.unpack('>'+'f'*N_angle,din[byte_beg:byte_end]) + + # 2 fortran tags + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + byte_end = byte_beg + 4*N_grid + if read_ind_latlon == 'ind': + col_ind=struct.unpack('>'+'i'*N_grid,din[byte_beg:byte_end]) + + # 2 fortran tags + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + byte_end = byte_beg + 4*N_grid + + row_ind=struct.unpack('>'+'i'*N_grid,din[byte_beg:byte_end]) + else: + col_ind=struct.unpack('>'+'f'*N_grid,din[byte_beg:byte_end]) + + # 2 fortran tags + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + byte_end = byte_beg + 4*N_grid + + row_ind=struct.unpack('>'+'f'*N_grid,din[byte_beg:byte_end]) + + if read_ind_latlon == 'latlon_id': + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + byte_end = byte_beg + 4*N_grid + tile_id = struct.unpack('>'+'f'*N_grid,din[byte_beg:byte_end]) + + + SMOS_data=np.full([N_out_fields,N_grid,N_angle],np.nan) + + for i in range(N_out_fields): + for j in range(N_angle): + byte_beg = byte_end; byte_end = byte_beg + 4*2 + + byte_beg = byte_end + byte_end = byte_beg + 4*N_grid + if (((i == 4 or i == 5 or i == 11 or i == 13) and \ + (data_product =='SCLF1C' or data_product == 'BWLF1C')) or \ + (i == 9 and data_product == 'SMUDP2') or \ + ((i == 4 or i == 9) and data_product == 'scaling')): + + tmp_data=struct.unpack('>'+'i'*N_grid,din[byte_beg:byte_end]) + + else: + + tmp_data=struct.unpack('>'+'f'*N_grid,din[byte_beg:byte_end]) + + SMOS_data[i,:,j]=tmp_data + + else: + SMOS_data=np.full([N_out_fields],np.nan) + col_ind=np.nan + row_ind=np.nan + inc_angle=np.nan + + return np.array(SMOS_data), np.array(inc_angle), np.array(col_ind), np.array(row_ind), asc_flag, \ + version, prep_version, start_time, end_time, N_grid + diff --git a/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/write_bin_SMOS_reg.py b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/write_bin_SMOS_reg.py new file mode 100644 index 0000000..dd71047 --- /dev/null +++ b/GEOSldas_App/util/inputs/obs_preproc/SMOS_preproc/src/readwrite/write_bin_SMOS_reg.py @@ -0,0 +1,251 @@ +import numpy as np +import sys +from os import path +import struct + +def write_bin_SMOS_reg(fname=None,colind=None,rowind=None,av_angle_bin=None, \ + data=None,asc_flag=None,version=None,version_prep=None,\ + start_time=None,end_time=None,overwrite=False,N_out_fields=None,\ + write_ind_latlon=None,data_product=None,tile_id=None,*args,**kwargs): + +# write "fortran sequential" tile tavg files (identical to LDASsa output) + +# optional input: +# overwrite = 0 -- do NOT overwrite existing files, print warning +# message, return +# overwrite = 1 -- overwrite existing files, print warning message + +# this py function is nearly identical to the matlab function that writes the Tb scaling parameters: +# ./util/inputs/obs_scaling_params/write_seqbin_file.m +# *except* this py function writes one additional field specifically for SMOS binary files + +# ------------------------------------------------------------------ + + #N_out_fields # 1 - Col-index, 0-based; + # 2 - Row-index, 0-based; + #OR (for nearest neighbout) + # 1 - Lon; + # 2 - Lat; + + #N_out_fields # 1 - Tbh; + # 2 - Tbv; + + # 3 - heterogeneity index Tbh + # 4 - heterogeneity index Tbv + + # 5 - # SMOS pixels in EASE grid pixel Tbh + # 6 - # SMOS pixels in EASE grid pixel Tbv + + # 7 - RA Tbh + # 8 - RA Tbv + + #=> repeated for T3 and T4 (9-16) + + #OR FOR SMUDP2: + + # 1 - SM + # 2 - ST + # 3 - opacity + # 4 - Tbh; + # 5 - Tbv; + + # 6 - SM RSTD + # 7 - ST RSTD + # 8 - opac RSTD + + # 9 - stdv in SM (grid cell averaging) + # 10 - # small SMOS pixels inside 1 EASE grid cell; + + # 11 - accumulated flag + + #NOT 11 - omega; scattering albedo + #NOT 12 - diff_albedos (om_H-om_V) + #NOT 13 - max_roughness + #NOT 14 - RSTD omega + #NOT 15 - RSTD diff_omega + #NOT 16 - RSTD max_roughness + + + # check dimensions + if data.shape[0] != N_out_fields: + sys.exit('ERROR: size of data incompatible with N_out_fields') + + # check if file exists + if path.isfile(fname): + if not overwrite: + sys.exit('RETURNING!!! -- NOT OVERWRITING EXISTING FILE '+fname) + return + else: + print('OVERWRITING '+fname) + else: + print('writing '+fname) + + N_grid= data.shape[1] + N_angle=1 + if (len(data.shape) == 3): + N_angle=data.shape[2] + data_org=data + if (N_angle != av_angle_bin.size): + sys.exit('ERROR in N_angle') + + if ( write_ind_latlon =='latlon_id' and len(args) == 14): + if tile_id.shape[0] != N_grid: + sys.exit('tile_id dimensions ??') + if tile_id.shape[1] > 1: + print('# subgridcells per gridcell: '+ str(tile_id.shape[1])) + + # open file + fout=open(fname,'wb') + # determine number of grid cells ; further check dimensions + + # write all records + # length of each reoord in bytes + fortran_tag= 3*4 + + # write output, '>' lead the format str force big endian order + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*3,asc_flag,version,version_prep)) + fout.write(struct.pack('>i',fortran_tag)) + + fortran_tag=5*4 + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*5, int(start_time.year),int(start_time.month),int(start_time.day),int(start_time.hour),int(start_time.minute))) + fout.write(struct.pack('>i',fortran_tag)) + + fortran_tag=5*4 + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*5, int(end_time.year),int(end_time.month),int(end_time.day),int(end_time.hour),int(end_time.minute))) + fout.write(struct.pack('>i',fortran_tag)) + + if not (data_product == 'scaling' and write_ind_latlon =='latlon_id' and len(args) == 14): + fortran_tag=2*4 + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*2,N_grid,N_angle)) + fout.write(struct.pack('>i',fortran_tag)) + else: + fortran_tag=3*4 + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*3,N_grid,N_angle,tile_id.shape[1])) + fout.write(struct.pack('>i',fortran_tag)) + + if N_grid >= 1: + fortran_tag=N_angle*4 + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'f'*N_angle,*av_angle_bin)) + fout.write(struct.pack('>i',fortran_tag)) + + fortran_tag=N_grid*4 + if write_ind_latlon == 'ind': + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*N_grid,np.round(colind))) + fout.write(struct.pack('>i',fortran_tag)) + + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*N_grid,np.round(rowind))) + fout.write(struct.pack('>i',fortran_tag)) + + elif write_ind_latlon == 'latlon': + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'f'*N_grid, *colind)) + fout.write(struct.pack('>i',fortran_tag)) + + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'f'*N_grid, *rowind)) + fout.write(struct.pack('>i',fortran_tag)) + + elif write_ind_latlon == 'latlon_id' and len(args) == 14: + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'f'*N_grid, *colind)) + fout.write(struct.pack('>i',fortran_tag)) + + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'f'*N_grid, *rowind)) + fout.write(struct.pack('>i',fortran_tag)) + + for i in range(tile_id.shape[1]): + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*N_grid, *tile_id[:][i])) + fout.write(struct.pack('>i',fortran_tag)) + else: + sys.exit('output-arguments do not line up') + + fortran_tag=N_grid*4 + for i in range(N_out_fields): + for jj in range(N_angle): + if N_angle > 1: + data=data_org[:,:,jj] + + if ((i == 4 or i == 5 or i == 12 or i == 13) and (data_product == 'SCLF1C' or data_product == 'BWLF1C')) \ + or (i == 9 and data_product == 'SMUDP2'): + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*N_grid, *data[i,:].round().astype('int'))) + fout.write(struct.pack('>i',fortran_tag)) + + elif (i == 4 or i == 9) and data_product == 'scaling': + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'i'*N_grid, *data[i,:].round().astype('int'))) + fout.write(struct.pack('>i',fortran_tag)) + else: + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'f'*N_grid, *data[i,:])) + fout.write(struct.pack('>i',fortran_tag)) + + else: + fortran_tag=N_angle*4 + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>'+'f'*N_angle,*av_angle_bin)) + fout.write(struct.pack('>i',fortran_tag)) + + fortran_tag=4 + if write_ind_latlon == 'ind': + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>i',0)) + fout.write(struct.pack('>i',fortran_tag)) + + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>i',0)) + fout.write(struct.pack('>i',fortran_tag)) + + elif write_ind_latlon == 'latlon': + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>f', 0.)) + fout.write(struct.pack('>i',fortran_tag)) + + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>f', 0.)) + fout.write(struct.pack('>i',fortran_tag)) + + elif write_ind_latlon == 'latlon_id' and len(args) == 14: + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>f', 0.)) + fout.write(struct.pack('>i',fortran_tag)) + + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>f', 0.)) + fout.write(struct.pack('>i',fortran_tag)) + + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>i', 0)) + fout.write(struct.pack('>i',fortran_tag)) + + else: + sys.exit('output-arguments do not line up') + + for i in range(N_out_fields): + for jj in range(N_angle): + if ((i == 4 or i == 5 or i == 12 or i == 13) and ( data_product == 'SCLF1C' or data_product == 'BWLF1C')) \ + or (i == 9 and data_product == 'SMUDP2'): + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>i', 0)) + fout.write(struct.pack('>i',fortran_tag)) + + elif (i == 4 or i == 9) and data_product == 'scaling': + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>i', 0)) + fout.write(struct.pack('>i',fortran_tag)) + else: + fout.write(struct.pack('>i',fortran_tag)) + fout.write(struct.pack('>f', -999.0)) + fout.write(struct.pack('>i',fortran_tag)) + + fout.close() diff --git a/GEOSldas_App/util/inputs/obs_scaling_params/write_seqbin_file.m b/GEOSldas_App/util/inputs/obs_scaling_params/write_seqbin_file.m index 1415372..1c6246a 100644 --- a/GEOSldas_App/util/inputs/obs_scaling_params/write_seqbin_file.m +++ b/GEOSldas_App/util/inputs/obs_scaling_params/write_seqbin_file.m @@ -13,6 +13,10 @@ % message, return % overwrite = 1 -- overwrite existing files, print warning message % +% this matlab function is nearly identical to the py function that writes the preprocessed SMOS obs: +% ./util/inputs/obs_preproc/SMOS_preproc/src/readwrite/write_bin_SMOS_reg.py +% *except* that the py function writes one additional field specifically for SMOS binary files +% % De Lannoy, 4 Oct 2010 % De Lannoy, 26 Sep 2012: added optional argument of tile_id % used to write scaling files, with ''latlon_id''. diff --git a/GEOSldas_App/util/shared/python/EASEv2.py b/GEOSldas_App/util/shared/python/EASEv2.py index 31cd439..b6f66eb 100644 --- a/GEOSldas_App/util/shared/python/EASEv2.py +++ b/GEOSldas_App/util/shared/python/EASEv2.py @@ -181,6 +181,10 @@ def EASEv2_ind2latlon(row=None,col=None,gridid=None,*args,**kwargs): ((((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) + + lat = np.atleast_1d(np.asarray(lat, dtype=np.float64)) + lon = np.atleast_1d(np.asarray(lon, dtype=np.float64)) + msk1= np.where(lon < - 180.0) lon[msk1]= lon[msk1] + 360.0 msk2= np.where(lon > 180.0) @@ -195,5 +199,218 @@ def EASEv2_ind2latlon(row=None,col=None,gridid=None,*args,**kwargs): lon[idx]= float('nan') return lat, lon + +# -------------------------------------------------------------------------------------------------------- + +# 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 + +def EASEv2_latlon2ind(lat=None,lon=None,gridid=None,*args,**kwargs): + lat = np.atleast_1d(np.asarray(lat, dtype=np.float64)) + lon = np.atleast_1d(np.asarray(lon, dtype=np.float64)) + # 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.') -# ============= EOF ============================================ + # Constants returned by EASE2_MAP_INFO.PRO + epsilon=1e-06 + 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_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*np.pi / 180.0 + lam= dlon*np.pi / 180.0 + sin_phi=np.sin(phi) + q=(1.0 - e2)*((sin_phi / (1.0 - e2*sin_phi*sin_phi)) - (1.0 / (2.0*map_eccentricity))*np.log((1.0 - map_eccentricity * sin_phi) / (1.0 + map_eccentricity * sin_phi))) + qp=1.0 - ((1.0 - e2) / (2.0*map_eccentricity)*np.log((1.0 - map_eccentricity) / (1.0 + map_eccentricity))) + if 'M' == projection: + x= map_equatorial_radius_m * kz * lam + y= (map_equatorial_radius_m * q) / (2.0*kz) + elif 'N' == projection: + tmp=qp - q + tmp[np.absolute(tmp) < epsilon]=0.0 + rho= map_equatorial_radius_m*np.sqrt(tmp) + x= rho * np.sin(lam) + y= -rho * np.cos(lam) + elif 'S' == projection: + tmp=qp + q + tmp[np.absolute(tmp) < epsilon]=0.0 + rho= map_equatorial_radius_m * np.sqrt(tmp) + x= rho * np.sin(lam) + y= rho * np.cos(lam) + + row=s0 - (y / map_scale_m) + col=r0 + (x / map_scale_m) + + if 'N' == projection: + idx=(lat < 0.0) + row[idx]=np.nan + col[idx]=np.nan + elif 'S' == projection: + idx=(lat > 0.0) + row[idx]=np.nan + col[idx]=np.nan + + # restrict indices to (zero-based) min and max indices + # (needed because of round-off error - GDL, 07 Jun 2013) + + row=np.maximum(0,np.minimum(rows - 1,np.round(row))) + col=np.maximum(0,np.minimum(cols - 1,np.round(col))) + + # assign max index if nan (unclear why. this is based on the matlab code) + row[(np.isnan(row))] = rows - 1 + col[(np.isnan(col))] = cols - 1 + return row.astype('int'),col.astype('int') + + +# ============= EOF ========================================================================================