From 28f4bb4d54d260db05350da1de864994dabf28ac Mon Sep 17 00:00:00 2001 From: Iwan Firmawan Date: Wed, 23 Oct 2024 10:42:02 +0700 Subject: [PATCH] Feature/1 update 2024 08 19 (#2) * [#1] Update scripts on libs * [#1] Update project settings * [#1] Update pattern settings - Rename MOD11C3 to MOD21C3 * [#1] Update execute_steps.bat * [#1] Update step 0000 * [#1] Update step 0101 * [#1] Update step 0303 * [#1] Update source folder to include mapping/data & mapping/output/maps * [#1] Update scripts all FinalMapOutput * [#1] Fix dir structure for source/mapping --- .gitignore | 8 +- FinalMapOutput-ESRI.py | 12 +- FinalMapOutput-POR-ESRI.py | 41 ++++--- STEP_0000_execute_all_steps.py | 1 + STEP_0101_read_hdf_create_LST_anom_netcdf.py | 35 +++--- STEP_0303_export_ranking_data_rasters.py | 112 ++++++------------ cdi_directory_settings.conf | 2 +- cdi_pattern_settings.conf | 2 +- cdi_project_settings.conf | 16 +-- execute_steps.bat | 5 +- libs/file_operations.py | 17 ++- libs/hdf_functions.py | 2 +- libs/statistics_operations.py | 33 +++--- libs/subgrid_calculations.py | 2 +- .../{MOD11C3_LST => MOD21C3_LST}/.gitkeep | 0 source/mapping/.gitkeep | 0 source/mapping/data/.gitkeep | 0 source/mapping/output/.gitignore | 4 + source/mapping/output/maps/.gitignore | 3 + 19 files changed, 136 insertions(+), 159 deletions(-) rename source/input_data/{MOD11C3_LST => MOD21C3_LST}/.gitkeep (100%) create mode 100644 source/mapping/.gitkeep create mode 100644 source/mapping/data/.gitkeep create mode 100644 source/mapping/output/.gitignore create mode 100644 source/mapping/output/maps/.gitignore diff --git a/.gitignore b/.gitignore index c31647c..d73c5fa 100644 --- a/.gitignore +++ b/.gitignore @@ -179,22 +179,22 @@ pyrightconfig.json .DS_Store # source files -source/input_data/MOD11C3_LST/* +source/input_data/MOD21C3_LST/* source/input_data/MOD13C2_NDVI/* source/input_data/CHIRPS/* source/input_data/soil_moisture/* source/map_data/* source/working_data/* source/output_data/GeoTiffs/* -source/mapping/* +source/mapping/data/* source/mapping/output/maps/* -!source/input_data/MOD11C3_LST/.gitkeep +!source/input_data/MOD21C3_LST/.gitkeep !source/input_data/MOD13C2_NDVI/.gitkeep !source/input_data/CHIRPS/.gitkeep !source/input_data/soil_moisture/.gitkeep !source/map_data/.gitkeep !source/working_data/.gitkeep !source/output_data/GeoTiffs/.gitkeep -!source/mapping/.gitkeep +!source/mapping/data/.gitkeep !source/mapping/output/maps/.gitkeep \ No newline at end of file diff --git a/FinalMapOutput-ESRI.py b/FinalMapOutput-ESRI.py index 64a2894..8b99f10 100644 --- a/FinalMapOutput-ESRI.py +++ b/FinalMapOutput-ESRI.py @@ -2,11 +2,11 @@ import datetime import arcpy ## Change these to match your directory structure. If you are missing directories you will need to create them -map_input_dir = r'D:/ProjectData/CDI/[base directory]/mapping/data/' -map_output_dir = r'D:/ProjectData/CDI/[base directory]/mapping/output/maps/' -tif_dir = r'D:/ProjectData/CDI/[base directory]/output_data/GeoTiffs/' +map_input_dir = r'./source/mapping/data/' +map_output_dir = r'./source/mapping/output/maps/' +tif_dir = r'./source/output_data/GeoTiffs/' reglyr = map_input_dir+'Percentile_colors_new.lyr' -region = '[Country name]' +region = 'eSwatini' ## Resolution can be changed here res = 150 @@ -36,6 +36,6 @@ print 'Creating png map' arcpy.mapping.ExportToPNG(mxd,png_dir+region+'_'+typ+'_'+tif_year+int_month+'.png',resolution=res) print 'Creating pdf map' - arcpy.mapping.ExportToPDF(mxd,pdf_dir+region+'_'+typ+'_'+tif_year+int_month+'.png',resolution=res) + arcpy.mapping.ExportToPDF(mxd,pdf_dir+region+'_'+typ+'_'+tif_year+int_month+'.pdf',resolution=res) print 'Creating jpg map' - arcpy.mapping.ExportToJPEG(mxd,jpg_dir+region+'_'+typ+'_'+tif_year+int_month+'.png',resolution=res) + arcpy.mapping.ExportToJPEG(mxd,jpg_dir+region+'_'+typ+'_'+tif_year+int_month+'.jpg',resolution=res) diff --git a/FinalMapOutput-POR-ESRI.py b/FinalMapOutput-POR-ESRI.py index 35c452e..067c7ea 100644 --- a/FinalMapOutput-POR-ESRI.py +++ b/FinalMapOutput-POR-ESRI.py @@ -2,30 +2,31 @@ import datetime import arcpy ## Change these to match your directory structure. If you are missing directories you will need to create them -map_input_dir = r'D:/Operations/Global_CDI/SA-Zimbabwe/mapping/data/' -map_output_dir = r'D:/Operations/Global_CDI/SA-Zimbabwe/mapping/output/maps/' -tif_dir = "D:\\Operations\\Global_CDI\\SA-Zimbabwe\\output_data\\GeoTiffs\\" +map_input_dir = r'./source/mapping/data/' +map_output_dir = r'./source/mapping/output/maps/' +tif_dir = r'./source/output_data/GeoTiffs/' reglyr = map_input_dir+'Percentile_colors_new.lyr' -region = 'Zimbabwe' +region = 'eSwatini' ## Resolution can be changed here res = 150 +map_remote_output_dir = r'//ndmc-webp01/WebArchive/GlobalCDI/'+region+'/' for typ in ['CDI','LST','NDVI','SPI','SM']: print 'Mapping '+typ arcpy.env.overwriteOutput = True - tif_list = glob.glob(tif_dir+typ+'\\*.tif') + tif_list = glob.glob(tif_dir+typ+'/*.tif') #in_tif = tif_list[-1] for in_tif in tif_list: tif_year = (in_tif[-10:])[0:4] - int_month = ((in_tif[-10:])[0:6])[-2:] - monthinteger = int(int_month) - month = datetime.date(1900, monthinteger, 1).strftime('%b') - if int(tif_year) > 2019: - print "Inside loop" - print in_tif + if tif_year>"2023": + + int_month = ((in_tif[-10:])[0:6])[-2:] + monthinteger = int(int_month) + month = datetime.date(1900, monthinteger, 1).strftime('%b') tif_layer = arcpy.mapping.Layer(in_tif) ## Change this to 10.x.mxd if not running ArcGIS 10.8 or later themxd = map_input_dir+region+'_'+typ+'-10.8.mxd' + #print themxd mxd = arcpy.mapping.MapDocument(themxd) df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0] arcpy.mapping.AddLayer(df, tif_layer,"BOTTOM") @@ -35,13 +36,17 @@ mapdate = (arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT","Date"))[0] mapdate.text = month+' '+tif_year png_dir = map_output_dir+typ+'/png/' - print png_dir pdf_dir = map_output_dir+typ+'/pdf/' jpg_dir = map_output_dir+typ+'/jpg/' - #print 'Creating png map' - print png_dir+region+'_'+typ+'_'+tif_year+int_month+'.png' + rpng_dir = map_remote_output_dir+typ+'/png/' + rpdf_dir = map_remote_output_dir+typ+'/pdf/' + rjpg_dir = map_remote_output_dir+typ+'/jpg/' + print 'Creating png map' arcpy.mapping.ExportToPNG(mxd,png_dir+region+'_'+typ+'_'+tif_year+int_month+'.png',resolution=res) - #print 'Creating pdf map' - arcpy.mapping.ExportToPDF(mxd,pdf_dir+region+'_'+typ+'_'+tif_year+int_month+'.png',resolution=res) - #$print 'Creating jpg map' - arcpy.mapping.ExportToJPEG(mxd,jpg_dir+region+'_'+typ+'_'+tif_year+int_month+'.png',resolution=res) + arcpy.mapping.ExportToPNG(mxd,rpng_dir+region+'_'+typ+'_'+tif_year+int_month+'.png',resolution=res) + print 'Creating pdf map' + arcpy.mapping.ExportToPDF(mxd,pdf_dir+region+'_'+typ+'_'+tif_year+int_month+'.pdf',resolution=res) + arcpy.mapping.ExportToPDF(mxd,rpdf_dir+region+'_'+typ+'_'+tif_year+int_month+'.pdf',resolution=res) + print 'Creating jpg map' + arcpy.mapping.ExportToJPEG(mxd,jpg_dir+region+'_'+typ+'_'+tif_year+int_month+'.jpg',resolution=res) + arcpy.mapping.ExportToJPEG(mxd,rjpg_dir+region+'_'+typ+'_'+tif_year+int_month+'.jpg',resolution=res) diff --git a/STEP_0000_execute_all_steps.py b/STEP_0000_execute_all_steps.py index 6379861..01d62e6 100644 --- a/STEP_0000_execute_all_steps.py +++ b/STEP_0000_execute_all_steps.py @@ -1,4 +1,5 @@ # -*- coding: utf-8 -*- +from netCDF4 import Dataset from STEP_0101_read_hdf_create_LST_anom_netcdf import main as step_0101 from STEP_0102_read_hdf_create_NDVI_anom_netcdf import main as step_0102 from STEP_0103_read_chirps_create_precip_netcdf_and_spi_netcdf import main as step_0103 diff --git a/STEP_0101_read_hdf_create_LST_anom_netcdf.py b/STEP_0101_read_hdf_create_LST_anom_netcdf.py index ff43f26..0521f4f 100644 --- a/STEP_0101_read_hdf_create_LST_anom_netcdf.py +++ b/STEP_0101_read_hdf_create_LST_anom_netcdf.py @@ -48,8 +48,11 @@ def __get_hdf_date(self, file_name): Returns: String of the year/month values in 'YYYYMM' format """ - (year, days) = self.__raw_file_match.match(file_name).groups() - test_date = date(int(year), 1, 1) + timedelta(days=int(days)) + if file_name.find('/') > -1: + file_name = file_name.split('/')[1] + match = self.__raw_file_match.match(file_name) + (year, month) = match.groups() + test_date = date(int(year), int(month), 1) return test_date.strftime("%Y%m") def __get_calendar_value(self, file_name): @@ -61,6 +64,8 @@ def __get_calendar_value(self, file_name): Returns: Integer value of the number of days since Jan 1, 1900 """ + if file_name.find('/') > -1: + file_name = file_name.split('/')[1] origin_date = date(1900, 1, 1) (year, days) = self.__raw_file_match.match(file_name).groups() test_date = date(int(year), 1, 1) + timedelta(days=(int(days) - 1)) @@ -139,14 +144,17 @@ def get_files_to_process(self, all_hdf=False): """ files = [] try: - raw_files = self.__fileHandler.get_raw_file_names('lst_hdf_regex') + raw_files = sorted(self.__fileHandler.get_raw_file_names('lst_hdf_regex'), reverse=True) if all_hdf: # include all HDF files files = raw_files else: # determine which HDF files have not been converted to NetCDF working_files = self.__fileHandler.get_working_file_names('lst_netcdf_regex') # parse out the year/days from the file name # for f in raw_files: - file_date = self.__get_hdf_date(f) + name = f + if name.find('/') > -1: + name = name.split('/')[1] + file_date = self.__get_hdf_date(name) test_file = "STEP_0101_LST_{}_{}.nc".format(self.__region, file_date) if test_file not in working_files: files.append(f) @@ -174,21 +182,16 @@ def create_lst_netcdf_file(self, file_name): # extract SubGrids of the required parameters # with HDFSubGrid(self.__bounds, raw_file_path, self.__hdf_group) as sg: - lst_day = sg.create_sub_grid('LST_Day_CMG') * 0.02 # data is scaled in the HDF file - lst_night = sg.create_sub_grid('LST_Night_CMG') * 0.02 # data is scaled in the HDF file + lst_day = sg.create_sub_grid('LST_Day').astype(float) * 0.02 # data is scaled in the HDF file + lst_night = sg.create_sub_grid('LST_Night').astype(float) * 0.02 # data is scaled in the HDF file qc_day = sg.create_sub_grid('QC_Day') qc_night = sg.create_sub_grid('QC_Night') # compute the LST delta # - qc_day_filter = np.logical_or(np.logical_and(qc_day > 16, qc_day < 64), qc_day > 79) - filtered_lst_day = ma.masked_where(qc_day_filter, lst_day) - - qc_night_filter = np.logical_or(np.logical_and(qc_night > 16, qc_night < 64), qc_night > 95) - filtered_lst_night = ma.masked_where(qc_night_filter, lst_night) - - delta = filtered_lst_day - filtered_lst_night - delta_filter = np.logical_or(filtered_lst_day == 0, filtered_lst_night == 0, delta <= 0) - lst_delta = np.round(ma.masked_where(delta_filter, delta).filled(self.__missing), 3) + filtered_lst_day = ma.masked_where(np.logical_or(qc_day < 16, lst_day == 0), lst_day) + filtered_lst_night = ma.masked_where(np.logical_or(qc_night < 16, lst_night == 0), lst_night) + delta = np.ma.clip(np.ma.subtract(filtered_lst_day, filtered_lst_night), -40.0, 40.0) + lst_delta = np.round(delta.filled(self.__missing), 3) # create the output file # out_properties = { @@ -253,7 +256,7 @@ def update_lst_anomaly_file(self): index = idx # loop thru the years and add the data to the NetCDF file # for y in month_anomalies: - lst_var[index] = y + lst_var[index] = y # np.where(self.__missing, y, np.multiply(y, -1)) # reverse so positive anomalies are dry index += 12 # increment 1 year except IOError: raise diff --git a/STEP_0303_export_ranking_data_rasters.py b/STEP_0303_export_ranking_data_rasters.py index 5837a12..ac8a38f 100644 --- a/STEP_0303_export_ranking_data_rasters.py +++ b/STEP_0303_export_ranking_data_rasters.py @@ -13,46 +13,34 @@ class NetCDFtoTIFF: """ This is the core processing class for executing GeoTiff export """ - def __init__(self, parameter, mode, cdi_date=None): self.__parameter = parameter self.cdi_date = cdi_date self.__config = ConfigParser() self.__mode = mode - self.__cdi_weights = self.__config.get("cdi_parameters", "weights") + self.__cdi_weights = self.__config.get('cdi_parameters', 'weights') def __enter__(self): - if ( - self.__parameter == "cdi" - or self.__cdi_weights[self.__parameter] > 0 - ): - print( - "Exporting Tiff(s) for {}...".format(self.__parameter.upper()) - ) - self.__working_dir = ( - self.__config.get("geotiff_dir").replace("\\", "/") - + "/" - + self.__parameter.upper() - ) - self.__output_dir = self.__config.get("output_dir").replace( - "\\", "/" - ) - self.__region = self.__config.get("region_name") - self.__bounds = self.__config.get("bounds") - self.__latitudes = self.__config.get("latitudes") - self.__longitudes = self.__config.get("longitudes") + if self.__parameter == 'cdi' or self.__cdi_weights[self.__parameter] > 0: + print("Exporting Tiff(s) for {}...".format(self.__parameter.upper())) + self.__working_dir = self.__config.get('geotiff_dir').replace("\\", '/') + '/' + self.__parameter.upper() + self.__output_dir = self.__config.get('output_dir').replace("\\", '/') + self.__region = self.__config.get('region_name') + self.__bounds = self.__config.get('bounds') + self.__latitudes = self.__config.get('latitudes') + self.__longitudes = self.__config.get('longitudes') self.__rows = len(self.__latitudes) self.__cols = len(self.__longitudes) - self.__file_patterns = self.__config.get("file_patterns") + self.__file_patterns = self.__config.get('file_patterns') self.__fileHandler = FileHandler( raw_data_dir=None, working_dir=self.__working_dir, - file_patterns=self.__file_patterns, + file_patterns=self.__file_patterns ) self.__times = None self.__data = None self.__transform = None - self.__projection = "+proj=latlong" + self.__projection = '+proj=latlong' self.__missing = -9999.0 # load the time(s) and data # self.__get_data() @@ -76,42 +64,25 @@ def __get_data(self): input_data_set = None # define the available file names for the sources # input_files = { - "lst": os.path.join( - self.__output_dir, - "STEP_0201_LST_anomaly_pct_rank_{}.nc".format(self.__region), - ), - "ndvi": os.path.join( - self.__output_dir, - "STEP_0202_NDVI_anomaly_pct_rank_{}.nc".format(self.__region), - ), - "spi": os.path.join( - self.__output_dir, - "STEP_0203_SPI_anomaly_pct_rank_{}.nc".format(self.__region), - ), - "sm": os.path.join( - self.__output_dir, - "STEP_0204_SM_pct_rank_{}.nc".format(self.__region), - ), - "cdi": os.path.join( - self.__output_dir, - "STEP_0302_CDI_pct_rank_{}.nc".format(self.__region), - ), + "lst": os.path.join(self.__output_dir, "STEP_0201_LST_anomaly_pct_rank_{}.nc".format(self.__region)), + "ndvi": os.path.join(self.__output_dir, "STEP_0202_NDVI_anomaly_pct_rank_{}.nc".format(self.__region)), + "spi": os.path.join(self.__output_dir, "STEP_0203_SPI_anomaly_pct_rank_{}.nc".format(self.__region)), + "sm": os.path.join(self.__output_dir, "STEP_0204_SM_pct_rank_{}.nc".format(self.__region)), + "cdi": os.path.join(self.__output_dir, "STEP_0302_CDI_pct_rank_{}.nc".format(self.__region)) } # define the NetCDF parameter names for each source # - input_parameters = self.__config.get("cdi_parameters", "names") - input_parameters["cdi"] = "cdi_wt_sum_pr" + input_parameters = self.__config.get('cdi_parameters', 'names') + input_parameters['cdi'] = "cdi_wt_sum_pr" # extract the time(s) and data # try: source = input_files[self.__parameter] source_parameter = input_parameters[self.__parameter] input_data_set = netcdf.open_dataset(source) - if self.__mode == "all": - self.__times = input_data_set.variables["time"][:] - self.__data = netcdf.extract_data( - input_data_set, source_parameter, -1 - ) + if self.__mode == 'all': + self.__times = input_data_set.variables['time'][:] + self.__data = netcdf.extract_data(input_data_set, source_parameter, -1) else: - all_times = input_data_set.variables["time"][:] + all_times = input_data_set.variables['time'][:] last = len(all_times) - 1 if self.cdi_date is not None: for idx, d in enumerate(all_times): @@ -123,9 +94,7 @@ def __get_data(self): self.cdi_date = all_times[last] # extract the data for the last CDI month # self.__times = [all_times[last]] - self.__data = [ - netcdf.extract_data(input_data_set, source_parameter, last) - ] + self.__data = [netcdf.extract_data(input_data_set, source_parameter, last)] except IOError: raise @@ -138,9 +107,7 @@ def __get_data(self): def __get_transformation(self): res = 0.05 try: - self.__transform = Affine.translation( - self.__longitudes[0] - res / 2, self.__latitudes[0] - res / 2 - ) * Affine.scale(res, -res) + self.__transform = Affine.translation(self.__longitudes[0] - res / 2, self.__latitudes[0] - res / 2) * Affine.scale(res, -res) except ValueError: raise except Exception: @@ -158,7 +125,7 @@ def create_date_string(time): """ origin_date = date(1900, 1, 1) data_date = origin_date + timedelta(days=time) - return data_date.strftime("%Y%m") + return data_date.strftime('%Y%m') def __export_geotiffs(self): """ @@ -171,23 +138,18 @@ def __export_geotiffs(self): # loop thru times and generate a GeoTiff for each date # for t, time in enumerate(self.__times): date_str = self.create_date_string(int(time)) - filename = os.path.join( - self.__working_dir, - "STEP_0303_{}_pct_rank_{}_{}.tif".format( - self.__parameter.upper(), self.__region, date_str - ), - ) + filename = os.path.join(self.__working_dir, "STEP_0303_{}_pct_rank_{}_{}.tif".format(self.__parameter.upper(), self.__region, date_str)) # create new GeoTiff # output = rasterio.open( filename, - "w", - driver="GTiff", + 'w', + driver='GTiff', width=self.__cols, height=self.__rows, count=1, dtype=rasterio.float32, crs=self.__projection, - transform=self.__transform, + transform=self.__transform ) # write the data to the image # output.write(self.__data[t].astype(rasterio.float32), 1) @@ -204,9 +166,7 @@ def main(args): """ This is the main entry point for the program """ - # print("PRINT", args) - # mode = str(args.times) - mode = "TEST" + mode = str(args.mode) script_start = datetime.now() try: # set the list of parameters to convert: cdi must be first # @@ -226,14 +186,10 @@ def main(args): print("Script execution: {}".format(script_end - script_start)) -if __name__ == "__main__": +if __name__ == '__main__': # set up the command line argument parser parser = ArgumentParser() - parser.add_argument( - "-t", - "--times", - default="all", - help="The times to export: latest or all. Default is latest", - ) + parser.add_argument("-m", "--mode", default="updates", + help="The times to export: latest or all. Default is updates") # execute the program with the supplied option main(parser.parse_args()) diff --git a/cdi_directory_settings.conf b/cdi_directory_settings.conf index 72bbc2b..a9d4e9b 100644 --- a/cdi_directory_settings.conf +++ b/cdi_directory_settings.conf @@ -1,6 +1,6 @@ { "raw_data_dirs": { - "lst_hdf": "./source/input_data/MOD11C3_LST", + "lst_hdf": "./source/input_data/MOD21C3_LST", "ndvi_hdf": "./source/input_data/MOD13C2_NDVI", "chirps_tif": "./source/input_data/CHIRPS", "fldas_data": "./source/input_data/soil_moisture" diff --git a/cdi_pattern_settings.conf b/cdi_pattern_settings.conf index ed71b4f..4fd0c82 100644 --- a/cdi_pattern_settings.conf +++ b/cdi_pattern_settings.conf @@ -1,6 +1,6 @@ { "file_patterns": { - "lst_hdf_regex": "MOD11C3\\.A((?:19|20)\\d\\d)(\\d\\d\\d)\\S+hdf", + "lst_hdf_regex": "MOD21C3\\.A((?:19|20)\\d\\d)(\\d\\d\\d)\\S+hdf", "ndvi_hdf_regex": "MOD13C2\\.A((?:19|20)\\d\\d)(\\d\\d\\d)\\S+hdf", "chirps_tif_regex": "c((?:19|20)\\d\\d)(0[1-9]|1[0-2])\\.tif", "fldas_data_regex": "FLDAS\\w+\\.A((?:19|20)\\d\\d)(0[1-9]|1[0-2])" diff --git a/cdi_project_settings.conf b/cdi_project_settings.conf index b538dc6..816008a 100644 --- a/cdi_project_settings.conf +++ b/cdi_project_settings.conf @@ -1,7 +1,7 @@ { -"region_name" : "PERU", - "bounds" : {"n_lat": -15.375, "s_lat": -22.675, "w_lon": 24.825, "e_lon": 33.375}, - "spi_periods": [1, 3, 9], +"region_name" : "Eswatini", + "bounds" : {"n_lat": -25.675, "s_lat": -27.825, "w_lon": 30.675, "e_lon": 32.825}, + "spi_periods": [3], "cdi_parameters": { "names": { "lst": "lst_anom_pct_rank", @@ -10,12 +10,12 @@ "sm": "RootZone2_SM_pct_rank" }, "weights" : { - "lst": 0.2, - "ndvi": 0.2, + "lst": 0.3, + "ndvi": 0.3, "spi": 0.4, - "sm": 0.2 + "sm": 0.0 } }, - "map_template": "PERU_template.qpt", - "map_project": "PERU_CDI.qgs" + "map_template": "eswatini_template.qpt", + "map_project": "eswatini_CDI.qgs" } diff --git a/execute_steps.bat b/execute_steps.bat index dd2df47..9533f2b 100644 --- a/execute_steps.bat +++ b/execute_steps.bat @@ -1,12 +1,9 @@ setlocal set PATH=%PATH%;C:\ProgramData\Anaconda3;C:\ProgramData\Anaconda3\scripts - -call activate cdi_37 D: -cd Operations\Global_CDI\[base directory]\scripts +cd Operations\Global_CDI\SA-eSwatini\scripts C:\ProgramData\Anaconda3\envs\cdi_37\python.exe STEP_0000_execute_all_steps.py -call conda deactivate endlocal \ No newline at end of file diff --git a/libs/file_operations.py b/libs/file_operations.py index b898f31..abf0667 100644 --- a/libs/file_operations.py +++ b/libs/file_operations.py @@ -25,23 +25,34 @@ def __validate_directories(self): except IOError: print("Error creating working directory; check permissions on parent directory.") - def get_raw_file_names(self, pattern): + def get_raw_file_names(self, pattern, root=None, subdir=None): """ This function reads the raw data directory to find available files to process Args: pattern (str): name of the file pattern (from the config) to use in the search + root (str): root directory to scan + subdir (str): name of the sub-directory if being used Returns: List of file names """ results = [] + if root is None: + root = self.__raw_data_dir try: # get the list fo files in the raw data directory # - for f in os.listdir(self.__raw_data_dir): + for f in os.listdir(root): + # check if we have a sub-directory # + full_path = os.path.join(root, f) + if os.path.isdir(full_path): + results = results + self.get_raw_file_names(pattern, full_path, f) # test the file name against the pattern # if re.search(r'{}'.format(self.__patterns[pattern]), f): # append the name to the results list # - results.append(str(f)) + if root != self.__raw_data_dir: + results.append('{}/{}'.format(subdir, f)) + else: + results.append(str(f)) except IOError: raise except Exception: diff --git a/libs/hdf_functions.py b/libs/hdf_functions.py index bed2dbb..32c6267 100644 --- a/libs/hdf_functions.py +++ b/libs/hdf_functions.py @@ -35,7 +35,7 @@ def extract_data(data_set, group, parameter, time=0): try: # retrieve the parameter values # data = data_set[group]['Data Fields'][parameter] - results = np.array(data[:]).astype(float) + results = np.array(data[:]) # return the values # if time >= 0: return results[time] diff --git a/libs/statistics_operations.py b/libs/statistics_operations.py index 5785eec..5d362b8 100644 --- a/libs/statistics_operations.py +++ b/libs/statistics_operations.py @@ -30,19 +30,20 @@ def compute_anomalies_from_files(self, files, parameter): data_set = netcdf.open_dataset(f) month_values.append(netcdf.extract_data(data_set, parameter, -1)) data_set.close() - month_values = ma.masked_equal(month_values, self.__missing) # mask out missing data + masked_values = ma.masked_equal(month_values, self.__missing) # mask out missing data + mask = np.where(np.mean(month_values, axis=0) == self.__missing, 1, 0) # compute the mean delta value # - month_mean = ma.average(month_values, axis=0) - masked_mean = ma.masked_equal(month_mean, self.__missing) + month_mean = ma.mean(masked_values, axis=0) + # masked_mean = ma.masked_equal(month_mean, self.__missing) # compute the standard deviation # - month_std = ma.std(month_values, axis=0, ddof=1) - masked_std = ma.masked_equal(month_std, self.__missing) + month_std = ma.std(masked_values, axis=0, ddof=1) + # masked_std = ma.masked_equal(month_std, self.__missing) # compute the anomaly for each month # anomalies = [] - for values in month_values: - month_anomaly = (values-masked_mean)/masked_std - masked_anomaly = month_anomaly.filled(self.__missing) - anomalies.append(masked_anomaly) + for values in masked_values: + month_anomaly = np.ma.true_divide(np.ma.subtract(values, month_mean), month_std) + masked_anomaly = np.ma.masked_where(mask, month_anomaly.filled(0.0)) + anomalies.append(month_anomaly.filled(self.__missing)) return anomalies except ValueError: raise @@ -97,11 +98,11 @@ def rank_parameter(self, values): # create place-holder for the year grids # ranked_data = [] # loop thru the sets and add the compare counts to the masks # - for j, base_data in enumerate(values): + for j, base_data in enumerate(masked_values): strict_ranks = np.zeros_like(base_data) weak_ranks = np.zeros_like(base_data) # loop thru again and compare the remaining sets against the base # - for i, compare_data in enumerate(values): + for i, compare_data in enumerate(masked_values): if i != j: # skip if the same index as the base # # increment the cells based on rank comparison # np.add(strict_ranks, ma.where(base_data > compare_data, 1, 0), out=strict_ranks) @@ -110,14 +111,10 @@ def rank_parameter(self, values): ranks = (strict_ranks + weak_ranks) * 0.5 # append the year to the output # ranked_data.append(ranks) - # mask ranked data # - masked_ranks = ma.array(ranked_data, mask=masked_values.mask) # divide sum by total years # - ranked_data = np.round(masked_ranks / len(values), 3) - # find the absolute maximum rank for scaling - scale = 1 / ma.amax(ma.amax(ranked_data, axis=0)) - # final scaling of the data # - final_ranks = ranked_data * scale + count = np.ma.masked_equal(ma.amax(ranked_data, axis=0) + 1, 0) + pct_data = np.round(np.ma.true_divide(np.ma.asarray(ranked_data), count), 3) + final_ranks = np.ma.masked_equal(pct_data, self.__missing) return final_ranks except ValueError: raise diff --git a/libs/subgrid_calculations.py b/libs/subgrid_calculations.py index 3bda22e..64d6ebe 100644 --- a/libs/subgrid_calculations.py +++ b/libs/subgrid_calculations.py @@ -17,7 +17,7 @@ def __init__(self, aoi, file_path, interpolate=False): self.first_root_y = 0 self.last_root_y = 0 self.columns = int(abs(aoi['e_lon'] - aoi['w_lon']) * 20) + 1 # 0.05 degree spacing = 20 columns/degree - self.rows = int(abs(aoi['n_lat'] - aoi['s_lat']) * 20) + 1 # 0.05 degree spacing = 20 rows/degree + self.rows = int(round(abs(aoi['n_lat'] - aoi['s_lat']) * 20, 0) + 1) # 0.05 degree spacing = 20 rows/degree def __enter__(self): try: diff --git a/source/input_data/MOD11C3_LST/.gitkeep b/source/input_data/MOD21C3_LST/.gitkeep similarity index 100% rename from source/input_data/MOD11C3_LST/.gitkeep rename to source/input_data/MOD21C3_LST/.gitkeep diff --git a/source/mapping/.gitkeep b/source/mapping/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/source/mapping/data/.gitkeep b/source/mapping/data/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/source/mapping/output/.gitignore b/source/mapping/output/.gitignore new file mode 100644 index 0000000..23b4323 --- /dev/null +++ b/source/mapping/output/.gitignore @@ -0,0 +1,4 @@ +# Ignore everything in output/ except .gitignore and maps +* +!.gitignore +!maps \ No newline at end of file diff --git a/source/mapping/output/maps/.gitignore b/source/mapping/output/maps/.gitignore new file mode 100644 index 0000000..2f1d800 --- /dev/null +++ b/source/mapping/output/maps/.gitignore @@ -0,0 +1,3 @@ +# Ignore everything in maps/ except .gitignore +* +!.gitignore \ No newline at end of file