diff --git a/lib/config/polly_global_config.json b/lib/config/polly_global_config.json index 4a62cbbb..5e4212eb 100644 --- a/lib/config/polly_global_config.json +++ b/lib/config/polly_global_config.json @@ -16,17 +16,38 @@ "flagMolDepolCali": false, "flagTransCor": true, "flagUseTheoreticalMDR": false, - - "MWRFolder": "", + "flagUseLatestGDAS": true, "dataFileFormat": "(?\\d{4})_(?\\d{2})_(?\\d{2})_\\w*_(?\\d{2})_(?\\d{2})_(?\\d{2})\\w*.nc", "gdas1Site": "", + + "meteorDataSource": "nc_cloudnet", + "meteo_folder": "/lacroshome/cloudnet/data/mindelo/calibrated/ecmwf", "meteorDataSource": "gdas1", "meteo_folder": "/data/level1a/model/gdas1/profiles", "AERONETSite": "", + "radiosondeType": 1, + "radiosondeSitenum": 14430, + "IWV_instrument": "mwr_cloudnet", + "maxIWVTLag": 0.0833, + "MWRFolder": "/data/level1b/cloudnetpy/products/SITENAME", + "angstrexp_NR": 1.0, + "angstrexp": 1.0, "max_height_bin": 3000, - "first_range_gate_indx": [251, 251, 251, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252], + "cloudScreenMode": 1, + "maxSigSlope4FilterCloud": 0.7e5, + "maxSigSlope4FilterCloud_NR": 0.7e5, + + "overlapCalMode": 2, + "overlapCorMode": 2, + "overlapSmoothBins": 8, + + "intNProfiles": 120, + "minIntNProfiles": 60, + + "init_depAng": 0, + "first_range_gate_indx": [251, 251, 251, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252], "first_range_gate_height": 3.75, "deltaT": 30, "dtCorModeLabel": ["Polynomial coefficients stores in netcdf file", "nonparalyzable correction", "User defined polynomial cofficients"], @@ -50,25 +71,33 @@ "bgCorRangeIndx": [10, 240], "mask_SNRmin": [0.5, 0.5, 0.5, 0.01, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5], "tempCorFunc": ["1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1"], - "init_depAng": 999, - "depol_cali_mode": 2, + + "depol_cali_mode": 1, + "TR": [1.04, 500, 1, 1, 1.01, 1500, 1, 1, 1, 1, 1, 1, 1, 1, 500, 1, 1, 1, 1, 500, 1], + "G": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], + "H": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], + "K": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], + "voldepol_error_355": [-999, -999, -999 ], + "voldepol_error_532": [-999, -999, -999 ], + "voldepol_error_1064": [-999, -999, -999 ], "depol_cal_time_fixed_p_start": [""], "depol_cal_time_fixed_p_end": [""], "depol_cal_time_fixed_m_start": [""], "depol_cal_time_fixed_m_end": [""], "maskDepCalAng": ["none", "none", "p", "p", "p", "p", "p", "p", "p", "p", "none", "n", "n", "n", "n", "n", "n", "n", "n"], + "depol_cal_minbin_532": 100, "depol_cal_maxbin_532": 300, "depol_cal_SNRmin_532": [1, 1, 1, 1], - "depol_cal_sigMax_532": [1500, 1500, 1500, 1500], + "depol_cal_sigMax_532": [2000, 2000, 2000, 2000], "rel_std_dplus_532": 0.2, "rel_std_dminus_532": 0.2, "depol_cal_segmentLen_532": 40, "depol_cal_smoothWin_532": 8, "depol_cal_minbin_355": 100, "depol_cal_maxbin_355": 300, - "depol_cal_SNRmin_355": [2, 2, 2, 2], - "depol_cal_sigMax_355": [1000, 1000, 1000, 1000], + "depol_cal_SNRmin_355": [1, 1, 1, 1], + "depol_cal_sigMax_355": [2000, 2000, 2000, 2000], "rel_std_dplus_355": 0.2, "rel_std_dminus_355": 0.2, "depol_cal_segmentLen_355": 40, @@ -81,73 +110,57 @@ "rel_std_dminus_1064": 0.2, "depol_cal_segmentLen_1064": 40, "depol_cal_smoothWin_1064": 8, - "isFR": [1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0], - "isNR": [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0], - "is532nm": [0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0], - "isRR": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - "is355nm": [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], - "is1064nm": [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], - "isTot": [1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0], - "isCross": [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], - "isParallel": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - "is387nm": [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], - "is407nm": [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], - "is607nm": [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0], - "is1058nm": [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0], + "isFR": [1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0], + "isNR": [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "is532nm": [0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "is355nm": [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "is1064nm": [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0], + "isTot": [1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "isCross": [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + "isRR": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + "isParallel": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "is387nm": [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "is407nm": [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "is607nm": [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "is1058nm": [], "minSNR_4_sigNorm": [], + "heightFullOverlap": [900, 900, 900, 900, 800, 800, 800, 800, 150, 150, 150, 150, 150, 800, 800,800, 900, 900, 900, 900, 800], + "channelTags": [], - "channelTag": ["FR-total-355 nm", "FR-cross-355 nm", "FR-387 nm", "FR-407 nm", "FR-total-532 nm", "FR-cross-532 nm", "FR-607 nm", "FR-total-1064 nm", "NR-total-532 nm", "NR-607 nm", "NR-total-355 nm", "NR-387 nm", "unknown"], + "channelTag": ["FR-total-355 nm", "FR-cross-355 nm", "FR-387 nm", "FR-407 nm", "FR-total-532 nm", "FR-cross-532 nm", "FR-607 nm", "FR-total-1064 nm", "NR-total-532 nm", "NR-607 nm", "NR-total-355 nm", "NR-387 nm", "DFOV","1058","1064s","none", "none", "none", "none", "none", "none"], + "minPC_fog": 60, - "TR": [0.898, 1086, 1, 1, 1.45, 778.8, 1, 1, 1, 1, 1, 1, 1], - "G": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], - "H": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], - "K": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], - "voldepol_error_355": [-999, -999, -999 ], - "voldepol_error_532": [-999, -999, -999 ], - "voldepol_error_1064": [-999, -999, -999], - "overlapCalMode": 1, - "overlapCorMode": 3, - "overlapSmoothBins": 8, + "saturate_thresh": 100, - "heightFullOverlap": [500, 500, 500, 500, 500, 500, 500, 500, 150, 150, 150, 150, 150], - "minSNR_4_sigNorm": [10], - "cloudScreenMode": 1, - "maxSigSlope4FilterCloud": 3e6, - "maxSigSlope4FilterCloud_NR": 4e5, - "intNProfiles": 120, - "minIntNProfiles": 90, - "meteorDataSource": "gdas1", - "flagUseLatestGDAS": true, - "radiosondeType": 1, - "radiosondeSitenum": 14430, - "minDecomLogDist355": 0.2, - "minDecomLogDist532": 0.2, - "minDecomLogDist1064": 0.2, - "maxDecomHeight355": 10000, - "maxDecomHeight532": 10000, - "maxDecomHeight1064": 10000, - "maxDecomThickness355": 1500, - "maxDecomThickness532": 1500, - "maxDecomThickness1064": 1500, - "decomSmoothWin355": 40, + + "minDecomLogDist355": 2, + "minDecomLogDist532": 2, + "minDecomLogDist1064": 2, + "maxDecomHeight355": 18000, + "maxDecomHeight532": 18000, + "maxDecomHeight1064": 25000, + "maxDecomThickness355": 6000, + "maxDecomThickness532": 3000, + "maxDecomThickness1064": 2000, + "decomSmoothWin355": 80, "decomSmoothWin532": 40, - "decomSmoothWin1064": 120, - "minRefThickness355": 700, - "minRefThickness532": 700, - "minRefThickness1064": 700, - "minRefDeltaExt355": 1, - "minRefDeltaExt532": 1, - "minRefDeltaExt1064": 1, + "decomSmoothWin1064": 30, + "minRefThickness355": 1000, + "minRefThickness532": 600, + "minRefThickness1064": 300, + "minRefDeltaExt355": 11, + "minRefDeltaExt532": 11, + "minRefDeltaExt1064":11, "refH_FR_355": [], "refH_FR_532": [], "refH_FR_1064": [], "refH_NR_355": [2500, 3000], "refH_NR_532": [2500, 3000], - "minRefSNR355": 10, + "minRefSNR355": 5, "minRefSNR532": 5, - "minRefSNR1064": 3, - "minRefSNR_NR_355": 10, - "minRefSNR_NR_532": 5, + "minRefSNR1064": 0.5, + "minRefSNR_NR_355": 3, + "minRefSNR_NR_532": 3, "LR355": 50, "LR532": 50, "LR1064": 50, @@ -170,30 +183,28 @@ "minLRConstrainFernald": 1, "maxLRConstrainFernald": 150, "minDeltaAOD": 0.01, - "minRamanRefSNR355": 50, - "minRamanRefSNR532": 20, - "minRamanRefSNR1064": 10, - "minRamanRefSNR387": 40, - "minRamanRefSNR607": 20, - "minRamanRefSNR_NR_355": 50, - "minRamanRefSNR_NR_532": 20, - "minRamanRefSNR_NR_387": 40, - "minRamanRefSNR_NR_607": 20, + "minRamanRefSNR355": 5, + "minRamanRefSNR532": 4, + "minRamanRefSNR1064": 0.5, + "minRamanRefSNR387": 3, + "minRamanRefSNR607": 2, + "minRamanRefSNR_NR_355": 2, + "minRamanRefSNR_NR_532": 4, + "minRamanRefSNR_NR_387": 2, + "minRamanRefSNR_NR_607": 2, "min_RR_RefSNR1058": 0, - "angstrexp": 0.9, - "angstrexp_NR": 0.9, + "LCMeanWindow": 50, - "LCMeanMinIndx": 70, + "LCMeanMinIndx": 150, "LCMeanMaxIndx": 1000, "LCCalibrationStatus": ["none", "Klett", "Raman", "Defaults", "History"], - "quasi_smooth_h": [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8], - "quasi_smooth_t": [10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10], - "IWV_instrument": "AERONET", - "maxIWVTLag": 0.0833, + "quasi_smooth_h": [10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10], + "quasi_smooth_t": [10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10], + "tTwilight": 0.0347, "hWVCaliBase": 120, "hWVCaliTop": 8000, - "minSNRWVCali": 10, + "minSNRWVCali": 1, "clear_thres_par_beta_1064": 1e-8, "turbid_thres_par_beta_1064": 2e-7, @@ -221,7 +232,7 @@ "xLim_Profi_WVMR": [0, 50], "xLim_Profi_RCS": [0.3, 1000], "xLim_Profi_LR": [0, 120], - "xLim_Profi_AE": [-1,2], + "xLim_Profi_AE": [-1,3], "xLim_beta_532_Poliphon": [0, 10], "yLim_beta_532_Poliphon": [0, 8000], "yLim_LC_355": [0, 1e14], @@ -231,6 +242,9 @@ "yLim_LC_607": [0, 3e14], "yLim_LC_355_NR": [0, 1e13], "yLim_LC_532_NR": [0, 2e13], + "yLim_LC_ratio_355_387": [0, 2], + "yLim_LC_ratio_532_607": [0, 2], + "yLim_WVConst": [0, 40], "yLim_FR_RCS": [0, 20000], "yLim_NR_RCS": [0, 3000], @@ -245,21 +259,24 @@ "yLim_Profi_DR": [0, 20000], "yLim_Profi_Bsc": [0, 20000], "yLim_Profi_WV_RH": [0, 7000], - "yLim_LC_ratio_355_387": [0, 2], - "yLim_LC_ratio_532_607": [0, 2], + "yLim_all_profiles_high_range": [0, 18000], + "yLim_all_profiles_low_range": [0, 5000], + + "yLim_depolConst_355": [0, 50], "yLim_depolConst_532": [0, 50], - "yLim_depolConst_1064": [0, 50], + "yLim_depolConst_1064": [0, 10], "yLim_cloudinfo": [0, 2000], - "yLim_all_profiles_high_range": [0, 18000], - "yLim_all_profiles_low_range": [0, 5000], + "yLim_overlap": [0, 3000], + "zLim_att_beta_355": [0, 15], "zLim_att_beta_532": [0, 5], "zLim_att_beta_1064": [0, 2], - "zLim_quasi_beta_355": [0, 10], - "zLim_quasi_beta_532": [0, 6], - "zLim_quasi_beta_1064": [0, 3], + "zLim_quasi_beta_355": [0, 10e-6], + "zLim_quasi_beta_532": [0, 6e-6], + "zLim_quasi_beta_1064": [0, 3e-6], "zLim_quasi_Par_DR_532": [0, 0.4], + "zLim_quasi_ANG": [-1,3], "zLim_FR_RCS_355": [1e-2, 30], "zLim_FR_RCS_387": [1e-2, 30], "zLim_FR_RCS_407": [1e-2, 35], @@ -290,12 +307,12 @@ "Data_Originator_email": "", "comment": "", + "partnerLabel": "", + "calibrationDB": "polly_calibration.db", "logbookFile": "", "logbookPath": "", "logbookFileName": "", "radiosondeFolder": "", - "calibrationDB": "polly_calibration.db", "imgFormat": "png", - "partnerLabel": "", "prodSaveList": ["overlap", "aerProfFR", "aerProfNR", "aerProfOC", "aerAttBetaFR", "aerAttBetaOC", "aerAttBetaNR", "WVMR_RH", "volDepol", "quasiV1", "quasiV2", "TC", "TCV2", "cloudinfo", "poliphon_one","RCS"] } diff --git a/lib/interface/concat_pollyxt_lvl0.py b/lib/interface/concat_pollyxt_lvl0.py index a23ceb77..20c54bc4 100644 --- a/lib/interface/concat_pollyxt_lvl0.py +++ b/lib/interface/concat_pollyxt_lvl0.py @@ -78,7 +78,8 @@ def main(): merging_flag = False merging_flag = concat_files() - get_pollyxt_logbook_files() + if device != "martha": + get_pollyxt_logbook_files() if merging_flag == True: sys.exit(0) # Indicates success/True @@ -94,7 +95,129 @@ def get_input_path(timestamp,device,raw_folder): print(input_path) return input_path +def get_input_path_martha(timestamp,device,raw_folder): + #raw_folder="/data/level0/martha/MARTHA_DATA" + YYYY=timestamp[0:4] + MM=timestamp[4:6] + #input_path = Path(raw_folder,device,"data_zip",f"{YYYY}{MM}") + input_path = Path(raw_folder) + print(input_path) + return input_path + ### start of function concat_pollyxt_files +def get_martha_files(): + ''' + This function locates multiple pollyxt level0 nc-zip files from one day measurements, + unzipps the files to output_path + and returns a list of files to be merged + and the title of the new merged nc-file + ''' + input_path = get_input_path_martha(timestamp,device,raw_folder) + path_exist = Path(input_path) + + if path_exist.exists() == True: + + ## set the searchpattern for the zipped-nc-files: + YYYY=timestamp[0:4] + MM=timestamp[4:6] + DD=timestamp[6:8] + + zip_searchpattern = str(YYYY)+'_'+str(MM)+'_'+str(DD)+'*_*[0-9].nc' + + polly_files = Path(r'{}'.format(input_path)).glob('{}'.format(zip_searchpattern)) + polly_zip_files_list0 = [x for x in polly_files if x.is_file()] + + + + ## convert type path to type string + polly_zip_files_list = [] + for file in polly_zip_files_list0: + polly_zip_files_list.append(str(file)) + + if len(polly_zip_files_list) < 1: + print('no files found!') + sys.exit() + else: + print(polly_zip_files_list) + + # Ensure the destination directory exists + Path(output_path).mkdir(parents=True, exist_ok=True) + + polly_files_list = [] + to_unzip_list = [] + for zip_file in polly_zip_files_list: + ## check for size of zip-files to ensure to exclude bad measurement files with wrong timestamp e.g. 19700101 + f_size = os.path.getsize(zip_file) + print(zip_file) + if f_size > 150000: + print(f_size) + print("filesize passes") + else: + print(f_size) + print("filesize too small, file will be skipped!") + continue ## go to next file + +# ## check if zipfile is a valid zip-file +# if not is_zipfile(zip_file): +# print(f"invalid zip-file: {zip_file}\nskipping file.") +# #polly_zip_files_list.remove(zip_file) +# continue +# else: +# pass + + #unzipped_nc = Path(zip_file).name + #unzipped_nc = Path(unzipped_nc) + #unzipped_nc = Path(output_path,unzipped_nc) + #print('here....') + #polly_files_list.append(unzipped_nc) + path = Path(zip_file) + polly_files_list.append(path) + + ### check if unzipped files already exists in outputfolder + #if path.is_file() == False: + # to_unzip_list.append(zip_file) + #if path.is_file() == True: + # os.remove(unzipped_nc) + # to_unzip_list.append(zip_file) + +# ## unzipping +# date_pattern = str(YYYY)+'_'+str(MM)+'_'+str(DD) +# if len(to_unzip_list) > 0: +# ## if working remotly on windows, copy zipped files first, than unzip +# if os_name.lower() == 'windows': +# print("\nCopy zipped files to local drive...") +# for zip_file in to_unzip_list: +# print(zip_file) +# shutil.copy2(Path(zip_file), Path(output_path) / Path(zip_file).name) +# print("\nUnzipping...") +# for zip_file in Path(output_path).iterdir(): +# if zip_file.is_file() and date_pattern in zip_file.stem and zip_file.suffix == '.zip': +# with ZipFile(zip_file, 'r') as zip_ref: +# print("unzipping "+str(zip_file)) +# zip_ref.extractall(output_path) +# print("Removing .zip file...") +# os.remove(zip_file) +# +# else: +# print("\nUnzipping...") +# for zip_file in to_unzip_list: +# with ZipFile(zip_file, 'r') as zip_ref: +# print("unzipping "+zip_file) +# zip_ref.extractall(output_path) + + + ## sort lists + polly_files_list.sort() + + print("\n"+str(len(polly_files_list))+" files found:\n") + print(polly_files_list) + print("\n") + + else: + print("\nNo data was found in {}. Correct path?\n".format(input_path)) + sys.exit() + return polly_files_list + def get_pollyxt_files(): ''' This function locates multiple pollyxt level0 nc-zip files from one day measurements, @@ -315,7 +438,10 @@ def checking_vars(): 'zenithangle' ] - polly_files_list = get_pollyxt_files() + if device == "martha": + polly_files_list = get_martha_files() + else: + polly_files_list = get_pollyxt_files() if len(polly_files_list) == 1: return polly_files_list @@ -591,7 +717,10 @@ def concat_files(): if len(sel_polly_files_list) == 1: print("\nOnly one file found. Nothing to merge!\n") - os.rename(sel_polly_files_list[0],Path(output_path,filestring)) + if device == "martha": + shutil.copy2(sel_polly_files_list[0],Path(output_path,filestring)) + else: + os.rename(sel_polly_files_list[0],Path(output_path,filestring)) return True else: # sel_polly_files_list = [ str(el) for el in sel_polly_files_list] @@ -625,10 +754,11 @@ def write_netcdf(ds: xarray.Dataset, out_file: Path) -> None: ds.close() - print("\ndeleting individual .nc files ...") - for el in sel_polly_files_list: - print(el) - os.remove(el) + if device != "martha": + print("\ndeleting individual .nc files ...") + for el in sel_polly_files_list: + print(el) + os.remove(el) destination_file = Path(output_path,filestring) if os.path.exists(destination_file): os.remove(destination_file) # Remove the existing destination file diff --git a/lib/interface/picassoProcV3.m b/lib/interface/picassoProcV3.m index 8d2a6104..2c0912bf 100644 --- a/lib/interface/picassoProcV3.m +++ b/lib/interface/picassoProcV3.m @@ -50,9 +50,14 @@ % processing report. % % HISTORY: -% - 2021-06-25: first edition by Zhenping +% - 2021-06-25: first edition by Zhenping Yin +% - 2023-06-06: Overlap function using Raman method was added by Cristofer Jimenez +% - 2023-06-14: POLIPHON method (step 1) added by Athena Floutsi +% - 2024-08-28: GHK formalism for depol calculation implemented by Moritz Haarig +% - 2025-02-19: Smoothing added into meteo profiles, read all meteo data for HR products and interpolate, recalculation of signals using mean shot_number by Cristofer Jimenez +% - 2025-03-14: Compute and save attenuated backscatter co and cross polarized by Cristofer Jimenez % -% .. Authors: - zhenping@tropos.de +% .. Authors: - zhenping@tropos.de, jimenez@tropos.de, floutsi@tropos.de, haarig@tropos.de global PicassoConfig global CampaignConfig @@ -508,6 +513,7 @@ end %% Polarization calibration +print_msg('Polarization calibration. \n', 'flagTimestamp', true); data.polCaliFac355=NaN; data.polCaliFac532=NaN; data.polCaliFac1064=NaN; @@ -535,6 +541,10 @@ %Taking the eta with lowest standard deviation [~, index_min] = min(data.polCali355Attri.polCaliEtaStd); data.polCaliEta355=data.polCali355Attri.polCaliEta(index_min); + print_msg('Depol cali 355 etas: \n'); + disp(data.polCali355Attri.polCaliEta); + print_msg('Depol Cali eta used355:\n'); + disp(data.polCaliEta355); else warning('Cross or total channel at 355 nm does not exist.'); data.polCaliEta355=NaN; @@ -565,6 +575,10 @@ %Taking the eta with lowest standard deviation [~, index_min] = min(data.polCali532Attri.polCaliEtaStd); data.polCaliEta532=data.polCali532Attri.polCaliEta(index_min); + print_msg('Depol cali 532 etas: \n'); + disp(data.polCali532Attri.polCaliEta); + print_msg('Depol Cali eta used532:\n'); + disp(data.polCaliEta532); else warning('Cross or total channel at 532 nm does not exist.'); data.polCaliEta532=NaN; @@ -595,6 +609,10 @@ %Taking the eta with lowest standard deviation [~, index_min] = min(data.polCali1064Attri.polCaliEtaStd); data.polCaliEta1064=data.polCali1064Attri.polCaliEta(index_min); + print_msg('Depol cali 1064 etas: \n'); + disp(data.polCali1064Attri.polCaliEta); + print_msg('Depol Cali eta used1064:\n'); + disp(data.polCaliEta1064); else warning('Cross or total channel at 1064 nm does not exist.') data.polCaliEta1064=NaN; @@ -756,14 +774,14 @@ %% Calculate molecular scattering properties for iGrp = 1:size(clFreGrps, 1) - [mBsc355(iGrp,:), mExt355(iGrp,:)] = rayleigh_scattering(355, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [mBsc387(iGrp,:), mExt387(iGrp,:)] = rayleigh_scattering(387, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [mBsc407(iGrp,:), mExt407(iGrp,:)] = rayleigh_scattering(407, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [mBsc532(iGrp,:), mExt532(iGrp,:)] = rayleigh_scattering(532, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [mBsc607(iGrp,:), mExt607(iGrp,:)] = rayleigh_scattering(607, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [mBsc1058(iGrp,:), mExt1058(iGrp,:)] = rayleigh_scattering(1058, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [mBsc1064(iGrp,:), mExt1064(iGrp,:)] = rayleigh_scattering(1064, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - number_density(iGrp,:) = number_density_at_pt(data.pressure(iGrp, :), data.temperature(iGrp, :)+ 273.17, 70, true); + [mBsc355(iGrp,:), mExt355(iGrp,:)] = rayleigh_scattering(355, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [mBsc387(iGrp,:), mExt387(iGrp,:)] = rayleigh_scattering(387, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [mBsc407(iGrp,:), mExt407(iGrp,:)] = rayleigh_scattering(407, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [mBsc532(iGrp,:), mExt532(iGrp,:)] = rayleigh_scattering(532, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [mBsc607(iGrp,:), mExt607(iGrp,:)] = rayleigh_scattering(607, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [mBsc1058(iGrp,:), mExt1058(iGrp,:)] = rayleigh_scattering(1058, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [mBsc1064(iGrp,:), mExt1064(iGrp,:)] = rayleigh_scattering(1064, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + number_density(iGrp,:) = number_density_at_pt(data.pressure(iGrp, :), data.temperature(iGrp, :)+ 273.15, 70, true); data.molBsc355(iGrp,:)= mBsc355(iGrp,:); data.molBsc532(iGrp,:)= mBsc532(iGrp,:); data.molBsc1064(iGrp,:)= mBsc1064(iGrp,:); @@ -3416,7 +3434,7 @@ trans387 = exp(-cumsum(mExt387(iGrp,:) .* [data.distance0(1), diff(data.distance0)])); trans407 = exp(-cumsum(mExt407(iGrp,:) .* [data.distance0(1), diff(data.distance0)])); - rhoAir = rho_air(data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17); + rhoAir = rho_air(data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15); [thisWVconst, thisWVconstStd, thisWVAttri] = pollyWVCali(data.height, ... sig387, bg387, sig407, E_tot_1064_IWV, E_tot_1064_cali, E_tot_1064_cali_std, ... @@ -3484,7 +3502,7 @@ % calculate saturated water vapor pressure es = saturated_vapor_pres(data.temperature(iGrp, :)); - rhoAir = rho_air(data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17); + rhoAir = rho_air(data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15); % calculate wvmr and rh data.wvmr(iGrp, :) = sig407 ./ sig387 .* trans387 ./ trans407 .* data.wvconstUsed; @@ -3528,6 +3546,11 @@ data.quality_mask_RH = 3 * ones(size(data.signal, 2), size(data.signal, 3)); ones_WV= ones(size(data.signal, 2), size(data.signal, 3)); +TimeM=floor(data.mTime(1)*24/3)*3/24-3/24:3/24:ceil(data.mTime(end)*24/3)*3/24+3/24; %timegrid to search gdas files (3h step) + +[TimeMg, HeightMg] = meshgrid(data.height, TimeM); %mehsgrids for 2d interpolation +[mTimeg, Heightg] = meshgrid(data.height, data.mTime); + if (sum(flag387FR) == 1) && (sum(flag407 == 1)) sig387 = squeeze(data.signal(flag387FR, :, :)); @@ -3556,7 +3579,7 @@ % read the meteorological data [temp, pres, ~, ~] = loadMeteor(... - mean(data.mTime), data.alt, ... + TimeM, data.alt, ... 'meteorDataSource', PollyConfig.meteorDataSource, ... 'gdas1Site', PollyConfig.gdas1Site, ... 'meteo_folder', PollyConfig.meteo_folder, ... @@ -3565,21 +3588,21 @@ 'radiosondeType', PollyConfig.radiosondeType); % repmat the array to matrix as the size of data.signal - temperature = repmat(transpose(temp), 1, length(data.mTime)); - pressure = repmat(transpose(pres), 1, length(data.mTime)); + temperature = transpose(interp2(TimeMg, HeightMg, temp, mTimeg, Heightg, 'linear')); + pressure = transpose(interp2(TimeMg, HeightMg, pres, mTimeg, Heightg, 'linear')); % calculate the molecule optical properties - [~, mExt387_highres] = rayleigh_scattering(387, transpose(pressure(:, 1)), transpose(temperature(:, 1)) + 273.17, 380, 70); - [~, mExt407_highres] = rayleigh_scattering(407, transpose(pressure(:, 1)), transpose(temperature(:, 1)) + 273.17, 380, 70); + [~, mExt387_highres] = rayleigh_scattering(387, pres, temp + 273.15, 380, 70); + [~, mExt407_highres] = rayleigh_scattering(407, pres, temp + 273.15, 380, 70); trans387 = exp(- cumsum(mExt387_highres .* [data.distance0(1), diff(data.distance0)])); trans407 = exp(- cumsum(mExt407_highres .* [data.distance0(1), diff(data.distance0)])); - TRANS387 = repmat(transpose(trans387), 1, length(data.mTime)); - TRANS407 = repmat(transpose(trans407), 1, length(data.mTime)); + TRANS387 = transpose(interp2(TimeMg, HeightMg, trans387, mTimeg, Heightg, 'linear')); + TRANS407 = transpose(interp2(TimeMg, HeightMg, trans407, mTimeg, Heightg, 'linear')); % calculate the saturation water vapor pressure - es = saturated_vapor_pres(temperature(:, 1)); - ES = repmat(es, 1, length(data.mTime)); + ES = saturated_vapor_pres(temperature); + - % rhoAir = rho_air(pressure(:, 1), temperature(:, 1) + 273.17); + % rhoAir = rho_air(pressure(:, 1), temperature(:, 1) + 273.15); % RHOAIR = repmat(rhoAir, 1, length(data.mTime)); % DIFFHeight = repmat(transpose([data.height(1), diff(data.height)]), 1, length(data.mTime)); @@ -4170,75 +4193,75 @@ %% attenuated backscatter print_msg('Start calculating attenuated backscatter.\n', 'flagTimestamp', true); -data.att_beta_355 = NaN(length(data.height), length(data.mTime)); +data.att_beta_355 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag355t) == 1) - data.att_beta_355 = squeeze(data.signal(flag355t, :, :)) .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed355; + data.att_beta_355 = squeeze(data.signal(flag355t, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed355; data.att_beta_355(:, data.depCalMask) = NaN; end -data.att_beta_532 = NaN(length(data.height), length(data.mTime)); +data.att_beta_532 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag532t) == 1) - data.att_beta_532 = squeeze(data.signal(flag532t, :, :)) .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed532; + data.att_beta_532 = squeeze(data.signal(flag532t, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed532; data.att_beta_532(:, data.depCalMask) = NaN; end -data.att_beta_1064 = NaN(length(data.height), length(data.mTime)); +data.att_beta_1064 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag1064t) == 1) - data.att_beta_1064 = squeeze(data.signal(flag1064t, :, :)) .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed1064; + data.att_beta_1064 = squeeze(data.signal(flag1064t, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed1064; data.att_beta_1064(:, data.depCalMask) = NaN; end -att_beta_387 = NaN(length(data.height), length(data.mTime)); +att_beta_387 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag387FR) == 1) - att_beta_387 = squeeze(data.signal(flag387FR, :, :)) .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed387; + att_beta_387 = squeeze(data.signal(flag387FR, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed387; att_beta_387(:, data.depCalMask) = NaN; end -att_beta_607 = NaN(length(data.height), length(data.mTime)); +att_beta_607 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag607FR) == 1) - att_beta_607 = squeeze(data.signal(flag607FR, :, :)) .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed607; + att_beta_607 = squeeze(data.signal(flag607FR, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed607; att_beta_607(:, data.depCalMask) = NaN; end -data.att_beta_OC_355 = NaN(length(data.height), length(data.mTime)); +data.att_beta_OC_355 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag355t) == 1) - data.att_beta_OC_355 = data.sigOLCor355 .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed355; + data.att_beta_OC_355 = data.sigOLCor355 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed355; data.att_beta_OC_355(:, data.depCalMask) = NaN; end -data.att_beta_OC_532 = NaN(length(data.height), length(data.mTime)); +data.att_beta_OC_532 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag532t) == 1) - data.att_beta_OC_532 = data.sigOLCor532 .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed532; + data.att_beta_OC_532 = data.sigOLCor532 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed532; data.att_beta_OC_532(:, data.depCalMask) = NaN; end -data.att_beta_OC_1064 = NaN(length(data.height), length(data.mTime)); +data.att_beta_OC_1064 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag1064t) == 1) - data.att_beta_OC_1064 = data.sigOLCor1064 .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed1064; + data.att_beta_OC_1064 = data.sigOLCor1064 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed1064; data.att_beta_OC_1064(:, data.depCalMask) = NaN; end -% att_beta_OC_387 = NaN(length(data.height), length(data.mTime)); +% att_beta_OC_387 = NaN(length(data.distance0), length(data.mTime)); % if (sum(flag387FR) == 1) -% att_beta_OC_387 = sigOLCor387 .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed387; +% att_beta_OC_387 = sigOLCor387 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed387; % att_beta_OC_387(:, data.depCalMask) = NaN; % end -% att_beta_OC_607 = NaN(length(data.height), length(data.mTime)); +% att_beta_OC_607 = NaN(length(data.distance0), length(data.mTime)); % if (sum(flag607FR) == 1) -% att_beta_OC_607 = sigOLCor607 .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed607; +% att_beta_OC_607 = sigOLCor607 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed607; % att_beta_OC_607(:, data.depCalMask) = NaN; % end -data.att_beta_NR_355 = NaN(length(data.height), length(data.mTime)); +data.att_beta_NR_355 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag355NR) == 1) - data.att_beta_NR_355 = squeeze(data.signal(flag355NR, :, :)) .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed355NR; + data.att_beta_NR_355 = squeeze(data.signal(flag355NR, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed355NR; data.att_beta_NR_355(:, data.depCalMask) = NaN; end -data.att_beta_NR_532 = NaN(length(data.height), length(data.mTime)); +data.att_beta_NR_532 = NaN(length(data.distance0), length(data.mTime)); if (sum(flag532NR) == 1) - data.att_beta_NR_532 = squeeze(data.signal(flag532NR, :, :)) .* repmat(transpose(data.height), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed532NR; + data.att_beta_NR_532 = squeeze(data.signal(flag532NR, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed532NR; data.att_beta_NR_532(:, data.depCalMask) = NaN; end @@ -4311,11 +4334,32 @@ end print_msg('Finish.\n', 'flagTimestamp', true); +%% Co (para) and cross (perp) polarized components in attenuated backscatter +print_msg('Start calculating co and cross attenuated backscatter.\n', 'flagTimestamp', true); + +data.att_beta_para_355 = NaN(length(data.height), length(data.mTime)); +data.att_beta_perp_355 = NaN(length(data.height), length(data.mTime)); +if (sum(flag355t) == 1) && (sum(flag355c) == 1) +data.att_beta_para_355=data.att_beta_355./(1+PollyConfig.TR(flag355t)*data.vdr355); +data.att_beta_perp_355=data.att_beta_para_355.*data.vdr355; +end +data.att_beta_para_532 = NaN(length(data.height), length(data.mTime)); +data.att_beta_perp_532 = NaN(length(data.height), length(data.mTime)); +if (sum(flag532t) == 1) && (sum(flag532c) == 1) +data.att_beta_para_532=data.att_beta_532./(1+PollyConfig.TR(flag532t)*data.vdr532); +data.att_beta_perp_532=data.att_beta_para_532.*data.vdr532; +end +data.att_beta_para_1064 = NaN(length(data.height), length(data.mTime)); +data.att_beta_perp_1064 = NaN(length(data.height), length(data.mTime)); +if (sum(flag1064t) == 1) && (sum(flag1064c) == 1) +data.att_beta_para_1064=data.att_beta_1064./(1+PollyConfig.TR(flag1064t)*data.vdr1064); +data.att_beta_perp_1064=data.att_beta_para_1064.*data.vdr1064; +end %% Quasi-retrieval (V1) print_msg('Start quasi-retrieval (V1).\n', 'flagTimestamp', true); % load meteorological data -[temperature, pressure, ~, ~, ~, thisMeteorAttri] = loadMeteor(mean(data.mTime), data.alt, ... +[temperature, pressure, ~, ~, ~, thisMeteorAttri] = loadMeteor(TimeM, data.alt, ... 'meteorDataSource', PollyConfig.meteorDataSource, ... 'gdas1Site', PollyConfig.gdas1Site, ... 'meteo_folder', PollyConfig.meteo_folder, ... @@ -4340,9 +4384,9 @@ % Rayleigh scattering %---------------achtung - [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.17, 380, 70); - mBsc355 = repmat(transpose(mBsc355), 1, length(data.mTime)); - mExt355 = repmat(transpose(mExt355), 1, length(data.mTime)); + [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.15, 380, 70); + mBsc355 = transpose(interp2(TimeMg, HeightMg, mBsc355, mTimeg, Heightg, 'linear')); + mExt355 = transpose(interp2(TimeMg, HeightMg, mExt355, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4365,10 +4409,10 @@ att_beta_532_qsi = smooth2(att_beta_532_qsi, PollyConfig.quasi_smooth_h(flag532t), PollyConfig.quasi_smooth_t(flag532t)); % Rayleigh scattering - [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.17, 380, 70); + [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70); %achtung - mBsc532 = repmat(transpose(mBsc532), 1, length(data.mTime)); - mExt532 = repmat(transpose(mExt532), 1, length(data.mTime)); + mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear')); + mExt532 = transpose(interp2(TimeMg, HeightMg, mExt532, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4392,9 +4436,9 @@ % Rayleigh scattering %achtung - [mBsc1064, mExt1064] = rayleigh_scattering(1064, pressure, temperature + 273.17, 380, 70); - mBsc1064 = repmat(transpose(mBsc1064), 1, length(data.mTime)); - mExt1064 = repmat(transpose(mExt1064), 1, length(data.mTime)); + [mBsc1064, mExt1064] = rayleigh_scattering(1064, pressure, temperature + 273.15, 380, 70); + mBsc1064 = transpose(interp2(TimeMg, HeightMg, mBsc1064, mTimeg, Heightg, 'linear')); + mExt1064 = transpose(interp2(TimeMg, HeightMg, mExt1064, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4421,8 +4465,8 @@ sig532CSm = smooth2(sig532C, PollyConfig.quasi_smooth_h(flag532c), PollyConfig.quasi_smooth_t(flag532c)); % Rayleigh scattering - [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.17, 380, 70); - mBsc532 = repmat(transpose(mBsc532), 1, length(data.mTime)); + [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70); + mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4431,7 +4475,8 @@ PollyConfig.G(flag532t),PollyConfig.G(flag532c), ... PollyConfig.H(flag532t),PollyConfig.H(flag532c), ... data.polCaliEta532); - data.qsiPDR532V1 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V1 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + %data.qsiPDR532V1 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V1 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + data.qsiPDR532V1 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) ./ (data.qsiBsc532V1 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; data.qsiPDR532V1((data.quality_mask_vdr_532 ~= 0) | (data.quality_mask_532 ~= 0)) = NaN; end else @@ -4445,14 +4490,15 @@ sig532CSm = smooth2(sig532C, PollyConfig.quasi_smooth_h(flag532c), PollyConfig.quasi_smooth_t(flag532c)); % Rayleigh scattering - [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.17, 380, 70); - mBsc532 = repmat(transpose(mBsc532), 1, length(data.mTime)); + [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70); + mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; vdr532Sm = pollyVDR2(sig532TSm, sig532CSm, PollyConfig.TR(flag532t), PollyConfig.TR(flag532c), data.polCaliFac532); - data.qsiPDR532V1 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V1 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + %data.qsiPDR532V1 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V1 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + data.qsiPDR532V1 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) ./ (data.qsiBsc532V1 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; data.qsiPDR532V1((data.quality_mask_vdr_532 ~= 0) | (data.quality_mask_532 ~= 0)) = NaN; end end @@ -4527,11 +4573,11 @@ att_beta_387_qsi = smooth2(att_beta_387_qsi, PollyConfig.quasi_smooth_h(flag387FR), PollyConfig.quasi_smooth_t(flag387FR)); % Rayleigh scattering - [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.17, 380, 70); - [~, mExt387] = rayleigh_scattering(387, pressure, temperature + 273.17, 380, 70); - mBsc355 = repmat(transpose(mBsc355), 1, length(data.mTime)); - mExt355 = repmat(transpose(mExt355), 1, length(data.mTime)); - mExt387 = repmat(transpose(mExt387), 1, length(data.mTime)); + [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.15, 380, 70); + [~, mExt387] = rayleigh_scattering(387, pressure, temperature + 273.15, 380, 70); + mBsc355 = transpose(interp2(TimeMg, HeightMg, mBsc355, mTimeg, Heightg, 'linear')); + mExt355 = transpose(interp2(TimeMg, HeightMg, mExt355, mTimeg, Heightg, 'linear')); + mExt387 = transpose(interp2(TimeMg, HeightMg, mExt387, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4551,11 +4597,11 @@ att_beta_607_qsi = smooth2(att_beta_607_qsi, PollyConfig.quasi_smooth_h(flag607FR), PollyConfig.quasi_smooth_t(flag607FR)); % Rayleigh scattering - [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.17, 380, 70); - [~, mExt607] = rayleigh_scattering(607, pressure, temperature + 273.17, 380, 70); - mBsc532 = repmat(transpose(mBsc532), 1, length(data.mTime)); - mExt532 = repmat(transpose(mExt532), 1, length(data.mTime)); - mExt607 = repmat(transpose(mExt607), 1, length(data.mTime)); + [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70); + [~, mExt607] = rayleigh_scattering(607, pressure, temperature + 273.15, 380, 70); + mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear')); + mExt532 = transpose(interp2(TimeMg, HeightMg, mExt532, mTimeg, Heightg, 'linear')); + mExt607 = transpose(interp2(TimeMg, HeightMg, mExt607, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4575,11 +4621,11 @@ att_beta_607_qsi = smooth2(att_beta_607_qsi, PollyConfig.quasi_smooth_h(flag607FR), PollyConfig.quasi_smooth_t(flag607FR)); % Rayleigh scattering - [mBsc1064, mExt1064] = rayleigh_scattering(1064, pressure, temperature + 273.17, 380, 70); - [~, mExt607] = rayleigh_scattering(607, pressure, temperature + 273.17, 380, 70); - mBsc1064 = repmat(transpose(mBsc1064), 1, length(data.mTime)); - mExt1064 = repmat(transpose(mExt1064), 1, length(data.mTime)); - mExt607 = repmat(transpose(mExt607), 1, length(data.mTime)); + [mBsc1064, mExt1064] = rayleigh_scattering(1064, pressure, temperature + 273.15, 380, 70); + [~, mExt607] = rayleigh_scattering(607, pressure, temperature + 273.15, 380, 70); + mBsc1064 = transpose(interp2(TimeMg, HeightMg, mBsc1064, mTimeg, Heightg, 'linear')); + mExt1064 = transpose(interp2(TimeMg, HeightMg, mExt1064, mTimeg, Heightg, 'linear')); + mExt607 = transpose(interp2(TimeMg, HeightMg, mExt607, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4601,8 +4647,8 @@ sig532CSm = smooth2(sig532C, PollyConfig.quasi_smooth_h(flag532c), PollyConfig.quasi_smooth_t(flag532c)); % Rayleigh scattering - [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.17, 380, 70); - mBsc532 = repmat(transpose(mBsc532), 1, length(data.mTime)); + [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70); + mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; @@ -4611,7 +4657,8 @@ PollyConfig.G(flag532t),PollyConfig.G(flag532c), ... PollyConfig.H(flag532t),PollyConfig.H(flag532c), ... data.polCaliEta532); - data.qsiPDR532V2 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V2 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + %data.qsiPDR532V2 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V2 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + data.qsiPDR532V2 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) ./ (data.qsiBsc532V2 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; data.qsiPDR532V2((data.quality_mask_vdr_532 ~= 0) | (data.quality_mask_532 ~= 0)) = NaN; end else @@ -4625,14 +4672,15 @@ sig532CSm = smooth2(sig532C, PollyConfig.quasi_smooth_h(flag532c), PollyConfig.quasi_smooth_t(flag532c)); % Rayleigh scattering - [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.17, 380, 70); - mBsc532 = repmat(transpose(mBsc532), 1, length(data.mTime)); + [mBsc532, ~] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70); + mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear')); data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1'); data.quasiAttri.meteorSource = thisMeteorAttri.dataSource; data.quasiAttri.timestamp = thisMeteorAttri.datetime; vdr532Sm = pollyVDR2(sig532TSm, sig532CSm, PollyConfig.TR(flag532t), PollyConfig.TR(flag532c), data.polCaliFac532); - data.qsiPDR532V2 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V2 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + %data.qsiPDR532V2 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) .* (data.qsiBsc532V2 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; + data.qsiPDR532V2 = (vdr532Sm + 1) ./ (mBsc532 .* (PollyDefaults.molDepol532 - vdr532Sm) ./ (data.qsiBsc532V2 .* (1 + PollyDefaults.molDepol532)) + 1) - 1; data.qsiPDR532V2((data.quality_mask_vdr_532 ~= 0) | (data.quality_mask_532 ~= 0)) = NaN; end end diff --git a/lib/interface/picasso_go24.sh b/lib/interface/picasso_go24.sh index 6f1790ae..f5d8782a 100755 --- a/lib/interface/picasso_go24.sh +++ b/lib/interface/picasso_go24.sh @@ -12,7 +12,9 @@ display_help() { echo " -d, --device specify device, e.g. pollyxt_lacros" echo " -c, --config_file specify Picasso configuration file, e.g.: ~/Pollynet_Processing_Chain/config/pollynet_processing_chain_config_rsd2_andi.json" # echo " -o, --output specify folder where to put merged nc-file, e.g.: ~/todo_filelist" + echo " -b, --base_folder specify folder where level0 files are located; default is set to polly-folder, can be switched to e.g. martha-folder" echo " --force_merging specify whether files will be merged independently if attributes have changed or not; default: false" + echo " --no_merging activate this flag to directly run picasso without merging to 24h-file" echo " --todolist write merged level0-file into todo-list (set to true or false); default: true" echo " --matlab specify location of matlab-executable; default is just set to: matlab" echo " --proc execute picasso-processing-chain (set to true or false); default: true" @@ -27,6 +29,7 @@ display_help() { ## initialize parameters FORCE_MERGING="false" MATLABEXEC="matlab" +BASE_FOLDER="/data/level0/polly" PICASSO_CONFIG_FILE="" PICASSO_DIR_interface="$( cd "$(dirname "$0")" ; pwd -P )" PICASSO_DIR="$(dirname "$(dirname "$PICASSO_DIR_interface")")" @@ -38,6 +41,8 @@ flagProc="true" flagFProc="false" flagDeleteMergedFiles="true" flagPlotonly="false" +flag_no_merging="false +" filename="" filesize="" @@ -75,6 +80,13 @@ while :; do shift 2 ;; + -b | --base_folder) + if [ $# -ne 0 ]; then + BASE_FOLDER="$2" + fi + shift 2 + ;; + --force_merging) if [ $# -ne 0 ]; then FORCE_MERGING="$2" @@ -110,6 +122,12 @@ while :; do flagPlotonly="true" shift 1 ;; + + --no_merging) + flag_no_merging="true" + shift 1 + ;; + -h | --help) display_help # Call your function @@ -178,6 +196,11 @@ main() { echo $DEVICE for DATE in ${DATE_LS[@]}; do echo $DATE + if [[ "$flag_no_merging" == "true" ]];then + echo "processing without merging" + process_without_merging $DEVICE $DATE + exit 1 + fi merging $DEVICE $DATE ## merging of level0 files if [[ "$flagWriteIntoTodoList" == "true" ]];then check_todo_list_consistency @@ -215,11 +238,16 @@ get_polly_filename() { local device=$1 local date=$2 # /data/level0/polly/pollyxt_lacros/data_zip/201907 - local YYYY=${date:0:4} + local YYYY=${date:0:4} local MM=${date:4:2} local DD=${date:6:2} - local input_path="/data/level0/polly/${device}/data_zip/${YYYY}${MM}" - local searchpattern="${YYYY}_${MM}_${DD}*_*[0-9].nc.zip" + if [[ "$DEVICE" == "martha" ]]; then + local input_path="${BASE_FOLDER}" + local searchpattern="${YYYY}_${MM}_${DD}*_*[0-9].nc" + else + local input_path="${BASE_FOLDER}/${device}/data_zip/${YYYY}${MM}" + local searchpattern="${YYYY}_${MM}_${DD}*_*[0-9].nc.zip" + fi local polly_files=`ls ${input_path}/${searchpattern}` # echo $polly_files # return ${polly_files} @@ -238,7 +266,11 @@ get_polly_filename() { process_history() { DEVICE=$1 DATE=$2 - local POLLY_FOLDER="/data/level0/polly/${DEVICE}" + if [[ "$DEVICE" == "martha" ]]; then + local POLLY_FOLDER="$BASE_FOLDER" + else + local POLLY_FOLDER="${BASE_FOLDER}/${DEVICE}" + fi local POLLY_TYPE=$DEVICE echo -e "\nSettings:\nPOLLY_FOLDER=$POLLY_FOLDER\nPOLLY_TYPE=$POLLY_TYPE\nPICASSO_CONFIG_FILE=$PICASSO_CONFIG_FILE\nSTART_DATE=$DATE\nEND_DATE=$DATE\n\n" @@ -268,7 +300,7 @@ merging() { mkdir -p $OUTPUT_FOLDER ## create folder if not existing, else skip echo "start merging... " - "$PY_FOLDER"python "$PICASSO_DIR_interface"/concat_pollyxt_lvl0.py -t $DATE -d $DEVICE -o $OUTPUT_FOLDER -f ${FORCE_MERGING^} + "$PY_FOLDER"python "$PICASSO_DIR_interface"/concat_pollyxt_lvl0.py -t $DATE -d $DEVICE -o $OUTPUT_FOLDER -r $BASE_FOLDER -f ${FORCE_MERGING^} exit_status=$? @@ -282,12 +314,25 @@ merging() { fi } +process_without_merging() { + DEVICE=$1 + DATE=$2 + process_history $DEVICE $DATE + delete_laserlogbookfile $DEVICE $DATE ## delete laserlogbook-file + delete_entry_from_todo_list $DEVICE $DATE ## delete entry from todo_list file + + +} + write_job_into_todo_list() { ## writing job to todo_list DEVICE=$1 DATE=$2 local OUTPUT_FOLDER=$TODO_FOLDER/$DEVICE/data_zip/${DATE:0:6} + #echo "$OUTPUT_FOLDER" local filename=$(get_polly_filename $DEVICE $DATE) + + #echo "$filename" # echo $filename already_in_list=0 if grep -q "$filename" $PICASSO_TODO_FILE diff --git a/lib/io/loadMeteor.m b/lib/io/loadMeteor.m index 747e1c6d..776a0565 100644 --- a/lib/io/loadMeteor.m +++ b/lib/io/loadMeteor.m @@ -54,6 +54,8 @@ % % HISTORY: % - 2021-05-22: first edition by Zhenping +% - 2025-02-19: modified by C. Jimenez (smoothing added into meteo +% profiles) % % .. Authors: - zhenping@tropos.de @@ -72,6 +74,8 @@ addParameter(p, 'flagReadLess', false, @islogical); addParameter(p, 'method', 'nearest', @ischar); addParameter(p, 'isUseLatestGDAS', true, @islogical); +addParameter(p, 'SmoothMeteo', 7, @isnumeric); + parse(p, mTime, asl, varargin{:}); meteorAttri.dataSource = cell(0); @@ -114,11 +118,11 @@ meteorAttri.datetime = cat(2, meteorAttri.datetime, attri.datetime); % interp the parameters - tempI = interpMeteor(altRaw, tempRaw, asl); - presI = interpMeteor(altRaw, presRaw, asl); - relhI = interpMeteor(altRaw, relhRaw, asl); - winsI = interpMeteor(altRaw, winsRaw, asl); - windI = interpMeteor(altRaw, windRaw, asl); + tempI = transpose(smooth(interpMeteor(altRaw, tempRaw, asl),p.Results.SmoothMeteo,'moving')); + presI = transpose(smooth(interpMeteor(altRaw, presRaw, asl),p.Results.SmoothMeteo,'moving')); + relhI = transpose(smooth(interpMeteor(altRaw, relhRaw, asl),p.Results.SmoothMeteo,'moving')); + winsI = transpose(smooth(interpMeteor(altRaw, winsRaw, asl),p.Results.SmoothMeteo,'moving')); + windI = transpose(smooth(interpMeteor(altRaw, windRaw, asl),p.Results.SmoothMeteo,'moving')); % concatenate the parameters tempQry = cat(1, tempQry, tempI); diff --git a/lib/io/pollySaveAttnBeta.m b/lib/io/pollySaveAttnBeta.m index 7e4b1c9c..9194b00b 100644 --- a/lib/io/pollySaveAttnBeta.m +++ b/lib/io/pollySaveAttnBeta.m @@ -11,8 +11,9 @@ function pollySaveAttnBeta(data) % - 2019-01-10: First Edition by Zhenping % - 2019-05-16: Extended the attributes for all the variables and comply with the ACTRIS convention. % - 2019-09-27: Turn on the netCDF4 compression. +% - 2025-03-14: co and cross polarized components were added to the saving list by Cristofer Jimenez % -% .. Authors: - zhenping@tropos.de +% .. Authors: - zhenping@tropos.de, jimenez@tropos.de missing_value = -999; @@ -38,8 +39,14 @@ function pollySaveAttnBeta(data) varID_height = netcdf.defVar(ncID, 'height', 'NC_FLOAT', dimID_height); varID_tilt_angle = netcdf.defVar(ncID, 'tilt_angle', 'NC_FLOAT', dimID_constant); varID_att_bsc_355 = netcdf.defVar(ncID, 'attenuated_backscatter_355nm', 'NC_FLOAT', [dimID_height, dimID_time]); +varID_att_bsc_para_355 = netcdf.defVar(ncID, 'attenuated_backscatter_para_355nm', 'NC_FLOAT', [dimID_height, dimID_time]); +varID_att_bsc_perp_355 = netcdf.defVar(ncID, 'attenuated_backscatter_perp_355nm', 'NC_FLOAT', [dimID_height, dimID_time]); varID_att_bsc_532 = netcdf.defVar(ncID, 'attenuated_backscatter_532nm', 'NC_FLOAT', [dimID_height, dimID_time]); +varID_att_bsc_para_532 = netcdf.defVar(ncID, 'attenuated_backscatter_para_532nm', 'NC_FLOAT', [dimID_height, dimID_time]); +varID_att_bsc_perp_532 = netcdf.defVar(ncID, 'attenuated_backscatter_perp_532nm', 'NC_FLOAT', [dimID_height, dimID_time]); varID_att_bsc_1064 = netcdf.defVar(ncID, 'attenuated_backscatter_1064nm', 'NC_FLOAT', [dimID_height, dimID_time]); +varID_att_bsc_para_1064 = netcdf.defVar(ncID, 'attenuated_backscatter_para_1064nm', 'NC_FLOAT', [dimID_height, dimID_time]); +varID_att_bsc_perp_1064 = netcdf.defVar(ncID, 'attenuated_backscatter_perp_1064nm', 'NC_FLOAT', [dimID_height, dimID_time]); varID_quality_mask_355 = netcdf.defVar(ncID, 'quality_mask_355nm', 'NC_BYTE', [dimID_height, dimID_time]); varID_quality_mask_532 = netcdf.defVar(ncID, 'quality_mask_532nm', 'NC_BYTE', [dimID_height, dimID_time]); varID_quality_mask_1064 = netcdf.defVar(ncID, 'quality_mask_1064nm', 'NC_BYTE', [dimID_height, dimID_time]); @@ -51,6 +58,12 @@ function pollySaveAttnBeta(data) netcdf.defVarFill(ncID, varID_att_bsc_355, false, missing_value); netcdf.defVarFill(ncID, varID_att_bsc_532, false, missing_value); netcdf.defVarFill(ncID, varID_att_bsc_1064, false, missing_value); +netcdf.defVarFill(ncID, varID_att_bsc_para_355, false, missing_value); +netcdf.defVarFill(ncID, varID_att_bsc_para_532, false, missing_value); +netcdf.defVarFill(ncID, varID_att_bsc_para_1064, false, missing_value); +netcdf.defVarFill(ncID, varID_att_bsc_perp_355, false, missing_value); +netcdf.defVarFill(ncID, varID_att_bsc_perp_532, false, missing_value); +netcdf.defVarFill(ncID, varID_att_bsc_perp_1064, false, missing_value); netcdf.defVarFill(ncID, varID_quality_mask_355, false, 1); netcdf.defVarFill(ncID, varID_quality_mask_532, false, 1); netcdf.defVarFill(ncID, varID_quality_mask_1064, false, 1); @@ -62,6 +75,12 @@ function pollySaveAttnBeta(data) netcdf.defVarDeflate(ncID, varID_att_bsc_355, true, true, 5); netcdf.defVarDeflate(ncID, varID_att_bsc_532, true, true, 5); netcdf.defVarDeflate(ncID, varID_att_bsc_1064, true, true, 5); +netcdf.defVarDeflate(ncID, varID_att_bsc_para_355, true, true, 5); +netcdf.defVarDeflate(ncID, varID_att_bsc_para_532, true, true, 5); +netcdf.defVarDeflate(ncID, varID_att_bsc_para_1064, true, true, 5); +netcdf.defVarDeflate(ncID, varID_att_bsc_perp_355, true, true, 5); +netcdf.defVarDeflate(ncID, varID_att_bsc_perp_532, true, true, 5); +netcdf.defVarDeflate(ncID, varID_att_bsc_perp_1064, true, true, 5); netcdf.defVarDeflate(ncID, varID_quality_mask_355, true, true, 5); netcdf.defVarDeflate(ncID, varID_quality_mask_532, true, true, 5); netcdf.defVarDeflate(ncID, varID_quality_mask_1064, true, true, 5); @@ -87,6 +106,12 @@ function pollySaveAttnBeta(data) netcdf.putVar(ncID, varID_att_bsc_355, single(fillmissing(data.att_beta_355, missing_value))); netcdf.putVar(ncID, varID_att_bsc_532, single(fillmissing(data.att_beta_532, missing_value))); netcdf.putVar(ncID, varID_att_bsc_1064, single(fillmissing(data.att_beta_1064, missing_value))); +netcdf.putVar(ncID, varID_att_bsc_para_355, single(fillmissing(data.att_beta_para_355, missing_value))); +netcdf.putVar(ncID, varID_att_bsc_para_532, single(fillmissing(data.att_beta_para_532, missing_value))); +netcdf.putVar(ncID, varID_att_bsc_para_1064, single(fillmissing(data.att_beta_para_1064, missing_value))); +netcdf.putVar(ncID, varID_att_bsc_perp_355, single(fillmissing(data.att_beta_perp_355, missing_value))); +netcdf.putVar(ncID, varID_att_bsc_perp_532, single(fillmissing(data.att_beta_perp_532, missing_value))); +netcdf.putVar(ncID, varID_att_bsc_perp_1064, single(fillmissing(data.att_beta_perp_1064, missing_value))); netcdf.putVar(ncID, varID_quality_mask_355, int8(fillmissing(data.quality_mask_355, 1))); netcdf.putVar(ncID, varID_quality_mask_532, int8(fillmissing(data.quality_mask_532, 1))); netcdf.putVar(ncID, varID_quality_mask_1064, int8(fillmissing(data.quality_mask_1064, 1))); @@ -185,6 +210,85 @@ function pollySaveAttnBeta(data) % netcdf.putAtt(ncID, varID_att_bsc_1064, 'bias_variable', 'att_beta_1064_bias'); netcdf.putAtt(ncID, varID_att_bsc_1064, 'comment', 'This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter.'); +% att_bsc_para_355 +netcdf.putAtt(ncID, varID_att_bsc_para_355, 'unit', 'sr^-1 m^-1'); +netcdf.putAtt(ncID, varID_att_bsc_para_355, 'unit_html', 'sr-1 m-1'); +netcdf.putAtt(ncID, varID_att_bsc_para_355, 'long_name', 'co polarized attenuated backscatter at 355 nm'); +netcdf.putAtt(ncID, varID_att_bsc_para_355, 'standard_name', 'att_beta_para_355'); +%netcdf.putAtt(ncID, varID_att_bsc_para_355, 'plot_range', PollyConfig.zLim_att_beta_355/1e6); +%netcdf.putAtt(ncID, varID_att_bsc_para_355, 'plot_scale', 'linear'); +netcdf.putAtt(ncID, varID_att_bsc_para_355, 'source', CampaignConfig.name); +%netcdf.putAtt(ncID, varID_att_bsc_para_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355); +% netcdf.putAtt(ncID, varID_att_bsc_para_355, 'error_variable', 'att_beta_355_error'); +% netcdf.putAtt(ncID, varID_att_bsc_para_355, 'bias_variable', 'att_beta_355_bias'); +netcdf.putAtt(ncID, varID_att_bsc_para_355, 'comment', 'This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter. It also make uses of the volume depolarization ratio to get the co polarized component'); + +% att_bsc_para_532 +netcdf.putAtt(ncID, varID_att_bsc_para_532, 'unit', 'sr^-1 m^-1'); +netcdf.putAtt(ncID, varID_att_bsc_para_532, 'unit_html', 'sr-1 m-1'); +netcdf.putAtt(ncID, varID_att_bsc_para_532, 'long_name', 'co polarized attenuated backscatter at 532 nm'); +netcdf.putAtt(ncID, varID_att_bsc_para_532, 'standard_name', 'att_beta_532'); +%netcdf.putAtt(ncID, varID_att_bsc_para_532, 'plot_range', PollyConfig.zLim_att_beta_532/1e6); +%netcdf.putAtt(ncID, varID_att_bsc_para_532, 'plot_scale', 'linear'); +netcdf.putAtt(ncID, varID_att_bsc_para_532, 'source', CampaignConfig.name); +%netcdf.putAtt(ncID, varID_att_bsc_para_532, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed532); +% netcdf.putAtt(ncID, varID_att_bsc_para_532, 'error_variable', 'att_beta_532_error'); +% netcdf.putAtt(ncID, varID_att_bsc_para_532, 'bias_variable', 'att_beta_532_bias'); +netcdf.putAtt(ncID, varID_att_bsc_para_532, 'comment', 'This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter. It also make uses of the volume depolarization ratio to get the co polarized component'); + +% att_bsc_para_1064 +netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'unit', 'sr^-1 m^-1'); +netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'unit_html', 'sr-1 m-1'); +netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'long_name', 'co polarized attenuated backscatter at 1064 nm'); +netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'standard_name', 'att_beta_1064'); +%netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'plot_range', PollyConfig.zLim_att_beta_1064/1e6); +%netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'plot_scale', 'linear'); +netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'source', CampaignConfig.name); +%netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed1064); +% netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'error_variable', 'att_beta_1064_error'); +% netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'bias_variable', 'att_beta_1064_bias'); +netcdf.putAtt(ncID, varID_att_bsc_para_1064, 'comment', 'This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter. It also make uses of the volume depolarization ratio to get the cross polarized component'); + + +% att_bsc_perp_355 +netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'unit', 'sr^-1 m^-1'); +netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'unit_html', 'sr-1 m-1'); +netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'long_name', 'cross polarized attenuated backscatter at 355 nm'); +netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'standard_name', 'att_beta_perp_355'); +%netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'plot_range', PollyConfig.zLim_att_beta_355/1e6); +%netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'plot_scale', 'linear'); +netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'source', CampaignConfig.name); +%netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355); +% netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'error_variable', 'att_beta_355_error'); +% netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'bias_variable', 'att_beta_355_bias'); +netcdf.putAtt(ncID, varID_att_bsc_perp_355, 'comment', 'This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter. It also make uses of the volume depolarization ratio to get the cross polarized component'); + +% att_bsc_perp_532 +netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'unit', 'sr^-1 m^-1'); +netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'unit_html', 'sr-1 m-1'); +netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'long_name', 'cross polarized attenuated backscatter at 532 nm'); +netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'standard_name', 'att_beta_532'); +%netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'plot_range', PollyConfig.zLim_att_beta_532/1e6); +%netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'plot_scale', 'linear'); +netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'source', CampaignConfig.name); +%netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed532); +% netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'error_variable', 'att_beta_532_error'); +% netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'bias_variable', 'att_beta_532_bias'); +netcdf.putAtt(ncID, varID_att_bsc_perp_532, 'comment', 'This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter. It also make uses of the volume depolarization ratio to get the cross polarized component'); + +% att_bsc_perp_1064 +netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'unit', 'sr^-1 m^-1'); +netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'unit_html', 'sr-1 m-1'); +netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'long_name', 'cross polarized attenuated backscatter at 1064 nm'); +netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'standard_name', 'att_beta_1064'); +%netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'plot_range', PollyConfig.zLim_att_beta_1064/1e6); +%netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'plot_scale', 'linear'); +netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'source', CampaignConfig.name); +%netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed1064); +% netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'error_variable', 'att_beta_1064_error'); +% netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'bias_variable', 'att_beta_1064_bias'); +netcdf.putAtt(ncID, varID_att_bsc_perp_1064, 'comment', 'This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter.'); + % quality_mask_355 netcdf.putAtt(ncID, varID_quality_mask_355, 'unit', ''); netcdf.putAtt(ncID, varID_quality_mask_355, 'long_name', 'quality mask for attenuated backscatter at 355 nm'); diff --git a/lib/io/pollySavePOLIPHON.m b/lib/io/pollySavePOLIPHON.m index 97988560..66001633 100644 --- a/lib/io/pollySavePOLIPHON.m +++ b/lib/io/pollySavePOLIPHON.m @@ -7,9 +7,12 @@ function pollySavePOLIPHON(data, POLIPHON1) % % HISTORY: % - 2023-06-26: first edition by Athena A. Floutsi -% - 2024-05-17: updates on startTime/endTime writing and on variable name +% - 2024-05-17: updates on startTime/endTime writing and on variable name +% - 2025-09-15: updates to fix bug that didn't allow proper time +% writting % % .. Authors: - floutsi@tropos.de + global PicassoConfig CampaignConfig PollyDataInfo PollyConfig missing_value = -999; @@ -28,7 +31,6 @@ function pollySavePOLIPHON(data, POLIPHON1) dimID_method = netcdf.defDim(ncID, 'method', 1); dimID_time = netcdf.defDim(ncID, 'time', length(data.mTime)); - % define variables varID_altitude = netcdf.defVar(ncID, 'altitude', 'NC_FLOAT', dimID_method); varID_longitude = netcdf.defVar(ncID, 'longitude', 'NC_FLOAT', dimID_method); @@ -36,7 +38,7 @@ function pollySavePOLIPHON(data, POLIPHON1) varID_startTime = netcdf.defVar(ncID, 'start_time', 'NC_DOUBLE', dimID_method); varID_endTime = netcdf.defVar(ncID, 'end_time', 'NC_DOUBLE', dimID_method); varID_height = netcdf.defVar(ncID, 'height', 'NC_FLOAT', dimID_height); -varID_time = netcdf.defVar(ncID, 'time', 'NC_FLOAT', dimID_time); +varID_time = netcdf.defVar(ncID, 'time', 'NC_DOUBLE', dimID_time); varID_aerBsc_klett_355 = netcdf.defVar(ncID, 'aerBsc_klett_355', 'NC_FLOAT', dimID_height); varID_aerBscStd_klett_355 = netcdf.defVar(ncID, 'uncertainty_aerBsc_klett_355', 'NC_FLOAT', dimID_height); @@ -163,6 +165,7 @@ function pollySavePOLIPHON(data, POLIPHON1) netcdf.putVar(ncID, varID_startTime, datenum_2_unix_timestamp(startTime)); netcdf.putVar(ncID, varID_endTime, datenum_2_unix_timestamp(endTime)); netcdf.putVar(ncID, varID_height, single(data.height)); +netcdf.putVar(ncID, varID_time, datenum_2_unix_timestamp(data.mTime)); netcdf.putVar(ncID, varID_aerBsc_klett_355, single(fillmissing(data.aerBsc355_klett(iGrp, :), missing_value))); netcdf.putVar(ncID, varID_aerBscStd_klett_355, single(fillmissing(data.aerBscStd355_klett(iGrp, :), missing_value))); diff --git a/lib/io/pollySaveProfiles_QC.m b/lib/io/pollySaveProfiles_QC.m index 1f0c7aad..0bb5d086 100644 --- a/lib/io/pollySaveProfiles_QC.m +++ b/lib/io/pollySaveProfiles_QC.m @@ -88,8 +88,10 @@ function pollySaveProfiles_QC(data) data.aerBsc532_klett(iGrp,(data.height <= PollyConfig.heightFullOverlap(6))) = data.aerBsc532_NR_klett(iGrp,(data.height <= PollyConfig.heightFullOverlap(6))); %IR FF data.aerBsc1064_raman(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; - data.aerExt1064_raman(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; - data.LR1064_raman(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; + %data.aerExt1064_raman(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; + %data.LR1064_raman(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; + data.aerExt1064_raman(iGrp,:) = missing_value; + data.LR1064_raman(iGrp,:) = missing_value; data.pdr1064_raman(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; data.pdr1064_klett(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; data.aerBsc1064_klett(iGrp,(data.height <= PollyConfig.heightFullOverlap(8))) = missing_value; diff --git a/lib/io/readMETncERA5.m b/lib/io/readMETncERA5.m new file mode 100644 index 00000000..d0cbc1d3 --- /dev/null +++ b/lib/io/readMETncERA5.m @@ -0,0 +1,136 @@ +function [alt, temp, pres, relh, wins, wind, fname] = readMETncERA5(tRange, site, folder, varargin) +% readMETnccloudnet read the cloudnet ecmwf netcdf file +% +% EXAMPLE: +% [alt, temp, pres, relh] = readMETnccloudnet(tRange, site, folder) +% +% INPUTS: +% tRange: 2-element array +% search range. +% site: char +% the location for gdas1. Our server will automatically produce the +% gdas1 products for all our pollynet location. You can find it in +% /lacroshome/cloudnet/data/model/gdas1 +% +% KEYWORDS: +% isUseLatestGDAS: logical +% whether to search the latest available GDAS profile (default: false). +% +% OUTPUTS: +% alt: array +% altitute (above ground) for each range bin. [m] +% temp: array +% temperature for each range bin. If no valid data, NaN will be +% filled. [C] +% pres: array +% pressure for each range bin. If no valid data, NaN will be filled. +% [hPa] +% rh: array +% relative humidity for each range bin. If no valid data, NaN will be +% filled. [%] +% wins: array +% wind speed [m/s] +% wind: array +% wind direction. [degree] +% fname: char +% filename (for legacy reasons the name is not changed). +% +% HISTORY: +% - 2023-05-13.:First implementation by martin-rdz +% +% .. Authors: - radenz@tropos.de + +p = inputParser; +p.KeepUnmatched = true; + +addRequired(p, 'tRange', @isnumeric); +addRequired(p, 'site', @ischar); +addRequired(p, 'folder', @ischar); + +parse(p, tRange, site, folder, varargin{:}); + +midTime = mean(tRange); + +[thisyear, thismonth, thisday, thishour, ~, ~] = datevec(midTime); +% /oceanethome/model/ecmwf/profiles/ecmwf/2023/20230512_neumayer_ecmwf.nc +fname = fullfile(folder, sprintf('%04d', thisyear), ... + sprintf('%04d%02d%02d_%s_era5-1-12.nc', ... + thisyear, thismonth, thisday, site)); + +disp(fname); + +% preallocate +alt = NaN; +pres = NaN; +temp = NaN; +relh = NaN; +uwd =NaN; +vwd = NaN; +wins = NaN; +wind =NaN; + +filenameList = dir(fname); +if numel(filenameList) == 0 + fprintf('ECMWF File (%s) does not exist.\n', fname); + return; +end +% Open the netCDF file +ncid = netcdf.open(fname, 'NOWRITE'); + + +% Get information about the netCDF file +nc_info = ncinfo(fname); +% Extract the variable names +var_names = {nc_info.Variables.Name}; +% Display the variable names +%disp(var_names); + +% Get the variable ID for the time variable and the variable of interest +%time_varid = netcdf.inqVarID(ncid, 'time'); +%data_varid = netcdf.inqVarID(ncid, 'variable_name'); + +% Read in the time variable and the variable of interest +time = ncread(fname, 'time'); +[~, ~, ~, hour, minute, second] = datevec(midTime); +hour_fraction = hour + minute/60 + second/3600; + +disp(hour_fraction); + +% Find the index of the time value closest to the specified time +[~, index] = min(abs(time - hour_fraction)); + +height = ncread(fname, 'height'); +alt = height(:, index); + +% Extract the corresponding data +data = ncread(fname, 'pressure'); +pres = data(:, index) ./ 100.; + +data = ncread(fname, 'temperature'); +temp = data(:, index) - 273.15; + +data = ncread(fname, 'rh'); +relh = data(:, index) .* 100; + +data = ncread(fname, 'uwind'); +uwd = data(:, index); + +data = ncread(fname, 'vwind'); +vwd = data(:, index); + +wins = sqrt(uwd.^2 + vwd.^2); +wind = atan2(vwd,uwd)*(180/pi); + +% Close the netCDF file +netcdf.close(ncid); + +%#[pres, alt, temp, relh, wind, wins] = ceilo_bsc_ModelSonde(gdas1file); + +%pres = NaN; +%alt = NaN; +%temp = NaN; +%relh = NaN; +%wind = NaN; +%wins = NaN; + +end diff --git a/lib/io/readMeteor.m b/lib/io/readMeteor.m index fb6a8110..f7427d3e 100644 --- a/lib/io/readMeteor.m +++ b/lib/io/readMeteor.m @@ -109,7 +109,7 @@ wind = NaN(size(temp)); wins = NaN(size(temp)); alt = alt * 1e3; % convert to [m] - temp = temp - 273.17; % convert to [\circC] + temp = temp - 273.15; % convert to [\circC] attri.dataSource = p.Results.meteorDataSource; attri.URL = ''; attri.datetime = datenum(0,1,0,0,0,0); @@ -189,6 +189,24 @@ attri.datetime = measTime; end +case 'era5' + [alt, temp, pres, relh, wins, wind, gdas1File] = readMETncERA5(measTime, ... + p.Results.gdas1Site, p.Results.meteo_folder, ... + 'isUseLatestGDAS', p.Results.isUseLatestGDAS); + + if isempty(alt) + alt = []; + temp = []; + pres = []; + relh = []; + wins = []; + wind = []; + else + attri.dataSource = p.Results.meteorDataSource; + attri.URL = gdas1File; + attri.datetime = measTime; + end + otherwise error('Unknown meteorological data source.\n%s\n', ... p.Results.meteorDataSource) @@ -217,7 +235,7 @@ [alt, ~, ~, temp, pres] = atmo(60, 0.03, 1); alt = alt * 1e3; pres = pres / 1e2; - temp = temp - 273.17; % convert to [\circC] + temp = temp - 273.15; % convert to [\circC] relh = NaN(size(temp)); wins = NaN(size(temp)); wind = NaN(size(temp)); diff --git a/lib/io/readSonde.m b/lib/io/readSonde.m index 79b78cf2..7656761e 100644 --- a/lib/io/readSonde.m +++ b/lib/io/readSonde.m @@ -135,7 +135,17 @@ relh = str2num(char(strData{1,5}(22:end))); wind = str2num(char(strData{1,6}(22:end))); wins = str2num(char(strData{1,7}(22:end))); - +case 4 % Meteor radiosonde standard file + + thisFilename = basename(file); + datetime = datenum(thisFilename(end-15:end-3), 'yyyymmdd_HHMMSS'); + + alt = ncread(file, 'altitude'); + temp = ncread(file, 'temperature'); + pres = ncread(file, 'pressure'); + relh = ncread(file, 'RH'); + wind = ncread(file, 'wind_direction'); + wins = ncread(file, 'wind_speed'); otherwise error('Unknown fileType %d', fileType); diff --git a/lib/misc/sondeSearch.m b/lib/misc/sondeSearch.m index 8a7a402d..9e037950 100644 --- a/lib/misc/sondeSearch.m +++ b/lib/misc/sondeSearch.m @@ -131,7 +131,35 @@ datestr(sondeTime(indxSondeFile), 'yyyymmdd HH:MM:SS')); end +case 4 % meteor radiosonde standard file + %% list all files + sondeFileList = listfile(sondeFolder, '_\d\d\d\d\d\d\d\d_\d\d\d\d\d\d'); + if isempty(sondeFileList) + warning(['No required radiosonde files was found in the sonde folder. ' ... + 'Please go to the folder below to have a look.\n%s'], sondeFolder); + return; + end + + %% parse the radiosonde time + sondeTime = NaN(size(sondeFileList)); + for iFile = 1:length(sondeFileList) + filenameISondeFile = basename(sondeFileList{iFile}); + sondeTime(iFile) = datenum(filenameISondeFile(end-17:end-3), 'yyyymmdd_HHMMSS');%(end-1:end-3) + end + + %% search the sonding file which is closest to the measurement time + deltaTime = abs(sondeTime - measurementTime); + [minDeltaTime, indxSondeFile] = min(deltaTime); + % determine whether the time lapse is out of range (max T diff: 1 day) + if minDeltaTime < datenum(0, 1, 1, 0, 0, 0) + sondeFile = sondeFileList{indxSondeFile}; + else + warning(['There was no sonde launching within 1 day.\n' ... + 'The measurement time: %s\nThe closest time of sonding: %s'], ... + datestr(measurementTime, 'yyyymmdd HH:MM:SS'), ... + datestr(sondeTime(indxSondeFile), 'yyyymmdd HH:MM:SS')); + end otherwise error('Unknown fileType %d', fileType); end diff --git a/lib/preprocess/pollyDTCor.m b/lib/preprocess/pollyDTCor.m index 54157719..9b2b56c3 100644 --- a/lib/preprocess/pollyDTCor.m +++ b/lib/preprocess/pollyDTCor.m @@ -34,7 +34,7 @@ % % HISTORY: % - 2021-05-16: first edition by Zhenping -% +% - 2025-02-19: edition by Cristofer Jimenez % .. Authors: - zhenping@tropos.de p = inputParser; @@ -56,6 +56,9 @@ reshape(mShots, size(mShots, 1), 1, size(mShots, 2)), ... [1, size(sigI, 2), 1]); % reshape mShots to the same dimensions of sigI +MShots_back = repmat(... + nanmedian(mShots,2),[1, size(sigI, 2), size(mShots, 2)]); %take median value over time to recalculate Photocount signal normalized to a common integration time. + %% Deadtime correction if p.Results.flagDeadTimeCorrection PCR = sigI ./ MShots * 150.0 ./ hRes; % convert photon counts to phton @@ -66,7 +69,7 @@ for iCh = 1:size(sigI, 1) PCR_Cor = polyval(p.Results.deadtime(iCh, end:-1:1), ... PCR(iCh, :, :)); - sigO(iCh, :, :) = PCR_Cor / (150.0 / hRes) .* MShots(iCh, :, :); + sigO(iCh, :, :) = PCR_Cor / (150.0 / hRes) .* MShots_back(iCh, :, :); end % nonparalyzable correction @@ -75,7 +78,7 @@ PCR_Cor = PCR(iCh, :, :) ./ ... (1.0 - p.Results.deadtimeParams(iCh) * 1e-3 * ... PCR(iCh, :, :)); - sigO(iCh, :, :) = PCR_Cor / (150.0 / hRes) .* MShots(iCh, :, :); + sigO(iCh, :, :) = PCR_Cor / (150.0 / hRes) .* MShots_back(iCh, :, :); end % user defined deadtime. @@ -85,7 +88,7 @@ for iCh = 1:size(sigI, 1) PCR_Cor = polyval(p.Results.deadtimeParams(iCh, end:-1:1), ... PCR(iCh, :, :)); - sigO(iCh, :, :) = PCR_Cor / (150.0 / hRes) .* MShots(iCh, :, :); + sigO(iCh, :, :) = PCR_Cor / (150.0 / hRes) .* MShots_back(iCh, :, :); end else warning(['User defined deadtime parameters were not found. ', ... diff --git a/lib/preprocess/pollyPreprocess.m b/lib/preprocess/pollyPreprocess.m index d5a5e0ff..fc2424a3 100644 --- a/lib/preprocess/pollyPreprocess.m +++ b/lib/preprocess/pollyPreprocess.m @@ -298,13 +298,14 @@ %% Height (first bin height correction) data.height = double((0:(size(data.signal, 2)-1)) * data.hRes * ... - cos(data.zenithAng / 180 * pi) + config.firstBinHeight); % [m] + cos(data.zenithAng / 180 * pi) + config.firstBinHeight * cos(data.zenithAng / 180 * pi)); % [m] data.alt = double(data.height + config.asl); % geopotential height % distance between range bin and system. -data.distance0 = double(data.height ./ cos(data.zenithAng / 180 * pi)); +%data.distance0 = double(data.height ./ cos(data.zenithAng / 180 * pi)); +data.distance0 = double((0:(size(data.signal, 2)-1)) * data.hRes + config.firstBinHeight); %% Temperature effect correction (for Raman signal) -if config.flagSigTempCor +if config.flagSigTempCor %attention this function might completely change your absolute values -- not recommended temperature = loadMeteor(mean(data.mTime), data.alt, ... 'meteorDataSource', config.meteorDataSource, ... 'gdas1Site', config.gdas1Site, ... @@ -314,7 +315,7 @@ 'radiosondeType', config.radiosondeType, ... 'method', 'linear', ... 'isUseLatestGDAS', config.flagUseLatestGDAS); - absTemp = temperature + 273.17; + absTemp = temperature + 273.15; for iCh = 1:size(data.signal, 1) leadingChar = config.tempCorFunc{iCh}(1); diff --git a/lib/qc/pollyOVLCalcRaman.m b/lib/qc/pollyOVLCalcRaman.m index 9b4cd346..632c89b4 100644 --- a/lib/qc/pollyOVLCalcRaman.m +++ b/lib/qc/pollyOVLCalcRaman.m @@ -90,8 +90,8 @@ if size(p.Results.aerBsc,1)>0 - [mBscRa, mExtRa] = rayleigh_scattering(Lambda_Ra, p.Results.pressure, p.Results.temperature + 273.17, 380, 70); - [~, mExtel] = rayleigh_scattering(Lambda_el, p.Results.pressure, p.Results.temperature + 273.17, 380, 70); % mBscel + [mBscRa, mExtRa] = rayleigh_scattering(Lambda_Ra, p.Results.pressure, p.Results.temperature + 273.15, 380, 70); + [~, mExtel] = rayleigh_scattering(Lambda_el, p.Results.pressure, p.Results.temperature + 273.15, 380, 70); % mBscel mExtRa=mean(mExtRa,1); mBscRa=mean(mBscRa,1); diff --git a/lib/retrievals/pollyLR.m b/lib/retrievals/pollyLR.m index 8585f072..2379945d 100644 --- a/lib/retrievals/pollyLR.m +++ b/lib/retrievals/pollyLR.m @@ -1,4 +1,4 @@ -function [aerLR, effRes, aerLRStd] = pollyLR(aerExt, aerBsc, varargin) +function [aerLR, aerLRStd, effRes] = pollyLR(aerExt, aerBsc, varargin) % POLLYLR calculate aerosol lidar ratio. % % USAGE: @@ -82,6 +82,6 @@ aerBscStd = p.Results.aerBscStd; end -aerLRStd = aerLR .* sqrt(aerExtStd.^2 ./ aerExt.^2 + aerBscStd.^2 ./ aerBsc.^2); +aerLRStd = real(aerLR .* sqrt(aerExtStd.^2 ./ aerExt.^2 + aerBscStd.^2 ./ aerBsc.^2)); end \ No newline at end of file diff --git a/lib/visualization/pollyDisplayOCProfiles.m b/lib/visualization/pollyDisplayOCProfiles.m index eb4b31ce..2e5a0144 100644 --- a/lib/visualization/pollyDisplayOCProfiles.m +++ b/lib/visualization/pollyDisplayOCProfiles.m @@ -131,9 +131,9 @@ function pollyDisplayOCProfiles(data) rcs1064 = transpose(smooth(rcs1064, smoothWin_1064)); % molecule signal - [molBsc355, molExt355] = rayleigh_scattering(355, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [molBsc532, molExt532] = rayleigh_scattering(532, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [molBsc1064, molExt1064] = rayleigh_scattering(1064, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); + [molBsc355, molExt355] = rayleigh_scattering(355, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [molBsc532, molExt532] = rayleigh_scattering(532, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [molBsc1064, molExt1064] = rayleigh_scattering(1064, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); molRCS355 = molBsc355 .* exp(- 2 * cumsum(molExt355 .* [data.distance0(1), diff(data.distance0)])); molRCS532 = molBsc532 .* exp(- 2 * cumsum(molExt532 .* [data.distance0(1), diff(data.distance0)])); molRCS1064 = molBsc1064 .* exp(- 2 * cumsum(molExt1064 .* [data.distance0(1), diff(data.distance0)])); diff --git a/lib/visualization/pollyDisplayProfiles.m b/lib/visualization/pollyDisplayProfiles.m index b3bbf9cc..e6ae3de8 100644 --- a/lib/visualization/pollyDisplayProfiles.m +++ b/lib/visualization/pollyDisplayProfiles.m @@ -134,9 +134,9 @@ function pollyDisplayProfiles(data) rcs1064 = transpose(smooth(rcs1064, smoothWin_1064)); % molecule signal - [molBsc355, molExt355] = rayleigh_scattering(355, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [molBsc532, molExt532] = rayleigh_scattering(532, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); - [molBsc1064, molExt1064] = rayleigh_scattering(1064, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.17, 380, 70); + [molBsc355, molExt355] = rayleigh_scattering(355, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [molBsc532, molExt532] = rayleigh_scattering(532, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); + [molBsc1064, molExt1064] = rayleigh_scattering(1064, data.pressure(iGrp, :), data.temperature(iGrp, :) + 273.15, 380, 70); molRCS355 = molBsc355 .* exp(- 2 * cumsum(molExt355 .* [data.distance0(1), diff(data.distance0)])); molRCS532 = molBsc532 .* exp(- 2 * cumsum(molExt532 .* [data.distance0(1), diff(data.distance0)])); molRCS1064 = molBsc1064 .* exp(- 2 * cumsum(molExt1064 .* [data.distance0(1), diff(data.distance0)])); diff --git a/lib/visualization/pypolly_display_3d_plots.py b/lib/visualization/pypolly_display_3d_plots.py index d2b18fae..8e5b7368 100644 --- a/lib/visualization/pypolly_display_3d_plots.py +++ b/lib/visualization/pypolly_display_3d_plots.py @@ -1819,8 +1819,10 @@ def pollyDisplayQR(nc_dict,config_dict, polly_conf_dict, saveFolder, q_param, q_ ## read from global config file yLim = polly_conf_dict['yLim_Quasi_Params'] - if q_param == "angexp" or q_param == "bsc_1064": + if q_param == "bsc_1064": zLim = polly_conf_dict['zLim_quasi_beta_1064'] + elif q_param == "angexp": + zLim = polly_conf_dict['zLim_quasi_ANG'] elif q_param == "bsc_532": zLim = polly_conf_dict['zLim_quasi_beta_532'] elif q_param == "par_depol_532": @@ -1890,6 +1892,7 @@ def pollyDisplayQR(nc_dict,config_dict, polly_conf_dict, saveFolder, q_param, q_ quasi_title = f'Quasi particle depolarization ratio ({q_version}) at 532 nm from {pollyVersion} at {location}' saveFilename = os.path.join(saveFolder,plotfile) + print(plotfile) ## fill time gaps in att_bsc matrix matrix, quality_mask = readout.fill_time_gaps_of_matrix(time, matrix, quality_mask) @@ -2017,6 +2020,12 @@ def pollyDisplayQR(nc_dict,config_dict, polly_conf_dict, saveFolder, q_param, q_ plt.close() ## write2donefilelist + if q_version == "V1": + prod_type_q = "" + elif q_version == "V2": + q = q_version.lower() + prod_type_q = f"_{q}" + readout.write2donefilelist_dict(donefilelist_dict = donefilelist_dict, lidar = pollyVersion, location = nc_dict['location'], @@ -2034,7 +2043,7 @@ def pollyDisplayQR(nc_dict,config_dict, polly_conf_dict, saveFolder, q_param, q_ GDAS_timestamp = f"{datetime.utcfromtimestamp(int(nc_dict['time'][0])).strftime('%Y%m%d')} 12:00:00", lidar_ratio = 50, software_version = version, - product_type = f'{prodtype}', + product_type = f'{prodtype}{prod_type_q}', product_starttime = datetime.utcfromtimestamp(int(nc_dict['time'][0])).strftime('%Y%m%d %H:%M:%S'), product_stoptime = datetime.utcfromtimestamp(int(nc_dict['time'][-1])).strftime('%Y%m%d %H:%M:%S') ) @@ -2067,7 +2076,7 @@ def pollyDisplay_Overlap(nc_dict,config_dict,polly_conf_dict,outdir,donefilelist ## read from global config file xLim = [-0.1, 1.1] - yLim = polly_conf_dict['yLim_att_beta_NR'] + yLim = polly_conf_dict['yLim_overlap'] partnerLabel = polly_conf_dict['partnerLabel'] imgFormat = polly_conf_dict['imgFormat'] diff --git a/lib/visualization/pypolly_display_profiles.py b/lib/visualization/pypolly_display_profiles.py index 5e6bf77c..ce842838 100644 --- a/lib/visualization/pypolly_display_profiles.py +++ b/lib/visualization/pypolly_display_profiles.py @@ -865,9 +865,10 @@ def pollyDisplay_profile_summary(nc_dict_profile,nc_dict_profile_NR,config_dict, 2022-09-01. First edition by Andi """ - if not nc_dict_profile : + if not nc_dict_profile: return + ## read from config file figDPI = config_dict['figDPI'] flagWatermarkOn = config_dict['flagWatermarkOn'] @@ -922,13 +923,16 @@ def pollyDisplay_profile_summary(nc_dict_profile,nc_dict_profile_NR,config_dict, param_dict = { "backscatter": - {"FR": ['aerBsc_raman_355','aerBsc_raman_532','aerBsc_raman_1064'], + {"FR": + ['aerBsc_raman_355','aerBsc_raman_532','aerBsc_raman_1064','aerBsc_RR_355','aerBsc_RR_532','aerBsc_RR_1064'], "NR":['aerBsc_raman_355','aerBsc_raman_532'] }, - "extinction": {"FR": ['aerExt_raman_355','aerExt_raman_532','aerExt_raman_1064'], + "extinction": {"FR": + ['aerExt_raman_355','aerExt_raman_532','aerExt_raman_1064','aerExt_RR_355','aerExt_RR_532','aerExt_RR_1064'], "NR": ['aerExt_raman_355','aerExt_raman_532'] }, - "lidarratio": {"FR": ['aerLR_raman_355','aerLR_raman_532','aerLR_raman_1064'], + "lidarratio": {"FR": + ['aerLR_raman_355','aerLR_raman_532','aerLR_raman_1064','aerLR_RR_355','aerLR_RR_532','aerLR_RR_1064'], "NR": ['aerLR_raman_355','aerLR_raman_532'] }, "angstroem": {"FR": ['AE_beta_355_532_Raman','AE_beta_532_1064_Raman','AE_parExt_355_532_Raman'], @@ -970,16 +974,19 @@ def pollyDisplay_profile_summary(nc_dict_profile,nc_dict_profile_NR,config_dict, if method == 'raman': param_dict = { "backscatter": - {"FR": ['aerBsc_raman_355','aerBsc_raman_532','aerBsc_raman_1064'], + {"FR": + ['aerBsc_raman_355','aerBsc_raman_532','aerBsc_raman_1064','aerBsc_RR_355','aerBsc_RR_532','aerBsc_RR_1064'], "NR": [] }, - "extinction": {"FR": ['aerExt_raman_355','aerExt_raman_532','aerExt_raman_1064'], + "extinction": {"FR": + ['aerExt_raman_355','aerExt_raman_532','aerExt_raman_1064','aerExt_RR_355','aerExt_RR_532','aerExt_RR_1064'], "NR": [] }, "lidarratio": {"FR": ['aerLR_raman_355','aerLR_raman_532','aerLR_raman_1064'], "NR": [] }, - "angstroem": {"FR": ['AE_beta_355_532_Raman','AE_beta_532_1064_Raman','AE_parExt_355_532_Raman'], + "angstroem": {"FR": + ['aerLR_raman_355','aerLR_raman_532','aerLR_raman_1064','aerLR_RR_355','aerLR_RR_532','aerLR_RR_1064'], "NR": [] }, "depolarization": {"FR": ['parDepol_raman_355','parDepol_raman_532','parDepol_raman_1064'], @@ -1031,9 +1038,13 @@ def plotting_procedure(col,param_dict,parameter,xlabel,xlim=[0,1],ylim=[0,1],sca if parameter == 'angstroem': color_ls = ['orange','magenta','black'] else: - color_ls = ['blue','green','red'] + color_ls = ['blue','green','red','purple','olive','lightsalmon'] - line_style = '-' + if 'RR' in p: + line_style = 'dashed' + else: + line_style = '-' + #line_style = '-' if parameter == 'wvmr': # ax2 = ax[col].secondary_xaxis('top') @@ -1322,13 +1333,16 @@ def pollyDisplay_profile_summary_QC(nc_dict_profile,config_dict,polly_conf_dict, param_dict = { "backscatter": {"Klett": ['aerBsc_klett_355','aerBsc_klett_532','aerBsc_klett_1064'], - "Raman":['aerBsc_raman_355','aerBsc_raman_532','aerBsc_raman_1064'] + "Raman":['aerBsc_raman_355','aerBsc_raman_532','aerBsc_raman_1064','aerBsc_RR_355','aerBsc_RR_532','aerBsc_RR_1064'] +# "Raman":['aerBsc_raman_355','aerBsc_raman_532','aerBsc_raman_1064'] }, "extinction": {"Klett": ['aerBsc_klett_355','aerBsc_klett_532','aerBsc_klett_1064'], - "Raman": ['aerExt_raman_355','aerExt_raman_532','aerExt_raman_1064'] + #"Raman": ['aerExt_raman_355','aerExt_raman_532','aerExt_raman_1064'] + "Raman": ['aerExt_raman_355','aerExt_raman_532','aerExt_raman_1064','aerExt_RR_355','aerExt_RR_532','aerExt_RR_1064'], }, "lidarratio": {"Klett": [], - "Raman": ['aerLR_raman_355','aerLR_raman_532','aerLR_raman_1064'] + "Raman": ['aerLR_raman_355','aerLR_raman_532','aerLR_raman_1064','aerLR_RR_355','aerLR_RR_532','aerLR_RR_1064'], + #"Raman": ['aerLR_raman_355','aerLR_raman_532','aerLR_raman_1064'] }, "angstroem": {"Klett": [], "Raman": ['AE_beta_355_532_Raman','AE_beta_532_1064_Raman','AE_parExt_355_532_Raman'] @@ -1359,9 +1373,13 @@ def plotting_procedure_QC(col,param_dict,parameter,xlabel,xlim=[0,1],ylim=[0,1], if parameter == 'angstroem': color_ls = ['orange','magenta','black'] else: - color_ls = ['blue','green','red'] + color_ls = ['blue','green','red','purple','olive','lightsalmon'] - line_style = '-' + + if 'RR' in p: + line_style = 'dashed' + else: + line_style = '-' if parameter == 'wvmr': # ax2 = ax[col].secondary_xaxis('top') diff --git a/lib/visualization/pypolly_readout.py b/lib/visualization/pypolly_readout.py index 480505b8..9ec83a18 100644 --- a/lib/visualization/pypolly_readout.py +++ b/lib/visualization/pypolly_readout.py @@ -278,7 +278,8 @@ def read_config(configfile): def read_excel_config_file(excel_file, timestamp, device): pd.set_option('display.width', 1500) pd.set_option('display.max_columns', None) - excel_file_ds = pd.read_excel(f'{excel_file}', engine='openpyxl',usecols = 'A:Z') + #excel_file_ds = pd.read_excel(f'{excel_file}', engine='openpyxl',usecols = 'A:Z') + excel_file_ds = pd.read_excel(f'{excel_file}', engine='openpyxl') print(excel_file) ## search for timerange for given timestamp filtered_device = excel_file_ds.loc[(excel_file_ds['Instrument'] == device)]