/*============================================================================== VHsuite.cpp: implements the major functions of VH suite package This class contains VIIRS-VH code for the following tasks 1) Merge GVI partial daily maps to a daily map 2) Make weekly composite map from GVI daily maps 3) Calculate ND file from GVI weekly composite map 4) Calculate SM files from a series of ND files 5) Calculate VH file from SM file 6) Create AVHRR climatology from SM files 7) Take sample of VH above files (sample file will be one file even if input is in multiple tiles (ex:500m resolution)) It also cant VH related other tasks, such as: 1) convert global Hydro Estimator data to HDF file 2) convert GVF to grib file example commands: export LD_LIBRARY_PATH="../libs/netcdf-4.3.2/lib/;../libs/hdf5-1.8.13/lib;../libs/hdf-4.2.7/lib" VHP/script/viirsVH_run_weekly.pl week2015001 JPSS1 ./bin/VHsuite.exe control/VHconfig.j01_daily_4km DailyMap j01 2017 330 2017 365 data/DAP1/VHP/bin/VHsuite.exe -S control_vhp/gvi.config.temp DailyMap NH 1989 100 1989 100 ./bin/VHsuite.exe -S control_vhp/VHconfig.NP_climat2016 Create-CLIMAT NP 2000 1 2000 52 ./bin/VHsuite.exe -S control/VHconfig.viirs_daily_500m DailyMap npp 2014 100 2014 100 -tile 19 19 1 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_500m CompositeMap npp 2014 2 2014 2 -tile 19 19 1 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_500m ND npp 2014 2 2014 2 -tile 19 19 1 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_500m SM npp 2014 2 2014 2 -tile 19 19 1 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_500m VH npp 2014 2 2014 2 -tile 19 19 1 ./bin/VHsuite.exe -S control/VHconfig.viirs_daily_4km DailyMap npp 2014 100 2014 100 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_4km MergeDailyMaps npp 2014 2 2014 2 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_4km CompositeMap npp 2014 2 2014 2 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_4km ND npp 2014 2 2014 2 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_4km SM npp 2014 2 2014 2 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_4km VH npp 2014 2 2014 2 ./bin/VHsuite.exe -S control/VHconfig.viirs_weekly_4km_climate Create-CLIMAT npp 2012 45 2012 45 ./bin10/VHsuite.exe -S Configs/VHP500m_weekly_CLIMAT.config Create-CLIMAT npp 2012 45 2012 45 ./bin7/VHsuite.exe -S Configs/VHP500m_weekly.config VH j01 2023 8 2023 8 > ../data/Log/Redo/VH_j01_2023.txt & yearDay=2014170 16km from 500m ./bin/VHsuite.exe SampleFile 32 8 8 "data/VIIRS/500m/daily/${yearDay}/VGVI.G500m.C01.npp.P${yearDay}*" \ data/16km/VIIRS_daily/VGVI.G16.C01.npp.P${yearDay}.nc "data/VIIRS/500m/daily_meta" yearWeek='2023022' satID='j02' ./bin5/VHsuite.exe SampleFile 32 8 8 "data/VIIRSB/500m/weekly/${yearWeek}/VGVI.G500m.C07.${satID}.P${yearWeek}*" \ data/16km/VIIRSB_weekly/VGVI.G16.C07.${satID}.P${yearWeek}.nc ./bin5/VHsuite.exe SampleFile 32 8 8 "data/${satID}/500m/weekly/${yearWeek}/VGVI.G500m.C07.${satID}.P${yearWeek}*" \ data/16km/VIIRS_weekly/VGVI.G16.C07.${satID}.P${yearWeek}.nc ./bin5/VHsuite.exe SampleFile 32 8 8 "data/${satID}/500m/VH/${yearWeek}/VGVI.G500m.C07.${satID}.P${yearWeek}*ND" \ data/16km/VIIRS_weekly/VGVI.G16.C07.${satID}.P${yearWeek}_ND.nc 16km from 4km daily ./bin/VHsuite.exe SampleFile 4 1 1 "data/VIIRS/4km/daily/VGVI_21Bands.G04.C01.npp.P${yearDay}*" \ data/16km/VIIRS_daily21Bands/VGVI_21Bands.G16.C07.npp.P${yearDay}.nc "data/VIIRS/4km/daily_meta" ./bin/VHsuite.exe AggregateFile 4 1 1 "data/VIIRS/4km/daily/VGVI_21Bands.G04.C01.npp.P${yearDay}*" \ data/16km/VIIRS_daily21Bands/VGVI_21Bands.G16.C07.npp.P${yearDay}.nc "data/VIIRS/4km/daily_meta" 1km ND from 500m ./bin/VHsuite.exe AggregateFile 2 8 8 "data/VIIRS/500m/VH/${yearWeek}/VGVI.G500m.C07.npp.P${yearWeek}.ND" \ data/VIIRS/1km/VH/VGVI.G1000m.C07.npp.P${yearWeek}.ND.nc "data/VIIRS/1km/VH_meta" 1km VH from 500m script/viirs_500m_to_1km.sh VH 2017 1 1 Below commands will extact active calibration for AVHRR from 1981 to 2014 ./bin/VHsuite.exe -S control_vhp/VHconfig.NP PrintAvhrrCalibration NC 1981 35 2014 45 -v 3 > \ data/Calibration/avhrr_cali_allyears.txt cat data/Calibration/avhrr_cali_allyears.txt | grep "Active" > data/Calibration/avhrr_cali_allyears_active.txt cat data/Calibration/avhrr_cali_allyears.txt | grep "Note" > data/Calibration/avhrr_cali_allyears_slope.txt // version 1.0 draft, 2011 v1.1, 1.2 2011, AVHRR-VHP output as HDF4 v1.3 April 2012, AVHRR-VHP; output as HDF4 or NETCDF; delivered to OSPO v1.4 Dec 2013, VIIRS-VH, SM output 8 weeks by one step before PDR v1.5 Jan 2015, VIIRS-VH, include procedure for NDE, v1.6 Feb 2015, VIIRS-VH, improvement on text file reader, avoid re-laod AVHRR calibration v1.7 -1.11, July 6 2015 VIIRS-VH modified for NDE requirements. v1.11 7/6/2015 final DAP of VIIRS-VH v1.11.1, 01/04/2015, Add function: (1) add function "AggregateFile" (2) fix bug on hdf5 code for outputing 3D array v2.1 Wenze 2021 : add J01 and J02 add quality_information, id, date_created v3.1 500m VHP: developer: Wei Guo Oct. 2023 - May 2025 - share configure files (one for daily, one for weekly;) - data of all satellites will go to the same folder /data - Daily Map : - Output folder: data/JPSS/500m/daily/yyyyjjj/ - File name: in NDE file name convention (long name, with start, end and creation time) example: data/JPSS/500m/daily/2025080/ VGVI500mDaily_r08c08_v3r1_j02_s202503210000000_e202503212359599_c202503272137494.nc - IP folder: data/JPSS/lists/ ( file lists for input SDR files) - Output minimum variables: use options in daily configure file, skip solar_azimuth, zensor_zenith, MOD band3, cloud Mask - Run once per day: need all granules are available, will take 90 minutes (120 minutes if do not skip azimuth and M3 band) 14GB /day - New function: create daily RGB image by C++ VHsuite.exe, for job monitoring will take 1 minute - Weekly Map: - output file as internal file, (short name) - will search the last daily map file for each tile - Run once per week: will take 60 minutes (100 minutes if do not skip azimuth and M3 band) 16 GB/ day (20 GB if do not skip azimuth and M3 band) - New function: create daily RGB image by C++ VHsuite.exe, for job monitoring will take 1 minute - ND - calculate row NDVI and BT by applying online calibration - output file as NDE file name convention, (long name, with start, end and creation time) - use "LandCover" of each tile, - will take 8 minutes (15 minutes if use 1km IGBP (reduce the reading) ===============================================================================*/ #include "VHsuite.h" #include "lib_grib.h" #include "AvhrrGvi.h" #include "ViirsGvi.h" #include "AdminVH.h" #include "MapTiles.h" #include "VHsuite_extra.cpp" int main(int argc, char * argv[]) { timeval time_begin; gettimeofday( & time_begin, nullptr); struct rusage resourceUsage; char * current_time = ctime( & time_begin.tv_sec); printf("%s start @ %s \n", argv[0], current_time); char inputCommand [LENGTH_1000]=""; for (int i = 0; i < argc; i++){ ADDSTRINGS3(LENGTH_1000, inputCommand, argv[i], " "); } printf("inputCommand= %s\n", inputCommand); char * pPostionArg[LENGTH_SDSNAME]; int N_positionArguments = 0; // below are parameters can be set by command line //perform tasks: char command[LENGTH_FILENAME]; char dir_excutable[LENGTH_FILENAME]; getDirectoryName(argv[0], dir_excutable); // ------------ decode the command line -------------------------------- int i = 1; char file_Configure[LENGTH_FILENAME]; while (i < argc) { char * s = argv[i]; if (strcmp(s, "-v") == 0) { checkNextArgv(s, i + 1, argc); m_global_messageLevel = atoi(argv[i + 1]); i += 2; } else if (strcmp(s, "-h") == 0) { SETSTRING(LENGTH_FILENAME, command, "tail -100 ", dir_excutable, "/../src_VH/README_VHsuite_v3r1.txt"); executeCommand(command); exit(SUCCEED); }else if (strcmp(s, "-S") == 0) { checkNextArgv(s, i + 1, argc); stringCopy(LENGTH_FILENAME, file_Configure, argv[i + 1]); i += 2; } else if (strcmp(s, ">") == 0) { break; //ignore flag stdout re-direct } else { pPostionArg[N_positionArguments] = s; N_positionArguments++; i++; } } //while //printf("After decoing the command: argc=%d, N_positionArguments=%d\n", argc, N_positionArguments); // for (int i = 0; i < N_positionArguments; i++) printf(" pPostionArg[ %d]= %s\n", i, pPostionArg[i]); // -------------------------- decode command ------------------------------ char * action = pPostionArg[0]; char ** p = pPostionArg; checkVerbose(argc, argv); if (strcmp(action, "SampleFile") == 0 || strcmp(action, "AggregateFile") == 0) { int sample = atoi(argv[2]); int tilesX = atoi(argv[3]); int tilesY = atoi(argv[4]); char * inputfilesPattern = argv[5]; char * filenameSampled = argv[6]; char dir_meta[LENGTH_SDSNAME]=""; if(N_positionArguments == 8) stringCopy( LENGTH_SDSNAME, dir_meta, argv[7]); CVHsuite * pVH = new CVHsuite; assert(pVH != nullptr); printf("action=%s sample=%d tilesX=%d tilesY=%d input='%s' output='%s' dir_meta='%s'\n", action, sample, tilesX, tilesY, inputfilesPattern, filenameSampled, dir_meta); int aggregate_method = (strcmp(action, "SampleFile") == 0) ? 0 : 1; pVH -> setDefaultPathsForAncillaryfiles(argv[0]); pVH -> create_sampled_VH_file(action, sample, tilesX, tilesY, inputfilesPattern, filenameSampled, dir_meta, aggregate_method); delete(pVH); }else if (strcmp(action, "Split1kmND") == 0 ) { char * filename = argv[2]; char * outputFolder = argv[3]; CVHsuite * pVH = new CVHsuite; assert(pVH != nullptr); pVH -> split_1km_ND_file(filename,outputFolder ); delete(pVH); }else if (strcmp(action, "RemoveOldFiles") == 0) { char * type = argv[2]; char * dir = argv[3]; int year = atoi(argv[4]); int jday = atoi(argv[5]); int daysToKeep = atoi(argv[6]); char * DIR_Temp = argv[7]; CVHsuite * pVH = new CVHsuite; assert(pVH != nullptr); if (strcmp(type, "daily") == 0) { pVH -> removeOldDailyMaps(dir, year, jday, daysToKeep, DIR_Temp); }else { int week = jday; int weeksToKeep = daysToKeep; pVH -> removeOldWeeklyVHfile(dir, year, week, weeksToKeep, DIR_Temp); } delete(pVH); }else if (strcmp(action, "GVF_nc2grib") == 0) { char ExecutableFolder[LENGTH_FILENAME]; getDirectoryName(argv[0], ExecutableFolder); char * filename = argv[1]; char * outputFilename = argv[2]; createGRIBfile(ExecutableFolder, filename, outputFilename); }else if (strcmp(action, "HE_convertToHDF") == 0) { char * dir = argv[2]; int year = atoi(argv[3]); int day = atoi(argv[4]); HydroEstimator_convert_ASCII_To_HDF(dir, year, day); }else if (strstr(action, "MapTiles") != nullptr) { CMapTiles objMapTile; objMapTile.mapTiles_main(N_positionArguments, pPostionArg); }else if (strstr(action, "adminVH") != nullptr) { CAdminVH objAdminVH; objAdminVH.setDefaultPathsForAncillary(argv[0]); objAdminVH.adminVH_readConfigure(file_Configure) ; objAdminVH.adminVH_main(N_positionArguments, pPostionArg); }else if (strcmp(action, "create_VH_Image") == 0) { //create image from VH file char *filename = p[1]; char *dataType = p[2]; char *filename_png =p[3]; CAdminVH objAdminVH; objAdminVH.setDefaultPathsForAncillary(argv[0]); objAdminVH.create_VH_Image(filename, dataType, filename_png); }else if (strcmp(action, "create_RGB_from500mTiles") == 0) { //create RGB image from VHP daily or weekly composite files // satID, inputFolder, filename_png, [ xTiles, yTiles] // VHsuite.exe create_RGB_from500mTiles j02 data/j02/500m/daily/2025060/ VH_daily_RGB_j02_2025060.png 8 8 // VHsuite.exe create_RGB_from500mTiles j02 data/j02/500m/weekly/2025052/ VH_weekly_RGB_j02_2025052.png char *satID = p[1]; char *folder = p[2]; char *filename_png =p[3]; int xTiles = 8; int yTiles = 8; if(N_positionArguments == 6){ xTiles = atoi(p[4]); yTiles = atoi(p[5]); } CAdminVH objAdminVH; objAdminVH.setDefaultPathsForAncillary(argv[0]); objAdminVH.create_RGB_from500mTiles(satID, folder, xTiles, yTiles,filename_png); }else if (strcmp(action, "Weekly_Job_Status") == 0) { //create RGB image from VHP daily or weekly composite files // VHsuite.exe Weekly_Job_Status 2025 10 data/Log/2025/W202510.log int year = atoi( p[1] ); int week = atoi( p[2] ); char *filename_text =p[3]; CAdminVH objAdminVH; objAdminVH.Monitor_VHPweeklyJobStatus( year, week, filename_text); } else if (strcmp(action, "convert_ViirsVCI_to_AvhrrSMN") == 0) { CVHsuite * pVH = new CVHsuite; assert(pVH != nullptr); char *satID = p[1]; int year = atoi( p[2] ); int week = atoi( p[3] ); char *file_ViirsVH = p[4]; char *file_AvhrrClimat = p[5]; char *file_AVHRR_SM = p[6]; pVH -> convert_ViirsVCI_to_AvhrrSMN(satID, year, week, file_ViirsVH, file_AvhrrClimat,file_AVHRR_SM); delete(pVH); } else if (strcmp(action, "createWeeklyPCF") == 0) { char *satID = pPostionArg[1]; int year = atoi(pPostionArg[2]); int week = atoi(pPostionArg[3]); char *dir_VHP = pPostionArg[4]; /* char *dir_500mDaily = pPostionArg[5]; char *dir_500mND = pPostionArg[6]; char *dir_working = pPostionArg[7]; createWeeklyPCF(satID, year,week,dir_VHP,dir_500mDaily,dir_500mND, dir_working); */ char *dir_500m = pPostionArg[5]; char *dir_working = pPostionArg[6]; createWeeklyPCF(satID, year,week,dir_VHP,dir_500m, dir_working); }else if (strcmp(action, "createDailyPCF") == 0) { //command: VHsuite.exe DailyMap satID year jday //command: VHsuite.exe createJpssDailyPCF satID year jday dir_VHP dailyWorkingFolder dirm_500m // VIIRS only CViirsGvi * objCViirsGvi = new CViirsGvi; assert(objCViirsGvi != nullptr); objCViirsGvi -> createDailyPCF(argc, argv); delete(objCViirsGvi); }else { //VHP tasks for operation printf("===VHP tasks for operation===\n"); CVHsuite * pVH = new CVHsuite; assert(pVH != nullptr); pVH -> processVhJob(argc, argv); delete(pVH); } printPeakMemory( &resourceUsage ); printStatus(inputCommand, "All Done", &time_begin); return (SUCCEED); } /*----------------------------------------------------------------------------- * Construction function for class CVHsuite; * create structure for configure ration * set default values for configuration parameters. -----------------------------------------------------------------------------*/ CVHsuite::CVHsuite() { m_pConfig = (struct_VHsuite_Config * ) new struct_VHsuite_Config; assert(m_pConfig != nullptr); memset(m_pConfig, 0, sizeof(struct_VHsuite_Config)); struct_VHsuite_Config * p = m_pConfig; //yyyymmddHHMMSSS stringCopy(LENGTH_20, p -> job_coverage_start, ""); stringCopy(LENGTH_20, p -> job_coverage_end, ""); stringCopy(LENGTH_FILENAME, p -> ScalingParameters, " 1, 0, 1, 0"); p -> SkipWater = 0; p -> Composite_Method_maxNDVI = 1; p -> PeriodsPerYear = 52; p -> DaysPerPeriod = 7; p -> FilterSize = 15; p -> ND_ProcessWater = 0; p -> output_cell_jday = 0; //update the last 8 weeks p -> smOutputOption = 1; p -> GVI_tiles_X = 1; p -> GVI_tiles_Y = 1; p -> bRunWithPCF= 0; p -> UseNDEfilename= 0; p -> ProcessOceanTile = 1; stringCopy(LENGTH_SDSNAME, p -> Prefix, "VGVI"); stringCopy(LENGTH_SDSNAME, p -> PrefixCLIMAT, "VGVI"); stringCopy(LENGTH_SDSNAME, p -> ResolutionString, "4 km"); p -> Resolution = getResolution(p -> ResolutionString, p -> sResolution_VHP, p -> sResolution_NDE); stringCopy(LENGTH_SDSNAME, p -> SatelliteID, "npp"); stringCopy(LENGTH_SDSNAME, p -> DataSource, "VGVI"); //is not empty only when run at NDE driven by PCF stringCopy(LENGTH_FILENAME, p -> DIR_Working, ""); stringCopy(LENGTH_FILENAME, p -> DIR_VH, "data/VH"); stringCopy(LENGTH_FILENAME, p -> DIR_ND, "data/VH"); stringCopy(LENGTH_FILENAME, p -> DIR_SM, "data/VH"); stringCopy(LENGTH_FILENAME, p -> DIR_PartialDailyMap, "data/dailyPart"); stringCopy(LENGTH_FILENAME, p -> DIR_DailyMap, "data/daily"); stringCopy(LENGTH_FILENAME, p -> DIR_WeeklyComposite, "data/weekly"); stringCopy(LENGTH_FILENAME, p -> DIR_Archive, "data/archive"); stringCopy(LENGTH_FILENAME, p -> DIR_Meta, "data/meta"); stringCopy(LENGTH_FILENAME, p -> DIR_Temp, "./"); stringCopy(LENGTH_FILENAME, p -> DIR_CLIMAT, "data/climate"); //viirsgvi.config.4km stringCopy(LENGTH_FILENAME, p -> File_Configure, ""); stringCopy(LENGTH_FILENAME, p -> File_PCF, ""); stringCopy(LENGTH_SDSNAME, p -> sActiveTile, ""); } /*----------------------------------------------------------------------------- * Deconstruction function for class CVHsuite; * release memory allocated -----------------------------------------------------------------------------*/ CVHsuite::~CVHsuite() { //printf("%s() delete structure 'm_pConfig'\n", __FUNCTION__); delete(m_pConfig); //printf("%s() done\n", __FUNCTION__); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: ReadConfigureFile() * FUNCTION: load configure parameters from a congigure file and save to a structure: * struct_VHsuite_Config m_pConfig * A configure file is an ASCII file which was designed to provide values to * parameters controlling the VHsuite.exe. It follows the followin rules: * - Each line in configure file will set value for a specific parameter, * in the form of "name=value". * - The order of parameter is not significant. * Parameters wll be recognized by name * - Spaces are allowed. * - comments are allowed. * comments are strings leading with "#". * - name of parameters are case sensitive. * Sapce and comments will be removed before read values for parameters * INPUT: * filename_Configure: file name of congfigure file * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::ReadConfigureFile(char * filename_Configure) { if (strlen(filename_Configure) == 0) return 0; printf("\nRead, configure, file=%s\n", filename_Configure); struct_VHsuite_Config * p = m_pConfig; char * s = readTextFileContent(filename_Configure); if (getStringOfTag(s, "working_directory", p -> DIR_Working) == SUCCEED) { checkDirectoryName(p -> DIR_Working, ACCESS_WRITE); char dir[LENGTH_FILENAME]; SETSTRING(LENGTH_FILENAME, dir, p -> DIR_Working, "/", "data"); SETSTRING(LENGTH_FILENAME, p -> DIR_VH, dir, "/", "VH"); SETSTRING(LENGTH_FILENAME, p -> DIR_ND, dir, "/", "VH"); SETSTRING(LENGTH_FILENAME, p -> DIR_SM, dir, "/", "VH"); SETSTRING(LENGTH_FILENAME, p -> DIR_PartialDailyMap, dir, "/", "dailyP"); SETSTRING(LENGTH_FILENAME, p -> DIR_DailyMap, dir, "/", "daily"); SETSTRING(LENGTH_FILENAME, p -> DIR_WeeklyComposite, dir, "/", "weekly"); SETSTRING(LENGTH_FILENAME, p -> DIR_Archive, dir, "/", "archive"); SETSTRING(LENGTH_FILENAME, p -> DIR_Meta, dir, "/", "meta"); } getStringOfTag(s, "File_Metadata_Regions", p -> File_Metadata_Regions); getStringOfTag(s, "File_Static_Metadata", p -> File_Static_Metadata); getStringOfTag(s, "File_IGBPLandTypes", p -> File_IGBPLandTypes); getStringOfTag(s, "File_AST_SurfaceTypes", p -> File_AST_SurfaceTypes); getStringOfTag(s, "File_AST_WaterFraction", p -> File_AST_WaterFraction); getStringOfTag(s, "DIR_GlobalCover", p -> DIR_GlobalCover); getStringOfTag(s, "File_preLaunchCalibration", p -> File_preLaunchCalibration); getStringOfTag(s, "File_postLaunchCalibration", p -> File_postLaunchCalibration); getStringOfTag(s, "File_EDF_NDVI", p -> File_EDF_NDVI); getStringOfTag(s, "File_EDF_BT4", p -> File_EDF_BT4); getStringOfTag(s, "job_coverage_start", p -> job_coverage_start); getStringOfTag(s, "job_coverage_end", p -> job_coverage_end); getStringOfTag(s, "ClimatFileExtension", p -> ClimatFileExtension); getStringOfTag(s, "ScalingParameters", p -> ScalingParameters); getIntegerOfTag(s, "SkipWater", p -> SkipWater); getIntegerOfTag(s, "ND_ProcessWater", p -> ND_ProcessWater); getIntegerOfTag(s, "CompositeMethodMaxNDVI", p -> Composite_Method_maxNDVI); getIntegerOfTag(s, "DaysPerPeriod", p -> DaysPerPeriod); p -> PeriodsPerYear = 366 / p -> DaysPerPeriod; getIntegerOfTag(s, "FilterSize", p -> FilterSize); getIntegerOfTag(s, "smOutputOption", p -> smOutputOption); getIntegerOfTag(s, "NumberOfTiles_X", p -> GVI_tiles_X); getIntegerOfTag(s, "NumberOfTiles_Y", p -> GVI_tiles_Y); getStringOfTag(s, "FilePrefix", p -> Prefix); getStringOfTag(s, "FilePrefixClimat", p -> PrefixCLIMAT); getStringOfTag(s, "Resolution", p -> ResolutionString); p -> Resolution = getResolution(p -> ResolutionString, p -> sResolution_VHP, p -> sResolution_NDE); getStringOfTag(s, "SatelliteID", p -> SatelliteID); getIntegerOfTag(s, "UseNDEfilename", p -> UseNDEfilename); getIntegerOfTag(s, "ProcessOceanTile", p -> ProcessOceanTile); //set defult folders // SETSTRING(LENGTH_FILENAME, p -> DIR_GranuleList, "data/", p -> SatelliteID,"/lists"); // SETSTRING(LENGTH_FILENAME, p -> DIR_DailyMap, "data/", p -> SatelliteID,"/daily"); getStringOfTag(s, "DataSource", p -> DataSource); getStringOfTag(s, "DIR_VH", p -> DIR_VH); getStringOfTag(s, "DIR_ND", p -> DIR_ND); getStringOfTag(s, "DIR_SM", p -> DIR_SM); getStringOfTag(s, "DIR_PartialDailyMap", p -> DIR_PartialDailyMap); getStringOfTag(s, "DIR_DailyMap", p -> DIR_DailyMap); getStringOfTag(s, "DIR_WeeklyComposite", p -> DIR_WeeklyComposite); getStringOfTag(s, "DIR_CLIMAT", p -> DIR_CLIMAT); getStringOfTag(s, "DIR_Archive", p -> DIR_Archive); getStringOfTag(s, "DIR_Meta", p -> DIR_Meta); getStringOfTag(s, "DIR_Temp", p -> DIR_Temp); free(s); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: showParameters() * FUNCTION: show Parameters -----------------------------------------------------------------------------*/ int CVHsuite::showParameters() { printf("%s(): begin\n", __FUNCTION__); struct_VHsuite_Config * p = m_pConfig; char fmtInt[LENGTH_20] = " %-25s= %d\n"; char fmtString[LENGTH_20]= " %-25s= '%s'\n"; printf(fmtInt, "SkipWater", p -> SkipWater); printf(fmtInt, "ND_ProcessWater", p -> ND_ProcessWater); printf(fmtInt, "CompositeMethodMaxNDVI", p -> Composite_Method_maxNDVI); printf(fmtInt, "DaysPerPeriod", p -> DaysPerPeriod); printf(fmtInt, "PeriodsPerYear", p -> PeriodsPerYear); printf(fmtInt, "FilterSize", p -> FilterSize); printf(fmtInt, "smOutputOption", p -> smOutputOption); printf(fmtInt, "NumberOfTiles_X", p -> GVI_tiles_X); printf(fmtInt, "NumberOfTiles_Y", p -> GVI_tiles_Y); printf(fmtInt, "UseNDEfilename", p -> UseNDEfilename); printf(fmtInt, "ProcessOceanTile", p -> ProcessOceanTile); printf(fmtString, "DataSource", p -> DataSource); printf(fmtString, "FilePrefix", p -> Prefix); printf(fmtString, "FilePrefixClimat", p -> PrefixCLIMAT); printf(fmtString, "ResolutionString", p -> ResolutionString); printf(fmtString, "sResolution_VHP", p -> sResolution_VHP); printf(fmtString, "sResolution_NDE", p -> sResolution_NDE); printf(fmtString, "SatelliteID", p -> SatelliteID); printf(fmtString, "DIR_VH", p -> DIR_VH); printf(fmtString, "DIR_ND", p -> DIR_ND); printf(fmtString, "DIR_SM", p -> DIR_SM); //printf(fmtString, "DIR_PartialDailyMap", p -> DIR_PartialDailyMap); printf(fmtString, "DIR_DailyMap", p -> DIR_DailyMap); printf(fmtString, "DIR_WeeklyComposite", p -> DIR_WeeklyComposite); printf(fmtString, "DIR_CLIMAT", p -> DIR_CLIMAT); printf(fmtString, "DIR_Archive", p -> DIR_Archive); //printf(fmtString, "DIR_Meta", p -> DIR_Meta); //printf(fmtString, "DIR_Temp", p -> DIR_Temp); printf("%s(): end\n", __FUNCTION__); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: setDefaultPathsForAncillaryfiles() * FUNCTION: set Default Paths For Ancillary files. * INPUT: * executable_file: the name of executable file (full or relative path), * i.e. the first argument of the command -----------------------------------------------------------------------------*/ void CVHsuite::setDefaultPathsForAncillaryfiles(char * executable_file) { char dir[LENGTH_FILENAME]; getDirectoryName(executable_file, dir); ADDSTRINGS(LENGTH_FILENAME, dir, "/../ancillary"); struct_VHsuite_Config * p = m_pConfig; SETSTRING(LENGTH_FILENAME, p -> File_Metadata_Regions, dir, "/", "regions_for_metadata.txt"); SETSTRING(LENGTH_FILENAME, p -> File_Static_Metadata, dir, "/", "vgvi_static_metadata.txt"); //SETSTRING(LENGTH_FILENAME, p -> File_IGBPLandTypes, dir, "/", "gigbp2_landtype_11.hdf"); //SETSTRING(LENGTH_FILENAME, p ->File_AST_SurfaceTypes, dir, "/", // "VIIRS-AST-IGBP17-GEO_v1r0_multi_s20230101__GVI.nc"); //SETSTRING(LENGTH_FILENAME, p ->File_AST_WaterFraction, dir, "/", "GlobalWaterFraction2020__GVI.nc"); //SETSTRING(LENGTH_FILENAME, p -> DIR_GlobalCover, dir, "/", "GlobalCover"); SETSTRING(LENGTH_FILENAME, p -> File_preLaunchCalibration, dir, "/", "AVHRR_calibration_prelaunch.txt"); SETSTRING(LENGTH_FILENAME, p -> File_postLaunchCalibration, dir, "/", "AVHRR_calibration_postlaunch.txt"); SETSTRING(LENGTH_FILENAME, p -> File_EDF_NDVI, dir, "/", "NVI_counts_ByLine_G04.hdf"); SETSTRING(LENGTH_FILENAME, p -> File_EDF_BT4, dir, "/", "BT4_counts_ByLine_G04.hdf"); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: processVhJob() * FUNCTION: processing a VH job. * PROCEDURE: * - decode command line, find the task name and start and end week (or jday), * - read configure file * - load ancillary data (if its required) such as * regions for metadata * Periods for GVI ( which satellite will be used for each period) * Period for Climatology (which satellite will be used for each period) * - determine the GVI tile to be processed * - loop for each week and each tile to process one of the following jobs: * 1) Merge GVI partial daily maps to a daily map * 2) Make weekly composite map from GVI daily maps * 3) Calculate ND file from GVI weekly composite map * 4) Calculate SM files from a series of ND files * 5) Calculate VH file from SM file * 6) Create AVHRR climatology from SM files * INPUT: * argc: the number of command arguments * argv: the strings of command arguments * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::processVhJob(int argc, char * argv[]) { timeval time_begin; gettimeofday( & time_begin, nullptr); //arguments will be: // action satelliteID [year1 week1 year2 week2] [-tile tileS tileE tileInterval] // [-S configure_filename] [-PCF PCF_filename] [-v verbose_level] //file_PCF will be common configure file, //we will loop for each week and tile, search the input files, create PCF, then run it char * pPostionArg[20]; int tileRange[3] = { 1, 9999, 1 }; int orbitRange[3]= { -1, -1}; struct_VHsuite_Config * p = m_pConfig; //WG2025 int cmd_SMoutput= -1; int N_positionArguments = vh_decodeCommand(argc, argv, p -> File_Configure, p -> File_PCF, tileRange, orbitRange, cmd_SMoutput, pPostionArg); char * action = pPostionArg[0]; if (strcmp(action, "DailyMap") == 0 ) { //command: VHsuite.exe DailyMap satID year jday int isVIIRS= (N_positionArguments==1 ) ? 1 : isVIIRS_satellite( pPostionArg[1]); if (isVIIRS ==1) { CViirsGvi * objCViirsGvi = new CViirsGvi; assert(objCViirsGvi != nullptr); objCViirsGvi -> createDailyMaps(argc, argv); delete(objCViirsGvi); } else { CAvhrrGvi * objCAvhrrGvi = new CAvhrrGvi; assert(objCAvhrrGvi != nullptr); objCAvhrrGvi -> createDailyMaps(argc, argv); delete(objCAvhrrGvi); } return (SUCCEED); } if (!((strcmp(action, "PrintAvhrrCalibration")) == 0 || (strcmp(action, "WeeklyTasks")) == 0 || (strcmp(action, "ND")) == 0 || (strcmp(action, "SM") == 0) || (strcmp(action, "VH") == 0) || (strcmp(action, "CompositeMap") == 0) || (strcmp(action, "Create-CLIMAT") == 0) )) { return (FAIL); } setDefaultPathsForAncillaryfiles(argv[0]); ReadConfigureFile(p -> File_Configure); ReadConfigureFile(p -> File_PCF); if (strcmp(p -> File_PCF, "") != 0) p -> bRunWithPCF = true; //WG2025 // if "-SMoutput" appeared in command line, use the command value for smOutputOption // i.e, override the default or 'smOutputOption' setting in configure file if(cmd_SMoutput >=0 ) p-> smOutputOption = cmd_SMoutput; showParameters(); //read data list used for GVI data // m_pGVI_dataPeriods = (Periods * ) new Periods; assert(m_pGVI_dataPeriods != nullptr); m_pGVI_dataPeriods -> yearPeriod1[0] = -1; int NsatPeriod = readDataPeriodsFromConfigureFile(m_pConfig -> File_Configure, "[Periods of GVI data used for VH]", m_pGVI_dataPeriods); m_pConfig -> SatellitePeriodsDefined = (NsatPeriod >0)? 1 :0; printf(" NsatPeriod=%d, ==> SatellitePeriodsDefined=%d\n", NsatPeriod,m_pConfig -> SatellitePeriodsDefined); //get satelliteID from command line int N = N_positionArguments; if (N >= 2) stringCopy(LENGTH_SDSNAME, m_pConfig -> SatelliteID, pPostionArg[1]); //get the range of tiles to processed int maxTiles = m_pConfig -> GVI_tiles_X * m_pConfig -> GVI_tiles_Y; if (tileRange[1] > maxTiles) tileRange[1] = maxTiles; //get the period to processing //for merge daily maps: ------------------------------------------------------- //this is number of command positonal arguments required for setting stat and end week (or jday) int NN = 6; { //for weekly tasks: ------------------------------------------------- int year1, week1, year2, week2; if (N >= NN) { year1 = atoi(pPostionArg[2]); week1 = atoi(pPostionArg[3]); year2 = atoi(pPostionArg[4]); week2 = atoi(pPostionArg[5]); } else { // week range is not provided, make it work for latest week getPreviousPeriod(m_pConfig -> DaysPerPeriod, year2, week2); year1 = year2; week1 = week2; if (strcmp(action, "VH") == 0) { week1 -= 7; adjustPeriod(year1, week1); } } printf(" action=%s, year1=%4d,week1=%3d year2=%4d week2=%4d tile=%4d %4d %4d\n", action, year1, week1, year2, week2, tileRange[0], tileRange[1], tileRange[2]); int year = year1; int week = week1; while (1) { processTiles4WeeklyTask(tileRange, action, year, week); if ((year * 1000 + week) >= (year2 * 1000 + week2)) break; week++; adjustPeriod(year, week); } //while } //end if for weekly tasks-------------------------------------------------- delete m_pGVI_dataPeriods; return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: processTiles4WeeklyTask() * FUNCTION: process Tiles for a Weekly task * PROCEDURE: * - loop each tile, perform the required task * INPUT: * tileRange: begin and end tile number * action: task to perform * year, week: the year and week number of the week being processed * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::processTiles4WeeklyTask(int * tileRange, const char * action, int year, int week) { int maxTiles = m_pConfig -> GVI_tiles_X * m_pConfig -> GVI_tiles_Y; printf("\n----------------------------------------------------------\n"); printf(" %s, ", __FUNCTION__); if (maxTiles > 1) printf(" tileRange=[%d,%d], ", tileRange[0], tileRange[1]); printf(" %s year=%d week=%02d\n", action, year, week); printf("----------------------------------------------------------\n"); for (int tile = tileRange[0]; tile <= tileRange[1]; tile += tileRange[2]) { setTagForActiveTile(tile, m_pConfig -> GVI_tiles_X, m_pConfig -> GVI_tiles_Y, m_pConfig -> sActiveTile); if( VHP_SkipOceanTile(m_pConfig -> sActiveTile+1, m_pConfig -> ProcessOceanTile) == true) continue; char task[LENGTH_FILENAME]; if (strcmp(action, "Create-CLIMAT") == 0) { snprintf(task, LENGTH_FILENAME, "Create '%s' file for %d week %d %s", action, year, week, m_pConfig -> sActiveTile); } else { snprintf(task, LENGTH_FILENAME, "Create '%s' file for week %d, tile='%s'", action, week, m_pConfig -> sActiveTile); } //if (maxTiles != 1) printf("task: %s, action:%s\n", task,action); if (strcmp(action, "MergeDailyMaps") == 0) { int jday1 = (week - 1) * m_pConfig -> DaysPerPeriod + 1; int jday2 = (week) * m_pConfig -> DaysPerPeriod; for (int jday = jday1; jday <= jday2; jday++) mergeDailyMaps(year, jday); } if (strcmp(action, "CompositeMap") == 0) composeWeeklyMap(task, year, week); if (strcmp(action, "ND") == 0) { int isVIIRS = isVIIRS_satellite(m_pConfig -> SatelliteID); if(isVIIRS ==1 ) create_ND_File_ForVIIRS(task, year, week); else create_ND_File_ForAVHRR(task, year, week); } if (strcmp(action, "SM") == 0) { printf("calling: create_SM_Files() ...\n"); create_SM_Files(task, year, week); } if (strcmp(action, "VH") == 0) create_VHI_Files(task, year, week); if (strcmp(action, "Create-CLIMAT") == 0) createClimatologyFiles(task, week); if (strcmp(action, "PrintAvhrrCalibration") == 0) { int jday = (week - 0.5) * m_pConfig -> DaysPerPeriod; loadCalibrationsForAVHRR(year, jday, & m_sAvhrrCalibration); } } //tile return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: create_ND_File_ForVIIRS() * FUNCTION: compute NDVI from weekly GVI for VIIRS only * PROCEDURE: * - search input file names (weekly composite file), * they can be provided by PCF file; if not, search specified folder to * get file with specified satelliteID, file type, year and week * - set buffers for intput files to read (setFileToRead) * - set for the output ND file (setOutputFile) * - create output file, * - create SDS in the output file * - write file attributes * - initialize buffer for metadata * - intiialieze buffers for SDS arries associated to input and output files * - read reflectance of VIIRS image band 1 and 2, BT of VIIRS image band 5 * - calculate NDVI and BT by apply calibration and BT correction if it is * required * - set QA flag based * set bit for dersert by IGBP landtype * set bit for landsea by landsea mask stored in weekly composite map * set bit for 'too cold' flag by BT * set bit for general quality based on value of NDVI and BT * - save to ND file * - calculate and save region averaged NDVI and BT to metadata file * - add result file names to PSF file * INPUT: * year, week: the year and week number of the week being processed * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::create_ND_File_ForVIIRS(char * task, int year, int week) { char inputVariables[LENGTH_FILENAME] = "reflectance_I1;reflectance_I2;temperature_I5;packed_cloud_mask;"; timeval time_begin; gettimeofday( & time_begin, nullptr); //input files: weekly composite of radiance count CVhFile fileCH; //output file will contain NDVI and BT4 CVhFile fileND; //will read SDS as a whole array int input_segments = 1; //open input weekly composite file, set the SDS to read, //read file attribute and SDS attributes, then close the file if (setInputVHfile( & fileCH, FILE_GVI, year, week, inputVariables, input_segments) == FAIL) { printf(" ND file was not produced for this week [%d week %d] because weekly composite file was not found.\n", year, week); printStatus(task, WARNING_ND_NO_WEEKLY_MAP); return FAIL; } enum FILE_SDSID_CH { ID_CH1, ID_CH2, ID_CH4, ID_CM }; //open output NDVI file, set the SDS to write, create file, //write file attribute and SDS attributes, then close the file fileND.InputFileList_initialize(LENGTH_FILENAME); fileND.InputFileList_add(fileCH.filename); if (setOutputVHfile( & fileND, FILE_ND, year, week, "NDVI;BT4;QA;", & fileCH) == FAIL) return FAIL; int NX= fileCH.m_Segment.NX; int NY= fileCH.m_Segment.ValidLines; long NP = fileCH.m_Segment.ValidPoints; //required buffers: // weekly map: 3 x int16 +1 x uint8 // ND : 2 x int16 +1 x uint8 float memory_required = NP / 1.0E6 * (7 + 5); printf(" Memory_required=%f MB for processing 'ND'\n", memory_required); enum FILE_SDSID_RAW { ID_NDVI, ID_BT4, ID_NDQA }; fileCH.processBuffer(ALLOCATE); fileND.processBuffer(ALLOCATE); //open file, read data arrays, then close the file; fileCH.processBuffer(OPEN_AND_READ); fileND.processBuffer(INITIALIZE); // for 'ND' procedure, we have 3 options to input information for 'desert/land/coast" // option1 (v2013) : if m_pConfig -> File_IGBPLandTypes is set, // read while array of 1km IGBPLandtype (gigbp2_landtype_11.hdf) // option2 (v2017): if(m_pConfig -> DIR_GlobalCover is set, // for 500m tiles, read "Global Cover" data for each tile ( saved in 64 files) // option3 (v2025): if(m_pConfig -> File_AST_SurfaceTypes and m_pConfig -> File_AST_WaterFraction is set // VIIRS_AST WaterFraction and //calcualte NDVI and BT4 ; uint8 * pNDQA = (uint8 * ) fileND.Buffers[ID_NDQA]; float NDVI, rawBT; uint8 * pIsLand = (uint8 * ) ALLOCATE(NP *2); uint8 * pDesert = pIsLand + NP; VHP_load_Land_desert( m_pConfig -> DIR_GlobalCover, m_pConfig -> File_AST_WaterFraction, m_pConfig -> File_AST_SurfaceTypes, m_pConfig -> sResolution_NDE, m_pConfig -> sActiveTile, NP, pIsLand,pDesert); for (long i = 0; i < NP; i++) { int isLand = pIsLand[i]; int isDesert= pDesert[i] ; //process land pixel only NDVI = -1; rawBT = -1; if (m_pConfig -> ND_ProcessWater ==1 || isLand == 1) { float R1 = fileCH.getValue(ID_CH1, i); float R2 = fileCH.getValue(ID_CH2, i); if (R1 > 0 && R2 > 0) { NDVI = (R2 - R1) / (R2 + R1); fileND.setValue(ID_NDVI, i, NDVI); } rawBT = fileCH.getValue(ID_CH4, i); if (rawBT > MISSING_FLOAT) { fileND.setValue(ID_BT4, i, rawBT); } } //set QA flags char qualityFlag = 0; //if both failed, set bit 0 as invalid (1) if (NDVI <= -1 && rawBT <= 0) qualityFlag |= 1; //set bit 2 as (1) for land pixel if (isLand == 1) qualityFlag |= 4; //set desert flag, bit 1 if (isDesert == 1) qualityFlag |= 2; //set bit 4: too cold if (rawBT > 0 && rawBT < 273.15) qualityFlag |= 16; pNDQA[i] = qualityFlag; } free(pIsLand); Box box; getYrangeWithLandPixels(pNDQA, NX, NY, m_pConfig -> sActiveTile, &box); fileND.checkInt16Outlier(pNDQA); fileND.saveVH_2NC( pNDQA ); fileCH.processBuffer(DEALLOCATE); fileND.processBuffer(DEALLOCATE); //printStatus(task, COMPLETED_NORMALLY, & time_begin); return (SUCCEED); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: create_ND_File_ForAVHRR() * FUNCTION: process NDVI and BT4 for a particular week for AVHRR only * PROCEDURE: * - search input file names (weekly composite file), * they can be provided by PCF file; if not, search specified folder to * get file with specified satelliteID, file type, year and week * - set buffers for intput files to read (setFileToRead) * - set for the output ND file (setOutputFile) * - create output file, * - create SDS in the output file * - write file attributes * - initialize buffer for metadata * - intiialieze buffers for SDS arries associated to input and output files * - load prelaunch calibration * - load postlaunch calibration * intepolate postlunch calibration to given week * - read reflectance of AVHRR band 1 and 2, BT of AVHRR band 4 * - calculate NDVI and BT by apply calibration and BT correction * - apply EDF correction if it is required * - set QA flag based * set bit for dersert by IGBP landtype * set bit for landsea by landsea mask stored in weekly composite map * set bit for 'too cold' flag by BT * set bit for general quality based on value of NDVI and BT * - save to ND file * - add result file names to PSF file * INPUT: * year, week: the week to be processed. -----------------------------------------------------------------------------*/ int CVHsuite::create_ND_File_ForAVHRR(char * task, int year, int week) { int jday = (week - 0.5) * m_pConfig -> DaysPerPeriod; int satNumber = loadCalibrationsForAVHRR(year, jday, & m_sAvhrrCalibration); char inputVariables[LENGTH_FILENAME] = "ch1_count;ch2_count;ch4_temperature;packed_cloud_mask;"; timeval time_begin; gettimeofday( & time_begin, nullptr); // read PCF file, and put file names in an array of strings . //input files: weekly composite of radiance count CVhFile fileCH; //output file will contain NDVI and BT4 CVhFile fileND; //will read SDS as a whole array int input_segments = 1; //open input weekly composite file, set the SDS to read, //read file attribute and SDS attributes, then close the file if (setInputVHfile( & fileCH, FILE_GVI, year, week, inputVariables, input_segments) == FAIL) return FAIL; enum FILE_SDSID_CH { ID_CH1, ID_CH2, ID_CH4, ID_CM }; //open output NDVI file, set the SDS to write, create file, //write file attribute and SDS attributes, then close the file fileND.InputFileList_initialize(LENGTH_FILENAME); fileND.InputFileList_add(fileCH.filename); if (setOutputVHfile( & fileND, FILE_ND, year, week, "NDVI;BT4;QA;", & fileCH) == FAIL) return FAIL; long NP = fileCH.m_Segment.ValidPoints; //required buffers: // weekly map: 3 x int16 +1 x uint8 // ND : 2 x int16 +1 x uint8 float memory_required = NP / 1.0E6 * (7 + 5); printf(" Memory_required=%f MB for processing 'ND'\n", memory_required); enum FILE_SDSID_RAW { ID_NDVI, ID_BT4, ID_NDQA }; fileCH.processBuffer(ALLOCATE); fileND.processBuffer(ALLOCATE); //open file, read data arrays, then close the file; fileCH.processBuffer(OPEN_AND_READ); fileND.processBuffer(INITIALIZE); fileND.loadIGBPLandTypeForOutputBuffer(m_pConfig -> File_IGBPLandTypes); //calcualte NDVI and BT4 ; uint8 * pcm = (uint8 * ) fileCH.Buffers[ID_CM]; uint8 * pNDQA = (uint8 * ) fileND.Buffers[ID_NDQA]; float NDVI, rawBT; for (long i = 0; i < NP; i++) { float count_CH1 = fileCH.getValue(ID_CH1, i); float count_CH2 = fileCH.getValue(ID_CH2, i); if (count_CH1 > 0 && count_CH2 > MISSING_FLOAT) { NDVI = m_sAvhrrCalibration.calculate_NDVI_For_AVHRR(count_CH1, count_CH2); fileND.setValue(ID_NDVI, i, NDVI); } rawBT = fileCH.getValue(ID_CH4, i); if (rawBT > MISSING_FLOAT) { float BT = m_sAvhrrCalibration.calculate_BT_For_AVHRR(satNumber, rawBT); fileND.setValue(ID_BT4, i, BT); } } //EDF correction for raw ND and BT, if need if ((year < 1989) || (year >= 1991 && year <= 1994) || (year >= 1999 && year <= 2000)) { EDF_CorrectLineByLine( & fileND, "NVI", year, week); EDF_CorrectLineByLine( & fileND, "BT4", year, week); } //set QA flags for (long i = 0; i < NP; i++) { NDVI = fileND.getValue(ID_NDVI, i); rawBT = fileND.getValue(ID_BT4, i); char landtype = fileND.m_pIGBPLandtype[i]; char qualityFlag = 0; //if both failed, set bit 0 as invalid (1) if (NDVI <= -1 && rawBT <= 0) qualityFlag |= 1; char landsea = pcm[i] & 12; //copy land and coast bits ( bit 2 and 3) qualityFlag |= landsea; // if(landsea ==0 ) qualityFlag |= 1; ; //for water, set as invalid //set desert flag, bit 1 if (landtype == 16 || landtype == 7) qualityFlag |= 2; //set bit 4: too cold if (rawBT > 0 && rawBT < 273.15) qualityFlag |= 16; pNDQA[i] = qualityFlag; } fileND.checkInt16Outlier(pcm); fileND.saveVH_2NC( pNDQA ); fileCH.processBuffer(DEALLOCATE); fileND.processBuffer(DEALLOCATE); //printStatus(task, COMPLETED_NORMALLY, & time_begin); return (SUCCEED); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: create_SM_Files() * FUNCTION: smooth the NVI or BT4 time series of N weeks. * The length of time series (or called filter size) is 15 weeks. * (it should be >15 and <52). * The output SDS are "SMN" or "SMT". multiple output SM files may be * re-written based on output option, by default, all SM files in the * time series will be updated except the first 7 weeks. * smoothing filter will be applied on each pixel. * We will use small buffers for input files (1/16 of full SDS size), * but the output buffer will be full SDS. * PROCEDURE: * - determine the start and end week number of time series by given * reference week * if the reference week is near start or end of dataset record, * adjument will be applied so that time series have enough valid * data file. for example, * the last week with available data is 2014 week 20, * filter size=15 and * reference week is 2014 week 18, * then, the time series will be * from [2014, 6] to [2014, 20] * If the reference week is in the middle of data set, the time * series will include the same number of weeks befor and after * the reference week, for example * the last week with available data is 2014 week 20, * filter size=15 and * reference week is 2014 week 6, * then, the time series will be * from [2013, week 51] to [2014, week 13] * - determine which SM file should be output * based on the optional parameter: smOutputOption * = 0 output SM file for the reference week * = 1 output SM files except the first 7 weeks. * = 2 output SM files except the last 7 weeks. * = 3 output all SM files * - search input file name (ND file) for each week in the time series, * if PCF is specified, we will read ND file names from PCF, * otherwise search files in specified folder to get file with * specified satelliteID, file type, year and week. * - check if there are enough input ND files * - set buffers for intput files to read (setFileToRead) * - set for the output SM files (setOutputFile) * - create output file, * - create SDS in the output file * - write file attributes * - initialize buffer for metadata * - intiialieze buffers for SDS arries associated to input and output files * - read QA flag (full SDS) from the first valid ND file * we will use the land and coast information * - loop for each input segment for input file * read ND and BT * apply smoothing filter for each pixel * set SMN and SMT to output arries * - set QA flag in SM file * copy bit 1,2,3 (desert,land,coast) from QA of fist input ND file * set bit for 'too cold' flag by BT * set bit for general quality based on value of SMN and SMt * - save to SM files * - add result file names to PSF file * RAM required: * It is determined by how many weeks' SM file will be output. * (The Input SDSs will be processed by segments and RAM required is . * for 4KM product, fileter size=15, input segement= 16), output last 8 weeks, * RAM required = 36MB (15 /16*2 +1 + 8 *2) ~= 700 MB * INPUT: * year, week: the reference week to be output * year2, week2: the last week of the time series * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::create_SM_Files(const char * task, int year, int week) { printf("%s(), year=%d, week=%d, begin......\n", __FUNCTION__, year, week); //timeval time_begin; //gettimeofday( & time_begin, nullptr); int input_segments = 1; //begin and ending week of the time series int year1, week1, year2, week2; //determine the time series coverage by configure file option printf("%s(), year=%d, week=%d, begin...\n", __FUNCTION__, year, week); getInputWeeksForSmooth(year, week, year1, week1, year2, week2); printf("%s(), year=%d, week=%d ==> year1= %d week1= %d, year2= %d week1= %d\n", __FUNCTION__, year, week, year1, week1, year2, week2); int FilterSize = m_pConfig -> FilterSize; //int DaysPerPeriod = m_pConfig->DaysPerPeriod; int smOutputOption = m_pConfig -> smOutputOption; printf("fileRaw: DigitalFilterMaxSize=%d\n", DigitalFilterMaxSize); char filenames_ND[DigitalFilterMaxSize][LENGTH_FILENAME]; // read PCF file, and put file names in an array of strings . //input files CVhFile fileRaw[DigitalFilterMaxSize]; //QA the first available input file CVhFile fileRawQA; //output file for the week being processed, CVhFile fileSM[DigitalFilterMaxSize]; year1 = year2; week1 = week2 - (FilterSize - 1); adjustPeriod(year1, week1); printf(" %s using: %04d week %02d, to %04d weeek %02d, FilterSize=%d\n", __FUNCTION__, year1, week1, year2, week2, FilterSize); printf(" FilterSize=%d\n", FilterSize); //for each week, search available ND file in the list, if found at K , set ND_fileAvailables[i]=K otherwise = -1 printf(" Check available ND files int this time series, \n"); printf(" If file exists, if yes, read file attributes (latlonRange,size of array):\n"); int timeSeries_years[DigitalFilterMaxSize]; int timeSeries_weeks[DigitalFilterMaxSize]; int inputFileStatus[DigitalFilterMaxSize]; int firstValidfile = -1; int NumberOfValidNDfiles = 0; //---CHECK inpput ND files ------------------------------------------------------------------- int missingFiles_uninterrupted = 0; int max_missingFiles_uninterrupted = 0; int yearTemp = year1; int weekTemp = week1; //check input ND files for (int i = 0; i < FilterSize; i++) { timeSeries_years[i] = yearTemp; timeSeries_weeks[i] = weekTemp; inputFileStatus[i] = FAIL; //use the internal file name convention ; stringCopy(LENGTH_FILENAME, filenames_ND[i], ""); getFilename_VHPstyle(filenames_ND[i], FILE_ND, yearTemp, weekTemp); inputFileStatus[i] = FAIL; //read ND file information, get array size, geo-location et al. if (FileExist(filenames_ND[i]) == true) { inputFileStatus[i] = fileRaw[i].setFileToRead(filenames_ND[i], "ND", yearTemp, weekTemp, "NDVI;", input_segments); } //counting: // (1) NumberOfValidNDfiles: number of files, // (2) max_missingFiles_uninterrupted: maximum of number of missing Files without uninterrupted // (3) the first valid file if (inputFileStatus[i] == FAIL) { missingFiles_uninterrupted++; if (max_missingFiles_uninterrupted < missingFiles_uninterrupted) max_missingFiles_uninterrupted = missingFiles_uninterrupted; } else { missingFiles_uninterrupted = 0; NumberOfValidNDfiles++; if (firstValidfile == -1) firstValidfile = i; } printf(" filenames_ND[%02d]=%s ", i, filenames_ND[i]); if (inputFileStatus[i] == FAIL) printf("FAIL \n"); else printf("OK \n"); weekTemp++; adjustPeriod(yearTemp, weekTemp); } //check if there are enough valid input files for smoothing int missingFiles = FilterSize - NumberOfValidNDfiles; printf(" Summary on input ND files:\n"); printf(" NumberOfValidNDfiles= %d\n", NumberOfValidNDfiles); printf(" missingFiles= %d\n", missingFiles); printf(" max_missingFiles_uninterrupted= %d\n", max_missingFiles_uninterrupted); if (max_missingFiles_uninterrupted > 3) { printf(" There are %d uninterrupted missing files in input time series (>3)\n", max_missingFiles_uninterrupted); printStatus(task, WARNING_COR_TOO_MUCH_UNITERRUPTED_MISSING_FILES); return FAIL; } if (missingFiles > 4) { printf(" There are %d missing files in input time series(>4)\n", missingFiles); printStatus(task, WARNING_COR_TOO_MUCH_MISSING_FILES); return FAIL; } //---CHECK inpput ND files ------------------------------------------------------------------- //determine which weeks' SM files are required to output int SmFilesNeedOutput[DigitalFilterMaxSize]; int SM_outputFiles = 0; int firstOutputFile = -1; int32 i, ix; float fmissing = MISSING_FLOAT; //this control the RAM use //due to limitation of RAM, we process input files by segments // but the output SDS will be writen as a whole array to utilize the compression feture of HDF. char DirectoryName[LENGTH_FILENAME]; getDirectoryName(fileRaw[firstValidfile].filename, DirectoryName); struct_Segment * pSegment = & fileRaw[firstValidfile].m_Segment; int segment_Points = pSegment -> NX * pSegment -> NY; //initialize smoothing filter: float filter_weight = 0.2; float MedianFilterSize = 5; // allocate 3 buffers with size FilterSize, for single point only float * filter_pTS = (float * ) ALLOCATE(FilterSize * 3 * 4); float * filter_pA = filter_pTS + FilterSize; float * filter_pTemp = filter_pTS + FilterSize * 2; initializeFilter(FilterSize, filter_weight, filter_pA); //create output SM files printf(" Create output SM files:\n"); printf(" smOutputOption=%d : ", smOutputOption); if (smOutputOption == 0) printf(" output one week only\n"); //No-Begin weeks if (smOutputOption == 1) printf(" output all weeks except the first 7 weeks\n"); //No-Tail Weeks if (smOutputOption == 2) printf(" output all weeks except the last 7 weeks\n"); //All weeks if (smOutputOption == 3) printf(" output all weeks \n"); //No Begin and Tail Weeks if (smOutputOption == 4) printf(" output all weeks except the first and last 7 weeks\n"); //only the initial and center week if (smOutputOption == 5) printf(" output only the initial and center week\n"); //open SM file if it is required for (i = 0; i < FilterSize; i++) { //check if SM file need to write SmFilesNeedOutput[i] = 0; if ( //specific week only (smOutputOption == 0 && timeSeries_years[i] == year && timeSeries_weeks[i] == week) || //No-Begin weeks (smOutputOption == 1 && (i >= 7)) || //No-Tail Weeks (smOutputOption == 2 && (i < (FilterSize - 7))) || //All weeks (smOutputOption == 3) || //No Begin and Tail Weeks (smOutputOption == 4 && (i >= 7) && (i < (FilterSize - 7))) || //only the initial and center week (smOutputOption == 5 && (i == (FilterSize / 2) || i == (FilterSize - 1)))) { SmFilesNeedOutput[i] = 1; SM_outputFiles++; if (firstOutputFile == -1) firstOutputFile = i; //save the list of input file names to meta data structure of output file ; fileSM[i].InputFileList_initialize(LENGTH_FILENAME * FilterSize); snprintf(fileSM[i].m_Meta.pInputFilenames, LENGTH_FILENAME * FilterSize, "Directory='%s'\n\nInput files:\n", DirectoryName); for (int k = 0; k < FilterSize; k++) { if (inputFileStatus[k] == SUCCEED) fileSM[i].InputFileList_add(getFileBaseName(fileRaw[k].filename)); } if (setOutputVHfile( & fileSM[i], FILE_SM, timeSeries_years[i], timeSeries_weeks[i], "SMN;SMT;QA;", & fileRaw[firstValidfile]) == FAIL) return FAIL; printf(" filenames_SM[%02d]=%s\n", i, fileSM[i].filename); } } //calculate memory for major variables struct_Segment * pS1 = & fileRaw[firstValidfile].m_Segment; struct_Segment * pS2 = & fileSM[firstOutputFile].m_Segment; float m1 = float(pS1 -> NX) * pS1 -> NY * NumberOfValidNDfiles / 1.0e6; //buffers input ND files float m2 = float(pS2 -> NX) * pS2 -> NY * SM_outputFiles / 1.0e6; //buffers for output files float memory_required = (m1 + m2) * 2; //both ND and SM are short integers printf(" Memory_required=%f MB for processing 'SM'\n", memory_required); printf(" input buffer size=[%d,%d], NumberOfValidNDfiles=%d , RAM=%f MB\n", pS1 -> NX, pS1 -> NY, NumberOfValidNDfiles, m1); printf(" output buffer size=[%d,%d], SM_outputFiles=%d RAM=%f MB\n", pS2 -> NX, pS2 -> NY, SM_outputFiles, m2); //allocate memory for (i = 0; i < FilterSize; i++) { if (inputFileStatus[i] == SUCCEED) fileRaw[i].processBuffer(ALLOCATE); if (SmFilesNeedOutput[i] == 1) { fileSM[i].processBuffer(ALLOCATE); fileSM[i].processBuffer(INITIALIZE); } } //read QA from first valid ND file, whole SDS will be read ; i = firstValidfile; fileRawQA.setFileToRead(filenames_ND[i], "ND", timeSeries_years[i], timeSeries_weeks[i], "QA;", 1); fileRawQA.processBuffer(ALLOCATE); fileRawQA.processBuffer(OPEN_AND_READ); uint8 * pNDQA = (uint8 * ) fileRawQA.Buffers[0]; Box box; int NP_valid = 0; //for SMN and SMT for (int id = 0; id < 2; id++) { char raw_sds_name[LENGTH_SDSNAME] = "NDVI"; if (id == 1) stringCopy(LENGTH_SDSNAME, raw_sds_name, "BT4"); char input_sds_name[LENGTH_SDSNAME] = ""; ADDSTRINGS3(LENGTH_SDSNAME, input_sds_name,raw_sds_name, ";"); printf(" processing '%s' ......\n", raw_sds_name); //open file for specific variable NDVI or BT4 for (i = 0; i < FilterSize; i++) { if (inputFileStatus[i] == SUCCEED) { if (setInputVHfile( & fileRaw[i], FILE_ND, timeSeries_years[i], timeSeries_weeks[i], input_sds_name, input_segments) == FAIL) continue; } } printf(" processing input files by %d segements\n", input_segments); for (int iSegment = 0; iSegment < input_segments; iSegment++) { int K0 = segment_Points * iSegment; fileRaw[firstValidfile].setSegment(iSegment); int NX = pSegment -> NX; int NY = pSegment -> ValidLines; long NP = pSegment -> ValidPoints; if( input_segments > 1 || id ==0){ NP_valid= getYrangeWithLandPixels(pNDQA +K0, NX, NY,m_pConfig -> sActiveTile, &box); } if(NP_valid == 0) { printf("skip reading ND for this segment, \n" ); continue; } // read a small section ( part of the segment), from line y1 to y2 in this segment int32 count [2]= {box.ny, box.nx}; int32 start [2]= {box.y1 + iSegment * pSegment -> NY ,box.x1}; int32 stride[2]= {1,1}; printf("start=[%d,%d]\ncount=[%d,%d]\n", start[0], start[1], count[0], count[1]); //setSegment for reading for (i = 0; i < FilterSize; i++) { if (inputFileStatus[i] == SUCCEED) { //fileRaw[i].setSegment(iSegment); //fileRaw[i].processBuffer(OPEN_AND_READ); fileRaw[i].open(fileRaw[i].filename, ACCESS_READ); fileRaw[i].IO_readSDS(raw_sds_name, start, stride, count, fileRaw[i].Buffers[0]); fileRaw[i].close(); } } char sType[LENGTH_FILENAME]; stringCopy(LENGTH_FILENAME, sType, fileRaw[firstValidfile].VARs[0].sds_name); if (showMessage(LEVEL2)) printf(" %s %s: segment=%d step 1: fill gaps by neighbor\n", task, sType, iSegment); //check each data field,fill gaps by neighbor , this is not need any more because the daily map has // already done. const int fill_gap_by_enighbor = 0; // set to 0 to disable below block if (fill_gap_by_enighbor == 1) { //fill gaps for (i = 0; i < FilterSize; i++) { if (inputFileStatus[i] == FAIL) continue; short * p0 = (short * ) fileRaw[i].Buffers[0]; short missing = fileRaw[i].VARs[0].FillValue; // fill gap along X direction, neighbor pixel only for (int iy = 0; iy < NY; iy++) { short * p = p0 + iy * NX; for (ix = (NX - 1); ix > 0; ix--) if (p[ix] == missing && p[ix - 1] != missing) p[ix] = p[ix - 1]; for (ix = 0; ix < (NX - 1); ix++) if (p[ix] == missing && p[ix + 1] != missing) p[ix] = p[ix + 1]; } if (showMessage(LEVEL2)) { print2Dsample(p0, NC_SHORT, fileRaw[i].filename, NX, NY, 8, 6); } } } if (showMessage(LEVEL2)){ printf(" %s %s: segment=%d step 2: Apply smoothing filter, NP=%ld, NX=%d,NY=%d\n", task, sType, iSegment, NP, pSegment -> NX, pSegment -> NY); } // smoothing time series ... int np= box.nx * box.ny; long NPmissing = 0; for (long k = 0; k < np; k++) { int y = k / box.nx; int L =K0 + (box.y1 + y) * NX + box.x1 + (k % box.nx); //skip water pixel if ((pNDQA[ L ] & 12) == 0) continue; // fill gap in the time series int totalBadValues = 0; float value; //check missing pixels for (i = 0; i < FilterSize; i++) { value = (inputFileStatus[i] == FAIL) ? MISSING_FLOAT : fileRaw[i].getValue(0, k); if (value <= MISSING_FLOAT) totalBadValues++; filter_pTS[i] = value; } if (totalBadValues >= (FilterSize - 2)) { NPmissing++; } else { //int minNPmissing=4; //if(totalBadValues > minNPmissing) printTimeSeries("inputTS",filter_pTS,FilterSize); FillGaps(FilterSize, filter_pTS, fmissing); //if(totalBadValues > minNPmissing) printTimeSeries("gapFilledTS",filter_pTS,FilterSize); MedianFilter(MedianFilterSize, FilterSize, filter_pTS, filter_pTemp); //if(totalBadValues > minNPmissing)printTimeSeries("MedianFilteredTS",filter_pTS,FilterSize); smoothTimeseriesByDigitalFilter(filter_weight, FilterSize, filter_pA, filter_pTS, filter_pTemp); //if(totalBadValues > minNPmissing){ // printTimeSeries("smoothedTS",filter_pTS,FilterSize); // exit(0); //} //save SMN and SMT of a pixel to buffers for (int i = 0; i < FilterSize; i++) { if (SmFilesNeedOutput[i] == 1) fileSM[i].setValue(id, L, filter_pTS[i]); } } } //k loop for pixels printf(" smoothing done for '%s' iSegment=%d (of %d) NP=%ld NPmissing=%ld\n", sType, iSegment + 1, input_segments, NP, NPmissing); } //Segment } //for SMN and SMT //set QA for (int i = 0; i < FilterSize; i++) { if (SmFilesNeedOutput[i] == 1) { // set QA flag uint8 * pSMQA = (uint8 * ) fileSM[i].Buffers[2]; for (long k = 0; k < fileSM[i].m_Segment.ValidPoints; k++) { //copy bit 1,2,3 (desert,land,coast) pSMQA[k] = (pNDQA[k] & 14); float SMN = fileSM[i].getValue(0, k); float SMT = fileSM[i].getValue(1, k); //set bit 1 if both SMN and SMT is invalid if (SMN <= MISSING_FLOAT && SMT <= MISSING_FLOAT) pSMQA[k] |= 1; //set bit 4 (too cold flag) if (SMT > MISSING_FLOAT && SMT <= 273) pSMQA[k] |= 16; } fileSM[i].saveVH_2NC( pSMQA ); fileSM[i].processBuffer(DEALLOCATE); } if (inputFileStatus[i] == SUCCEED) fileRaw[i].processBuffer(DEALLOCATE); } free(filter_pTS); fileRawQA.processBuffer(DEALLOCATE); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: composeWeeklyMap() * FUNCTION: make weekly GVI composite from daily maps * PROCEDURE: * - find the input daily maps * if the files with VHP style name exist, use it * otherwise, search the files in NDE style file name * - open the first daily map file available * - copy the file attributs, add file attributes for weekly composite file * - read SDS attibutes for each SDS in daily map * - create SDS with the same SDS attributes as daily map * - add two SDS to weekly composite map "cell_jday" and "ValidDaysForCH1" * - allocate buffers for daily and weekly composite. * - Major steps: * step 1: loop each day, read daily R1,R2 and update weekly * R1,R2, JDAY and ValidDaysForCH1 * step 2: write weekly R1,R2, JDAY and ValidDaysForCH1 to file * step 3: loop for other variables, * read variable, update it by JDAY, and write to file * save weekly composite for this variable to NetCDF file * - clear allocated memory * - add output file name to PSF file * INPUT: * year, week: the year and week number of the week being processed * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::composeWeeklyMap(char * task, int year, int week) { timeval time_begin; gettimeofday( & time_begin, nullptr); int jday1, jday2; getJdayRangeForPeriod(week, & jday1, & jday2); //found input Daily Maps by search directory int Ndays = m_pConfig -> DaysPerPeriod; //char *filenames_DailyMaps[Ndays]; //allocateMemory((VOIDP *)filenames_DailyMaps,Ndays,LENGTH_FILENAME); char ** filenames_DailyMaps = allocateCharArray(Ndays); char filename[LENGTH_FILENAME]; int total_inputFiles = 0; int firstValidJday = -1; //search from the file printf(" Search daily maps in directory '%s'\n", m_pConfig -> DIR_DailyMap); for (int jday = jday1; jday <= jday2; jday++) { //set default for file names as "" (not found) stringCopy(LENGTH_FILENAME, filenames_DailyMaps[jday - jday1], ""); int month, dayInMonth; getMonth(year, jday, & month, & dayInMonth); //search file in VHP file name convention (same as the AVHRR VHP files) VHP_getInputFileName( filename, FILE_DAILY, year, jday); if (FileExist(filename)) { stringCopy(LENGTH_FILENAME, filenames_DailyMaps[jday - jday1], filename); total_inputFiles++; if (firstValidJday < 0) firstValidJday = jday; } if (showMessage(LEVEL4)) printf(" jday=%d, file='%s' \n", jday, filenames_DailyMaps[jday - jday1] ); } printf(" %s(): total_inputFiles =%d, jday=[%d,%d], firstValidJday=%d \n", __FUNCTION__, total_inputFiles, jday1, jday2, firstValidJday); int NfilesRequired = int(m_pConfig -> DaysPerPeriod * 0.8); // for weekly, total_inputFiles < 5 if (total_inputFiles < NfilesRequired) { printf(" Making weekly composite map requires '%d' input daily maps\n", NfilesRequired); printf(" There are %d daily maps for this week.\n", total_inputFiles); printStatus(task, WARNING_WEEKLYGVI_NO_ENOUGH_DAILY_MAPS); freePointers((VOIDP * ) filenames_DailyMaps, Ndays); return (FAIL); } printf(" Making weekly composite map use the following input daily maps:\n"); for (int jday = jday1; jday <= jday2; jday++) { printf(" jday=%d, file='%s'\n", jday, filenames_DailyMaps[jday - jday1]); } int isAVHRR = isNoaaAvhrr(m_pConfig -> SatelliteID); //printf(" SatelliteID=%s, isAVHRR=%d\n", m_pConfig -> SatelliteID, isAVHRR); int jday = (week - 0.5) * m_pConfig -> DaysPerPeriod; if (isAVHRR == 1) loadCalibrationsForAVHRR(year, jday, & m_sAvhrrCalibration); // find output file name for weekly composite char filename_weeklyComposite[LENGTH_FILENAME] = ""; getFilename_VHPstyle(filename_weeklyComposite,FILE_GVI , year, week); // open outputfile // write file attributes // process ch1,ch2,solzen, senzen //extract options of the command line, and put the position parameters in //file info for daily and weekly maps CVhFile GVIdaily, GVIweekly; int ID_JDAY = -1; int ID_VALIDDAYSFORCH1 = -1; int ID_CM = -1; int num_datasets, num_global_attrs, sds_id, rank, num_attrs, dims[8]; structVAR * pVar; GVIdaily.setFileFormat(filenames_DailyMaps[firstValidJday - jday1]); GVIweekly.setFileFormat(filename_weeklyComposite); //check GVI files in this week GVIweekly.InputFileList_initialize(LENGTH_FILENAME * total_inputFiles); GVIweekly.open(filename_weeklyComposite, ACCESS_WRITE); if(showMessage(LEVEL1)) printf(" add filenames of Daily maps to weekly Composite map file\n"); for (int jday = jday1; jday <= jday2; jday++) { char * filename = filenames_DailyMaps[jday - jday1]; if (strlen(filename) == 0) continue; if(showMessage(LEVEL1)) printf(" read, %s\n", filename); GVIweekly.InputFileList_add(filename); }; //iFile //printf(" There are %d valid input files.\n%s\n\n", // GVIweekly.m_Meta.InputFiles, GVIweekly.m_Meta.pInputFilenames); printf(" read the first daily map file, copy the file attributs, read SDS attibutes \n"); int N_arraysToCopy = 0; for (int jday = jday1; jday <= jday2; jday++) { char * filename = filenames_DailyMaps[jday - jday1]; if (strlen(filename) == 0) continue; if(showMessage(LEVEL1)) printf(" read, %s\n", filename); if (jday != firstValidJday) continue; //get sds_names, size, copy file attributes to weekly map GVIdaily.open(filename, ACCESS_READ); GVIdaily.readSpatialRange( & GVIdaily.m_SDSlatlonRange); GVIdaily.IO_FileInfo( & num_datasets, & num_global_attrs); if(showMessage(LEVEL1)) printf(" num_datasets=%d, num_global_attrs=%d\n", num_datasets, num_global_attrs); ; //check if the last SDS is the special one "PLATE_CARREE" witch rank =0 for (int i = 0; i < num_datasets; i++) { int varid = i; pVar = & GVIdaily.VARs[varid]; GVIdaily.IO_SDgetinfo(varid, pVar -> sds_name, & rank, dims, & pVar -> sds_type, & num_attrs); if(showMessage(LEVEL1)) printf(" check rank: %-30s, rank= %d\n", pVar -> sds_name, rank); if (rank <= 1) continue; if (strcmp(pVar -> sds_name, "latitude") == 0 || strcmp(pVar -> sds_name, "longitude") == 0) continue; N_arraysToCopy++; } setWeeklyFileAttributes( & GVIweekly, m_pConfig -> SatelliteID, "Weekly Composite Map", year, week, GVIdaily.m_SDSlatlonRange); ID_JDAY = N_arraysToCopy - 1; ID_VALIDDAYSFORCH1 = N_arraysToCopy; ID_CM = N_arraysToCopy + 1; int dim_found = 0; if(showMessage(LEVEL1)) printf(" %s(), N_arraysToCopy= %d\n", __FUNCTION__, N_arraysToCopy); for (int i = 0; i < N_arraysToCopy; i++) { int varid = (i == (N_arraysToCopy - 1)) ? ID_CM : i; pVar = & GVIdaily.VARs[varid]; GVIdaily.IO_SDgetinfo(i, pVar -> sds_name, & rank, dims, & pVar -> sds_type, & num_attrs); if(showMessage(LEVEL1)) printf(" read, SDS=%d %s\n", i, pVar -> sds_name); if (rank == 2 && dim_found == 0) { GeoRange SDS_latlonRange; GVIdaily.readSpatialRange( & SDS_latlonRange); GVIweekly.setSDS(N_arraysToCopy + 2, dims[1], dims[0], SDS_latlonRange, 1, "weeklyMap"); GVIdaily.setSDS(N_arraysToCopy + 2, dims[1], dims[0], SDS_latlonRange, 1, "dailyMap"); dim_found = 1; } stringCopy(LENGTH_SDSNAME, GVIweekly.VARs[varid].sds_name, pVar -> sds_name); GVIweekly.VARs[varid].sds_type = pVar -> sds_type; sds_id = GVIdaily.IO_selectSDS(pVar -> sds_name); GVIdaily.readScaleParameters(sds_id, pVar); GVIdaily.readScaleParameters(sds_id, & GVIweekly.VARs[varid]); } GVIdaily.close(); }; //iFile if(showMessage(LEVEL1)) printf(" read the first daily map file, copy the file attributs(end)\n"); //printf(" There are %d valid input files.\n%s\n\n", // GVIweekly.m_Meta.InputFiles, GVIweekly.m_Meta.pInputFilenames); //add two SDS to weekly map pVar = & GVIweekly.VARs[ID_JDAY]; defineVariable(pVar, m_pConfig->output_cell_jday, "cell_jday", "Day", 1, m_pConfig -> DaysPerPeriod, 1, 1, 0, NC_BYTE, 0); pVar = & GVIweekly.VARs[ID_VALIDDAYSFORCH1]; defineVariable(pVar, 1, "ValidDaysForCH1", "Day", 1, m_pConfig -> DaysPerPeriod, 0, 1, 0, NC_BYTE, 0); // below: in case we read in daily Maps with azimuth angle but do not output in weekly composite map const int OutputAzimuth=0; if(OutputAzimuth ==0){ for (int iVariable = 0; iVariable < N_arraysToCopy -1; iVariable++) { structVAR * pVar = & GVIweekly.VARs[iVariable]; if(strstr(pVar -> sds_name,"azimuth") != nullptr || strstr(pVar -> sds_name,"reflectance_M3") != nullptr) { printf(" skip output this SDS: '%s'\n", pVar -> sds_name); pVar->NeedOutput= false; } } } //create SDS GVIweekly.createSDS(); GVIweekly.enddef(); //printf(" There are %d valid input files.\n%s\n\n", // GVIweekly.m_Meta.InputFiles, GVIweekly.m_Meta.pInputFilenames); PixelsCounter poll; initializePixelsCounter( & poll); int nx = GVIdaily.m_Segment.NX; int ny = GVIdaily.m_Segment.ValidLines; long np = long(nx) * ny; int ID_Band1 = 0; int ID_Band2 = 1; //we need 4 int16 buffers and 4 uint8 buffers = (np *12) bytes float memory_required = np / 1.0e6 * (4 * 2 + 4); int16 * pData = (int16 * ) ALLOCATE(np * 2 * 4); int16 * pDailyR1 = pData; int16 * pDailyR2 = pData + np; int16 * pWeeklyR1 = pData + np * 2; int16 * pWeeklyR2 = pData + np * 3; uint8 * pWeeklyJDAY = (uint8 * ) ALLOCATE(np * 4); uint8 * pWeeklyValidDaysForCH1 = pWeeklyJDAY + np; uint8 * pWeeklyCM = pWeeklyJDAY + np * 2; uint8 * pDailyCM = pWeeklyJDAY + np * 3; int16 FillValue = GVIdaily.VARs[ID_Band1].FillValue; printf(" Memory_required=%f MB for processing 'Weekly Composite', ", memory_required); printf(" array size= [%d %d], ", nx, ny); printf(" need 4 int16 buffers and 4 uint8 buffers\n"); //printf(" Initialize weekly R1 and R2 as FillValue= %d\n", FillValue); //initialize output for weekly R1 and R2 uint8 missing_CM = 1; for (int i = 0; i < np; i++) { *(pWeeklyR1 + i) = FillValue; *(pWeeklyR2 + i) = FillValue; *(pWeeklyJDAY + i) = 0; *(pWeeklyValidDaysForCH1 + i) = 0; *(pWeeklyCM + i) = missing_CM; } //step 1: loop each day, read daily R1,R2 and update weekly R1,R2, JDAY and ValidDaysForCH1 printf(" SkipWater =%d, ", m_pConfig -> SkipWater); printf(" Composite_Method_maxNDVI=%d\n", m_pConfig -> Composite_Method_maxNDVI); printf(" Step 1: update weekly R1,R2, JDAY, ValidDaysForCH1,CM \n"); for (int jday = jday1; jday <= jday2; jday++) { char * filename = filenames_DailyMaps[jday - jday1]; if (strlen(filename) == 0) continue; int scaled_day = jday - jday1 + 1; long pointsUpdated = 0; GVIdaily.open(filename, ACCESS_READ); //file does not exist; if (GVIdaily.m_fileID < 0) continue; printf(" read, GVIdaily, R1,R2,CM; "); GVIdaily.readWholeSDS(GVIdaily.VARs[ID_Band1].sds_name, pDailyR1); GVIdaily.readWholeSDS(GVIdaily.VARs[ID_Band2].sds_name, pDailyR2); GVIdaily.readWholeSDS(GVIdaily.VARs[ID_CM ].sds_name, pDailyCM); GVIdaily.close(); float fDailyR1, fDailyR2, daily_NDVI, R1, R2, weekly_NDVI; for (int i = 0; i < np; i++) { uint8 cm = * (pDailyCM +i); int isLandPixel = cm & 12; if(pWeeklyCM[i] == missing_CM) pWeeklyCM[i] = pDailyCM [i]; if(m_pConfig -> SkipWater ==1 && !isLandPixel ) continue; fDailyR1 = * (pDailyR1 + i); if (fDailyR1 <= 0) continue; fDailyR2 = * (pDailyR2 + i); if (fDailyR2 <= 0) continue; daily_NDVI = (isAVHRR == 1) ? m_sAvhrrCalibration.calculate_NDVI_For_AVHRR(fDailyR1, fDailyR2) : (fDailyR2 - fDailyR1) / (fDailyR2 + fDailyR1); if (daily_NDVI <= -1) continue; pWeeklyValidDaysForCH1[i]++; R1 = * (pWeeklyR1 + i); if (R1 > 0) { R2 = * (pWeeklyR2 + i); weekly_NDVI = (isAVHRR == 1) ? m_sAvhrrCalibration.calculate_NDVI_For_AVHRR(R1, R2) : (R2 - R1) / (R2 + R1); //printf("isAVHRR=%d daily_NDVI=%f weekly_NDVI=%f, %f %f %f %f, %d\n", isAVHRR, // daily_NDVI,weekly_NDVI,fDailyR1, fDailyR2, R1,R2,m_pConfig->Composite_Method_maxNDVI); if (m_pConfig -> Composite_Method_maxNDVI == 1 && (daily_NDVI <= weekly_NDVI)) continue; if (m_pConfig -> Composite_Method_maxNDVI == 0 && (fDailyR1 >= R1)) continue; } //this pixel is selected *(pWeeklyR1 + i) = * (pDailyR1 + i); *(pWeeklyR2 + i) = * (pDailyR2 + i); *(pWeeklyCM + i) = * (pDailyCM + i); *(pWeeklyJDAY + i) = scaled_day; pointsUpdated++; } poll.pointsUpdated += pointsUpdated; printf(" day= %3d pointsUpdated= %ld\n", jday, pointsUpdated); } //loop jday //step 2: write weekly R1,R2, JDAY and ValidDaysForCH1 to file int minCH1=0; GVIweekly.countMissingPixels( pWeeklyR1, pWeeklyCM, minCH1); GVIweekly.addStatistics(); printf(" Step 2: save weekly R1,R2, JDAY, ValidDaysForCH1 and CM to weekly file, \n"); GVIweekly.save1SDS( & GVIweekly.VARs[ID_Band1], (VOIDP) pWeeklyR1); GVIweekly.save1SDS( & GVIweekly.VARs[ID_Band2], (VOIDP) pWeeklyR2); if(m_pConfig->output_cell_jday ==1) GVIweekly.save1SDS( & GVIweekly.VARs[ID_JDAY ], (VOIDP) pWeeklyJDAY); GVIweekly.save1SDS( & GVIweekly.VARs[ID_VALIDDAYSFORCH1], (VOIDP) pWeeklyValidDaysForCH1); GVIweekly.save1SDS( & GVIweekly.VARs[ID_CM], (VOIDP) pWeeklyCM); //step 3: loop for other variables, // read variable, update it by JDAY, and write to file //for other band, read each other variable for section ROI and copy to daily map ; printf(" Step 3: process other bands: \n"); //process other variables for (int iVariable = 2; iVariable < N_arraysToCopy - 1; iVariable++) { structVAR * pVar = & GVIweekly.VARs[iVariable]; char * sds_name = pVar->sds_name; if(pVar->NeedOutput == false) continue; uint8 * pWeeklyGviUint8 = (uint8 * ) pWeeklyR1; uint8 * pDailyGviUint8 = (uint8 * ) pDailyR1; int16 * pWeeklyGviInt16 = pWeeklyR1; int16 * pDailyGviInt16 = pDailyR1; setArray(pVar -> sds_type, pWeeklyR1, np, pVar -> FillValue); int dailyFiles = 0; printf(" read, SDS= %-20s, jday= ", sds_name ); for (int jday = jday1; jday <= jday2; jday++) { char * filename = filenames_DailyMaps[jday - jday1]; if (strlen(filename) == 0) continue; int scaled_day = jday - jday1 + 1; GVIdaily.open(filename, ACCESS_READ); if (GVIdaily.m_fileID <= 0) continue; printf("%d, ", jday); int status = GVIdaily.readWholeSDS(sds_name, (VOIDP) pDailyR1); GVIdaily.close(); if (status == FAIL) continue; dailyFiles++; if (pVar -> sds_type == NC_CHAR || pVar -> sds_type == NC_BYTE) { for (int i = 0; i < np; i++) { if ( * (pWeeklyJDAY + i) != scaled_day) continue; if ( * (pDailyGviUint8 + i) == pVar -> FillValue) continue; *(pWeeklyGviUint8 + i) = * (pDailyGviUint8 + i); } } else { for (int i = 0; i < np; i++) { if ( * (pWeeklyJDAY + i) != scaled_day) continue; *(pWeeklyGviInt16 + i) = * (pDailyGviInt16 + i); } } } //day printf("\n"); if (dailyFiles > 0) GVIweekly.save1SDS(pVar, (VOIDP) pWeeklyR1); } GVIweekly.saveLatLonArray(); GVIweekly.close(); //free buffer free(pData); free(pWeeklyJDAY); freePointers((VOIDP * ) filenames_DailyMaps, Ndays); printNameOfSavedFile(__FUNCTION__,GVIweekly.filename, & time_begin); return (SUCCEED); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: getInputWeeksForSmooth() * FUNCTION: get begin and end week of input time series data for smoothing * INPUT: * year,week: the reference week * OUTPUT: * year1, week1: the begin week of the time series * year2, week2: the end week of the time series * LIMITATION/RULE: * The time series should include the reference week. * We will check the computer time and get the latest week number to update * the end of availble data record. * usually, the time series for smoothing is the N weeks around reference * week(N=filetsize), and the reference week is at the center of time series. * If the reference week is near the end of availble data record, we will * adjusted the time series, so that * (1) make the reference week close to center of time serious as possible * (2) have enough valid weeks for smoothing. * For example: for AVHRR, the available data record covers from 1981 week 35 * to present. if the reference week is 1981 week 36, and filetsize= 15, then, * the time sereis for smoothing will be set as * from 1981 week 35 to 1981 week 49. * For the last 8 weeks, the time series will be the 15 weeks which ends at * the latest week. -----------------------------------------------------------------------------*/ void CVHsuite::getInputWeeksForSmooth(int year, int week, int & year1, int & week1, int & year2, int & week2) { /* if(m_pConfig->smooth_usingPeriod[0]>0){ year1= m_pConfig->smooth_usingPeriod[0]; week1= m_pConfig->smooth_usingPeriod[1]; year2= year1; week2= week2+ m_pConfig->FilterSize; adjustPeriod(year2,week2); return; } */ //======================================================== // New addition by ASSISTT for testing // ======================================================= char* env_ending_week = std::getenv("ENDING_WEEK"); char* env_ending_year = std::getenv("ENDING_YEAR"); if (env_ending_week != nullptr) { printf("**** ADDITIONAL ASSISTT code for testing******"); int ending_week_val = atoi(env_ending_week); // Default to the input year if ENDIING_YEAR explicitly provided int ending_year_val = year; if (env_ending_year != nullptr){ ending_year_val = atoi(env_ending_year); } // set the end of the 15 week window using the ENV variables year2 = ending_year_val; week2 = ending_week_val; // check if year2 and week2 falls within 1 to 52 weeks in a year adjustPeriod(year2, week2); // check if year2 and week2 falls within 1 to 52 weeks in a year adjustPeriod(year2, week2); // calculate the start of the window by strictly counting backward year1 = year2; week1 = week2 - (m_pConfig->FilterSize - 1); // handle cross year if the start week drop below 1 adjustPeriod(year1, week1); printf(" [ENV OVERRIDE] ENDING_YEAR=%d, ENDING_WEEK=%d detected.\n", ending_year_val, ending_week_val); printf(" Overriding window: %04d week %02d to %04d week %02d\n", year1, week1, year2, week2); return; // exit and bypass the original VHP centering logic } // ==== end of ASSISTT additional code for testing ========================= // Below is special procedure to avoid there are too many missing files // (weeks) in time series to be smoothed /* 1985 week 1 to 1985 week 8 1994 week 37 to 1995 week 3 2004 week 11 to 2004 week 24 */ if ((strcmp(m_pConfig -> DataSource, "VHP") == 0 || strcmp(m_pConfig -> DataSource, "GVIX") == 0) && m_pConfig -> FilterSize == 15) { year1 = year; year2 = year; //for known begin or end of a period with data ; if (year == 1984 && week >= (52 - 7) && week <= 52) { week2 = 49; week1 = week2 - 14; return; } if (year == 1985 && week <= (9 + 7)) { week1 = 9; week2 = week1 + 14; return; } if (year == 1994 && week >= (36 - 7)) { week2 = 36; week1 = week2 - 14; return; } if (year == 1995 && week <= (4 + 7)) { week1 = 4; week2 = week1 + 14; return; } if (year == 2004 && week >= (10 - 7) && week <= 10) { week2 = 10; week1 = week2 - 14; return; } if (year == 2004 && (week >= 11 && week < (25 + 7))) { week1 = 25; week2 = week1 + 14; return; } } //find the start year and week from available data list //(input through configure file: vh.config) int startYearWeek = m_pGVI_dataPeriods -> yearPeriod1[0]; int startYear = startYearWeek / 100; int startWeek = startYearWeek % 100; int FilterSize = m_pConfig -> FilterSize; //year0,week0 is temporary use, int year0 = startYear; int week0 = startWeek + FilterSize / 2; adjustPeriod(year0, week0); //determine input periods: year1,week1,year2,week2 //at the begin of time series if ((year * 100 + week) <= (year0 * 100 + week0)) { //for begin of the AVHRR time series year1 = year; week1 = startWeek; year2 = year; week2 = startWeek + FilterSize - 1; adjustPeriod(year2, week2); } else { year1 = year; week1 = week - FilterSize / 2; year2 = year; week2 = week + FilterSize / 2; adjustPeriod(year1, week1); adjustPeriod(year2, week2); //check if the week is near the end of the AVHRR time series int latest_Year, latest_Week; getPreviousPeriod(m_pConfig -> DaysPerPeriod, latest_Year, latest_Week); if ((year2 * 100 + week2) > (latest_Year * 100 + latest_Week)) { year2 = latest_Year; week2 = latest_Week; year1 = year2; week1 = week2 - (FilterSize - 1); adjustPeriod(year1, week1); } } } /*----------------------------------------------------------------------------- * FUNCTION_NAME: loadCalibrationsForAVHRR() * FUNCTION: load calibration for a particular week from pre-launch and * post-launch calibration files * PROCEDURE: * - find the satellite ID for the given week * - load pre-launch calibration for this satellite * - load most recent post-launch calibration for this satellite and * interpolate to the give jday * INPUT: * year,week: the particular week * OUTPUT: * pCali: the pointer to structure CAvhrrCalibration which holds the * calibration parameters. * RETURN: * satellite number -----------------------------------------------------------------------------*/ int CVHsuite::loadCalibrationsForAVHRR(int year, int jday, CAvhrrCalibration * pCali) { char * satID; int satNumber; satID = getSatelliteID(FILE_DAILY, year, jday, & satNumber); if (satNumber <= 0) { printf("WARNING: Failed to load calibration due to invalid satellite %d jday=%d: satNumber=%d\n", year, jday, satNumber); } else { //read pre-launch calibration file pCali -> readPrelaunchCalibration(m_pConfig -> File_preLaunchCalibration, satID); //read post-launch calibration file pCali -> getPostLaunchCalibrationForAVHRR(m_pConfig -> File_postLaunchCalibration, satNumber, year, jday); //show calibration parameters //if(showMessage(LEVEL2)) pCali -> showCalibrationString(); } return satNumber; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: getSatelliteID() * FUNCTION: Get the satellite ID for a given ween based on GVI data available. * INPUT: * year: day of year * jday: day-of-year * OUTPUT: * satNumber: the satellite number of data used for this week * RETURN: * satellite ID of the data used for this week * if there is no data for this week, then, satNumber=1 and return nullptr * NOTE: * GVI data available is a table described which satellite we should use * for each period. * This is described in configure file in the section titled as * "[Periods of GVI data used for VH]" * example: * [Periods of GVI data used for VH] * #satID satNumber yearWeek1 yearWeek2 * NC 07 198135 198449 * NF 09 198509 198845 * NH 11 198846 199436 * NJ 14 199504 200052 * NL 16 200101 200401 * NL 16 200405 200410 * NL 16 200425 200428 * NL 16 200430 200523 * NN 18 200524 201252 * NP 19 201301 399999 * The above information was loaded when loading the configure file by function * readDataPeriodsFromConfigureFile(), and save in structure * Periods* m_pGVI_dataPeriods; * The following shows ID and number for each satellite * satellite, ID, number * NOAA7 NC, 7 * NOAA9 NF, 9 * NOAA11 NH, 11 * NOAA14 NJ, 14 * NOAA16 NL, 16 * NOAA17 NM, 17 * NOAA18 NN, 18 * NOAA19 NP, 19 * Soumi NPP npp * New feature: 20240909: if the above configure section "[Periods of GVI data used for VH]" is not defined in configure file, VHsuite will use the same satellite specified on command line -----------------------------------------------------------------------------*/ char * CVHsuite::getSatelliteID(int fileType, int year, int period, int * satNumber) { // when configure section "[Periods of GVI data used for VH]" is not defined in congifure file, //use the same satellite specified on command line if(m_pConfig -> SatellitePeriodsDefined ==0){ *satNumber = 999; return m_pConfig-> SatelliteID; } int myJday = (fileType == FILE_DAILY) ? period : (period - 0.5) * m_pConfig -> DaysPerPeriod; int week = (myJday - 1) / 7 + 1; int yrwk = year * 100 + week; Periods * p = m_pGVI_dataPeriods; for (int i = 0; i < p -> N_validPeriods; i++) { if ((yrwk >= p -> yearPeriod1[i] && yrwk <= p -> yearPeriod2[i])) { * satNumber = p -> satNumber[i]; //printf(" getSatelliteID(): satelliteID found, i=%d satID= %s yrwk=%d, %d %d\n", // i, p->satID[i], yrwk, p->yearPeriod1[i],p->yearPeriod2[i]); return p -> satID[i]; } } ; //not found * satNumber = -1; printf(" getSatelliteID(): satelliteID not found for %d week %d\n", year, week); return nullptr; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: isNoaaAvhrr() * FUNCTION: check if it is NOAA satellie by satellite ID * INPUT: * SatelliteID: a string representing satellite ID, such as "NP", "npp" * RETURN: * 1, if it is NOAA satellite * 0, if not * NOTE: * ID for known NOAA satellites are * The following shows ID and number for each satellite -----------------------------------------------------------------------------*/ int CVHsuite::isNoaaAvhrr(const char * SatelliteID) { char AVHRR4GVI[LENGTH_FILENAME] = "NC,NF,NH,NJ,NK,NL,NM,NN,NP,M1,M2,M3"; if (strstr(AVHRR4GVI, SatelliteID) != nullptr) return 1; return 0; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: EDF_HistoToEDF() * FUNCTION: calculate EDF values from histogram of data * PROCEDURE: * - calculate the realaive histogrm by dividing the histogram value by total * number of pixels * - calulate the EDF ( accumulative histogram for bin 0 to i) * INPUT: * Histo: array of histogram of data. * NBINS: number of bins * OUTPUT: * EDF: EDF values * EDF[i] was defined as Sum(Histogram(0:i])/ total(Histogram) -----------------------------------------------------------------------------*/ int CVHsuite::EDF_HistoToEDF(float * Histo, float * EDF, int32 NBINS) { //calculat total valid pixels int i; float np = 0; for (i = 0; i < NBINS; i++) np += Histo[i]; if (np <= 0) { //printf("ERROR in EDF_HistoToEDF(): np =%d <=0\n", np); return (-1); } // calculate normalized histogram for (i = 0; i < NBINS; i++) Histo[i] /= np; //calculate EDF EDF[0] = Histo[0]; for (i = 1; i < NBINS; i++) EDF[i] = EDF[i - 1] + Histo[i]; return (0); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: EDF_Adjustment() * FUNCTION: calculate adjustment amount based on EDF difference betweeen * particualr data set and bench mark data set * NOTE: * For each EDF bin of NDVI, calculate adjustment amount of scaled NDVI. * For a given NDVI value (NDVIi), find corresponding CDF value, * then find the NDVI (NDVIj) with the same EDF in N-years benchmark * data field, the difference between NDVIj and NDVIi is the adjustment * amount of NDVI. * adjustment = (NDVIj - NDVIi) * INPUT: * NBINS: the number of bins * EDF: the EDF values of particular dataset * EDF_benchMark: the EDF values of bench mark data set * dNDVIperBin: interval of NDVI per statistic bins * OUTPUT: * Adjustment: the amount of NDVI adjusted -----------------------------------------------------------------------------*/ void CVHsuite::EDF_Adjustment(int32 NBINS, float * EDF, float * EDF_benchMark, float dNDVIperBin, float * Adjustment) { float edf; Adjustment[0] = 0; for (int i = 1; i < NBINS; i++) { edf = EDF[i]; for (int j = 0; j < (NBINS - 1); j++) { if (edf >= EDF_benchMark[j] && edf < EDF_benchMark[j + 1]) { Adjustment[i] = (j - i) * dNDVIperBin; break; } } } } /*----------------------------------------------------------------------------- * FUNCTION_NAME: EDF_CorrectLineByLine() * FUNCTION:EDF correction on NDVI by latitude line * EDF correction is also called CDF * EDF: Empirical Distribution Function * CDF: Cumulative Distribution Function * This function calculate EDF from the NDVI field (read from file) * Then adjust NDVI by EDF-benchmark * EDF-benchmark was calculated based on NDVI data of selected years * for each week. Detail instruction was written by Marco Vargas as below: *[Marco Vargas said:] *As per our meeting with Felix this afternoon, here a list of some of the things *that we agreed to implement: *1. EDF correction of NDVI weekly images 4 KM resolution (before smoothing) * for the following years: 1982, 1983, 1984, 1985, 1986, 1987, 1988, * 1991, 1992, 1993, 1994, 1999 and 2000. *2. To build the standard or benchmark EDFs use the NDVI from the following * years: 1989, 1990, 1995, 1996, 1997. *3. The NDVI threshold for applying the correction is 0.01 * The BT threshold for applying the correction is 4.0 K * For example, for given NDVI (value= NDVI1), find the EDF value at NDVI1 * from the EDF profile of this week, then find the NDVI_standard value * which have same EDF value in the EDF profile for the benchmark data set. * ajustment = NDVI_standard -NDVI1 //no correction * (if ajustment < threshold) ajustment=0 ; * NDVI_new= NDVI1 + ajustment * Similar correction will applied to BT. *4. A single benchmark EDF will be used for all latitudes between 75N - 60N. * The statistical sample for this benchmark EDF will be composed by all * the valid NDVI pixels from those latitudes from the five images that * conform the particular benchmark week. *5. A single benchmark EDF will be used for all latitudes between 30S to 55S. * The statistical sample for this benchmark EDF will be composed by all * the valid NDVI pixels from those latitudes from the five images that * conform the particular benchmark week. *6. Later it will be decided how to fill the 1994 data gap -----------------------------------------------------------------------------*/ int CVHsuite::EDF_CorrectLineByLine(CVhFile * pFileND, const char * sInputType, int32 year, int32 week) { if (showMessage(LEVEL1)) { printf(" EDF_CorrectLineByLine() for '%s' %4d week %2d ......\n", sInputType, year, week); } enum FILE_SDSID_RAW { ID_NDVI, ID_BT4 }; char * fnameStandardEDF = m_pConfig -> File_EDF_NDVI; int id_input = ID_NDVI; int EDFbins; //marco: for NDVI, if adjustment >0.01, then apply adjustment. //marco: for BT, if adjustment >4K, then apply adjustment. float AdjustmentThreshold = 0.01; if (strcmp(sInputType, "BT4") == 0) { fnameStandardEDF = m_pConfig -> File_EDF_BT4; id_input = ID_BT4; AdjustmentThreshold = 4; } if (showMessage(LEVEL1)) { printf(" EDF_CorrectLineByLine: "\ "fnameStandardEDF='%s', id_input=%d AdjustmentThreshold=%f\n", fnameStandardEDF, id_input, AdjustmentThreshold); } int i; int NX = pFileND -> m_Segment.NX; int NY = pFileND -> m_Segment.ValidLines; float32 MinMaxInterval[3]; //read NDVIcounts_BenchMark from HDF file with resolution of 4km int sd_id = SDstart(fnameStandardEDF, DFACC_READ); if (sd_id < 0) { ERROR_OPEN_FILE(fnameStandardEDF); } int32 start[2] = { 0, 0 }; int32 stride1D[1] = { 1 }, edge1D[1] = { 3 }; long sds_id; char sds_name[LENGTH_SDSNAME]; snprintf(sds_name, LENGTH_SDSNAME, "%s_MinMaxInterval", sInputType); sds_id = SDselect(sd_id, SDnametoindex(sd_id, sds_name)); SDreaddata(sds_id, start, stride1D, edge1D, MinMaxInterval); SDendaccess(sds_id); float zmin = MinMaxInterval[0]; float zmax = MinMaxInterval[1]; float dNDVIperBin = MinMaxInterval[2]; EDFbins = int((zmax - zmin) / dNDVIperBin); //allocate buffer float * pBuffer = (float * ) ALLOCATE(EDFbins * 5 * 4); //histogram for benchmark float * Histo_Benchmark = pBuffer; //EDF profile from benchmark float * EDF_benchMark = pBuffer + EDFbins; //(6 years SMN) //histogam of SMN of the week being processed. float * Histo = pBuffer + EDFbins * 2; //EDF profile from the week being processed. float * EDF = pBuffer + EDFbins * 3; //adjusted amount of Scaled NDVI for each NDVI level float * Adjustment = pBuffer + EDFbins * 4; int16 * NDVIcounts_BenchMark = (int16 * ) ALLOCATE(NY * EDFbins * 2); //read NVIcounts_ByLine_Gxx.hdf snprintf(sds_name, LENGTH_SDSNAME, "%s_Counts_Week%02d", sInputType, week); sds_id = SDselect(sd_id, SDnametoindex(sd_id, sds_name)); int sample = 3616 / NY; int32 stride2D[2] = { sample, 1 }, edge2D[2] = { NY, EDFbins }; SDreaddata(sds_id, start, stride2D, edge2D, NDVIcounts_BenchMark); SDendaccess(sds_id); SDend(sd_id); if (showMessage(LEVEL1)) { printf(" read counts_BenchMark: sample=%d, "\ "%s[%d,%d]\n", sample, sds_name, NY, EDFbins); printf(" min=%f, max=%f, interval=%f\n", zmin, zmax, dNDVIperBin); printf(" read,%s OK\n", fnameStandardEDF); } //Bench mark has been read int32 ix, iLevel, j1, j2; // this is the latitude range which data is used to calcualte EDF // consistent to Java VH // pos1 and pos2 are the start and end pixels in the data array which data // will be used to calculate EDF profile float lat_North = pFileND -> m_SDSlatlonRange.maxLat; float lat_South = pFileND -> m_SDSlatlonRange.minLat; float dLatitudePerLine = (lat_North - lat_South) / NY; int y1 = int((lat_North - 60.) / dLatitudePerLine); int y2 = int((lat_North - (-30.)) / dLatitudePerLine); printf(" 60N: y1=%d, 30S: y2=%d\n", y1, y2); //get the data range and scaling parameters for scaling SMN data to integer int32 NP_adjusted = 0; for (int y = y1; y <= y2; y++) { j1 = y; j2 = y; if (y == y1) { j1 = 0; j2 = y1; } if (y == y2) { j1 = y2; j2 = NY - 1; } int n = EDFbins * 5; for (i = 0; i < n; i++) pBuffer[i] = 0; for (i = 0; i < EDFbins; i++) { for (int j = j1; j <= j2; j++) Histo_Benchmark[i] += NDVIcounts_BenchMark[i + j * EDFbins]; } EDF_HistoToEDF(Histo_Benchmark, EDF_benchMark, EDFbins); int pos1 = j1 * NX; int pos2 = (j2 + 1) * NX - 1; //statistics on histogram for (ix = pos1; ix <= pos2; ix++) { float value = pFileND -> getValue(id_input, ix); if (value <= MISSING_FLOAT) continue; iLevel = int((value - zmin) / dNDVIperBin); if (iLevel < 0 || iLevel >= EDFbins) continue; Histo[iLevel]++; } EDF_HistoToEDF(Histo, EDF, EDFbins); EDF_Adjustment(EDFbins, EDF, EDF_benchMark, dNDVIperBin, Adjustment); if (showMessage(LEVEL4) && (y == y1 || y == y2)) { printf("type=%s range[%f,%f],interval=%f Bins=%d\n", sInputType, zmin, zmax, dNDVIperBin, EDFbins); printf("y=%d\n", y); printf("EDF=["); for (i = 0; i < (EDFbins - 1); i++) printf("%f,", EDF[i]); printf("%f]\n", EDF[EDFbins - 1]); printf("EDF_benchMark=["); for (i = 0; i < (EDFbins - 1); i++) printf("%f,", EDF_benchMark[i]); printf("%f]\n", EDF_benchMark[EDFbins - 1]); printf("Adjustment=["); for (i = 0; i < (EDFbins - 1); i++) printf("%f,", Adjustment[i]); printf("%f]\n", Adjustment[EDFbins - 1]); } for (ix = pos1; ix <= pos2; ix++) { float value = pFileND -> getValue(id_input, ix); if (value <= MISSING_FLOAT) continue; //printf("xxx %s ix=%d value=%f zmin=%f dNDVIperBin=%f\n,", // sInputType, ix,value,zmin, dNDVIperBin); iLevel = int((value - zmin) / dNDVIperBin); if (iLevel < 0 || iLevel >= EDFbins) continue; if (Adjustment[iLevel] < AdjustmentThreshold) continue; //printf("xxx %s ix=%d value=%f iLevel=%d Adjustment=%f\n,", // sInputType, ix,value,iLevel, Adjustment[iLevel]); value += Adjustment[iLevel]; pFileND -> setValue(id_input, ix, value); NP_adjusted++; } //ix } //y free(pBuffer); free(NDVIcounts_BenchMark); printf(" EDF_CorrectLineByLine() for '%s' finished. NP_adjusted=%d\n", sInputType, NP_adjusted); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: createAvhrrClimate() * FUNCTION: create Climatology file for a given week; * climatology includes: * [MXN, MNN, meanSMN, stddevSMN] and * [MXT, MNT, meanSMT, stddevSMT] * RAM requested for 4km product: 72MB *(8 +1 +3 *2) = 1.08GB * PROCEDURE: * - determine the data used (periods for each satellite, * defined in configure file by section titled as * "[Periods of AVHRR data used for GVI climatology]"); * - read the first file and get the file attributes (SDS size); * - create output file and 8 SDS; * - allocate buffers (8 outputs, 1 input); * - open each file * - loop for each input segement * read input * aggregate the 8 output varables * - close input file * - write output * - close output file * - free buffers * INPUT: * filename_PCF: configure file * week: week -----------------------------------------------------------------------------*/ int CVHsuite::createClimatologyFiles(char * task, int week) { printf("Calculate: 'CLIMAT' week=%2d\n", week); //determine the data used (periods for each satellite) Periods MaxMinPeriods; readDataPeriodsFromConfigureFile(m_pConfig -> File_Configure, "[Periods of data used for GVI climatology]", & MaxMinPeriods); int input_segments = 1; CVhFile fileSM, fileClimat; /* open the first smoothing file and load SDS parameters, it was required this file must exist!*/ char *satID = MaxMinPeriods.satID[0]; int yearweek1 = MaxMinPeriods.yearPeriod1[0]; int year1 = yearweek1 / 100; int week1 = yearweek1 % 100; printf("open the first smoothing file and load SDS parameters, "); printf("it was required this file must exist, year=%d week=%d\n", year1, week1); stringCopy(LENGTH_20, m_pConfig -> SatelliteID, satID); if (setInputVHfile( & fileSM, FILE_SM, year1, week1, "SMN;SMT;QA;", input_segments) == FAIL) return FAIL; printf("%s(), read array size from first SM file= %s\n", __FUNCTION__,fileSM.filename); char sds_names_CLIMAT[LENGTH_FILENAME] = "NDVI_MAX;NDVI_MIN;NDVI_MEAN;NDVI_STDDEV;BT_MAX;BT_MIN;BT_MEAN;BT_STDDEV;QA;Samples;"; printf(" nx=%d ny=%d\n", fileSM.m_SDScolumns, fileSM.m_SDSrows); printf(" fileSM.m_SDSlatlonRange.minLat=%f\n", fileSM.m_SDSlatlonRange.minLat); long NP = fileSM.m_Segment.ValidPoints; /* memory_required: Input: SM: 2 x int16 +1 x uint8 Output: CLM: 9 x int16 +1 x uint8 internal: 3 xfloat */ float memory_required = NP / 1.0E6 * (5 + 19 + 12); printf("Total Memory_required=%f MB for processing 'Making climatology'\n", memory_required); //get the input file names int SM_files = 0; const int maxYears = 100; int SM_years[maxYears]; char SM_satIDs[maxYears][LENGTH_20]; char SM_filenames[maxYears][LENGTH_FILENAME]; memset(&SM_years[0], 0, sizeof(int)); memset(SM_satIDs[0], 0, LENGTH_20 * maxYears); memset(SM_filenames[0], 0, LENGTH_FILENAME * maxYears); fileClimat.InputFileList_initialize(LENGTH_FILENAME * maxYears); snprintf(fileClimat.m_Meta.pInputFilenames, LENGTH_FILENAME , "Input files:\n" ); for (int i= 0; i < MaxMinPeriods.N_validPeriods; i++){ char *satID = MaxMinPeriods.satID[i]; int yearweek1 = MaxMinPeriods.yearPeriod1[i]; int yearweek2 = MaxMinPeriods.yearPeriod2[i]; int year1= yearweek1 / 100; int year2= yearweek2 / 100; for (int year = year1; year <= year2; year++) { int32 yearweek = year * 100 + week; if (yearweek < yearweek1 || yearweek > yearweek2) continue; stringCopy(LENGTH_20, m_pConfig -> SatelliteID, satID); char *filename = SM_filenames[SM_files]; VHP_getInputFileName( filename, FILE_SM, year, week); if (FileExistWarning(filename) != 1) continue; printf("%s(), SM file found: file= %s\n", __FUNCTION__,filename); fileClimat.InputFileList_add(filename); stringCopy(LENGTH_20, SM_satIDs[ SM_files ] , satID); SM_years[ SM_files ] = year; SM_files++; } } //create output file and add file attributes stringCopy(LENGTH_20, m_pConfig -> SatelliteID, "npp+j01+j02"); if (setOutputVHfile( & fileClimat, FILE_CLIMAT, 0, week, sds_names_CLIMAT, & fileSM) == FAIL) return FAIL; enum FILE_SDSID_CLIMAT { ID_MXN, ID_MNN, ID_MEAN_ND, ID_STDDEV_ND, ID_MXT, ID_MNT, ID_MEAN_BT, ID_STDDEV_BT, ID_CLIMATQA, ID_SAMPLE }; printf(" fileSM.processBuffer(ALLOCATE);\n"); fileSM.processBuffer(ALLOCATE); printf(" fileClimat.processBuffer(ALLOCATE);\n"); fileClimat.processBuffer(ALLOCATE); float * pClimatDataBuffer = (float * ) ALLOCATE(NP * 3 * 4); float * pmean = pClimatDataBuffer; float * pstddev = pClimatDataBuffer + NP; float * pNpixels = pClimatDataBuffer + NP * 2; /* read data and calculate climatology, then save to climatology file */ printf("Create Climatology file for week %d ......\n", week); fileClimat.processBuffer(INITIALIZE); uint8 * pSmoothQA = (uint8 * ) fileSM.Buffers[2]; // for SMN and SMT for (int id = 0; id <= 1; id++) { char sds_name[LENGTH_20]="SMN"; if(id == 1) stringCopy(LENGTH_20, sds_name, "SMT"); short FillValue = fileSM.VARs[id].FillValue; int id_max = (id == 0) ? ID_MXN : ID_MXT; int id_min = (id == 0) ? ID_MNN : ID_MNT; short * pSM = (short * ) fileSM.Buffers[id]; short * pmax = (short * ) fileClimat.Buffers[id_max]; short * pmin = (short * ) fileClimat.Buffers[id_min]; for (long j = 0; j < NP; j++) { pmean[j] = 0; pstddev[j] = 0; pNpixels[j] = 0; pmax[j] = FillValue; pmin[j] = FillValue; } /* loop for all posible years */ for (int i= 0; i < SM_files; i++){ /* open file for SMN or SMT */ stringCopy(LENGTH_20, m_pConfig -> SatelliteID, SM_satIDs[i]); if (setInputVHfile( & fileSM, FILE_SM, SM_years[i], week, "SMN;SMT;QA;", input_segments) == FAIL) continue; printf("%s(), reading, SM file= %s, sds_name= %s\n", __FUNCTION__,fileSM.filename, sds_name); if(fileSM.readWholeSDS(fileSM.filename, sds_name, pSM) == FAIL) continue; if(i==0 && id==0) { fileSM.readWholeSDS(fileSM.filename, "QA", pSmoothQA); printf("%s(), reading, SM file= %s, sds_name= %s\n", __FUNCTION__,fileSM.filename, "QA"); } for (long j = 0; j < NP; j++) { int16 smn = pSM[j]; if (smn != FillValue) { if (pmax[j] == FillValue) { pmax[j] = smn; pmin[j] = smn; } else { if (pmax[j] < smn) pmax[j] = smn; if (pmin[j] > smn) pmin[j] = smn; } float value = fileSM.getValue(id, j); pmean[j] += value; pstddev[j] += (value * value); pNpixels[j]++; } } }// iSat //compute mean and stddev: int id_mean = (id == 0) ? ID_MEAN_ND : ID_MEAN_BT; int id_stddev = (id == 0) ? ID_STDDEV_ND : ID_STDDEV_BT; for (long j = 0; j < NP; j++) { if (pmax[j] == FillValue) pmin[j] = FillValue; if (pNpixels[j] > 0) { float m = pmean[j] / pNpixels[j]; fileClimat.setValue(id_mean, j, m); float sigma = pstddev[j] / pNpixels[j] - m * m; sigma = (sigma > 0) ? sqrt(sigma) : 0; fileClimat.setValue(id_stddev, j, sigma); } } } //id , for SMN and SMT int16 * pClimatNP = (int16 * ) fileClimat.Buffers[ID_SAMPLE]; uint8 * pClimatQA = (uint8 * ) fileClimat.Buffers[ID_CLIMATQA]; for (long k = 0; k < NP; k++) { //copy bit 1,2,3 (desert,land,coast) pClimatQA[k] = (pSmoothQA[k] & 14); pClimatNP[k] = pNpixels[k]; int VaraiblesMissing = 0; for (int id = 0; id < 8; id++) { if (fileClimat.getValue(id, k) <= MISSING_FLOAT) VaraiblesMissing++; } //set bit 1 if all variables are invalid if (VaraiblesMissing == 8) pClimatQA[k] |= 1; //set bit 4 (too cold flag) float maxT = fileClimat.getValue(ID_MXT, k); if (maxT > MISSING_FLOAT && maxT <= 273) pClimatQA[k] |= 16; } fileClimat.saveVH_2NC( pClimatQA ); fileClimat.processBuffer(DEALLOCATE); fileSM.processBuffer(DEALLOCATE); free(pClimatDataBuffer); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: getFilename_VHPstyle() * FUNCTION: get VHP style file name for a VH file * short file name: File name is as the same style as AVHRR VHP, example as * VGVI.G500m.C01.npp.P2014001_r01c01.nc for FILE_DAILY * VGVI.G500m.C07.npp.P2014001.S001.E007_r01c01.nc for FILE_GVI * VGVI.G500m.C01.npp.P2014001.ND_r01c01.nc for FILE_ND * VGVI.G500m.C01.npp.P2014001.SM_r01c01.nc for FILE_SM * VGVI.G500m.C01.npp.P2014001.VH_r01c01.nc for FILE_VH * where, rRRcCC is row and column number for a tile (for 1km or * finer resolution) * VGVI.G04.C01.npp.P2014001.nc for FILE_DAILY * VGVI.G04.C07.npp.P2014001.S001.E007.nc for FILE_GVI * VGVI.G04.C01.npp.P2014001.ND.nc for FILE_ND * VGVI.G04.C01.npp.P2014001.SM.nc for FILE_SM * VGVI.G04.C01.npp.P2014001.VH.nc for FILE_VH * INPUT: * fileType: output file type * year, week: output year and week * OUTPUT: * filename_output: output file name -----------------------------------------------------------------------------*/ void CVHsuite::getFilename_VHPstyle(char * filename_output, int fileType, int year, int week) { char fileTag[LENGTH_SDSNAME], sFileType[LENGTH_SDSNAME]; char * folder = getFolder(fileType, fileTag, sFileType); char NewFolder[LENGTH_FILENAME] = ""; stringCopy(LENGTH_FILENAME, NewFolder, folder); char * satID; int satNumber; satID = getSatelliteID(fileType, year, week, & satNumber); if (strlen(m_pConfig -> sActiveTile) > 0) { //if(m_pConfig -> bRunWithPCF ==0){ char subFolder[LENGTH_FILENAME]=""; snprintf(subFolder, LENGTH_FILENAME, "%s/%4d/%4d%03d", satID, year, year, week); SETSTRING(LENGTH_FILENAME, NewFolder, folder, "/", subFolder); //} } char fname1[LENGTH_FILENAME]; int jday1, jday2; char sFileTypeTag[LENGTH_SDSNAME]; stringCopy(LENGTH_SDSNAME, sFileTypeTag, ""); int ndays = (fileType == FILE_DAILY) ? 1 : m_pConfig -> DaysPerPeriod; snprintf(fname1, LENGTH_FILENAME, "%s.%s.C%02d.%s.P%04d%03d", m_pConfig -> Prefix, m_pConfig -> sResolution_VHP, ndays, satID, year, week); switch (fileType) { case FILE_CLIMAT: snprintf(fname1, LENGTH_FILENAME, "%s.%s.C%02d.CLIMAT.P%02d", m_pConfig -> PrefixCLIMAT, m_pConfig -> sResolution_VHP, ndays, week); if (strlen(m_pConfig -> sActiveTile) > 0) { char file1[LENGTH_2000]; //WG2025 on 20250720 YAN CHANGED the climatology file TO *.nc : // /data/data043/DATA_LUO/JPSS/500m/climate_v2025_long snprintf(file1, LENGTH_2000, "P%03d/%s%s", week, fname1, m_pConfig -> sActiveTile); SETSTRING(LENGTH_FILENAME, filename_output, folder, "/", file1); } else { SETSTRING(LENGTH_FILENAME, filename_output, folder, "/", fname1); } if(strlen(m_pConfig -> ClimatFileExtension) ==0) stringCopy(LENGTH_20,m_pConfig -> ClimatFileExtension,".nc"); ADDSTRINGS(LENGTH_FILENAME, filename_output, m_pConfig -> ClimatFileExtension); return; case FILE_GVI: getJdayRangeForPeriod(week, & jday1, & jday2); snprintf(sFileTypeTag, LENGTH_SDSNAME, ".S%03d.E%03d", jday1, jday2); break; case FILE_ND: stringCopy(LENGTH_SDSNAME, sFileTypeTag, ".ND"); break; case FILE_SM: stringCopy(LENGTH_SDSNAME, sFileTypeTag, ".SM"); break; case FILE_VH: stringCopy(LENGTH_SDSNAME, sFileTypeTag, ".VH"); break; case FILE_DAILY: stringCopy(LENGTH_SDSNAME, sFileTypeTag, ""); break; default: break; } //snprintf(filename_output,LENGTH_FILENAME,"%s/%s%s%s.nc",NewFolder,fname1,sFileTypeTag,m_pConfig->sActiveTile ); SETSTRING(LENGTH_FILENAME, filename_output, NewFolder, "/", fname1); ADDSTRINGS3(LENGTH_FILENAME, filename_output, sFileTypeTag, m_pConfig -> sActiveTile); ADDSTRINGS(LENGTH_FILENAME, filename_output, ".nc"); //printf("%s(): filename_output=%s\n",__FUNCTION__,filename_output ); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: getFolder() * FUNCTION: get directory name for a VH file type * INPUT: * fileType: output file type * RETURN: * directory name -----------------------------------------------------------------------------*/ char * CVHsuite::getFolder(int fileType, char * fileTag, char * sFileType) { char FileTypeStrings[7][30] = { "DAILY", "GVI", "ND", "SM", "VH", "GVF", "CLIMAT" }; char FileTags[7][30] = { "DAILY_MAP", "WEEKLY_COMPOSITE", "ND_FILE", "SM_FILE", "VH_FILE", "GVF_FILE", "CLIMATE_FILE" }; stringCopy(LENGTH_SDSNAME, fileTag, FileTags[fileType]); stringCopy(LENGTH_SDSNAME, sFileType, FileTypeStrings[fileType]); if (fileType == FILE_DAILY) return m_pConfig -> DIR_DailyMap; if (fileType == FILE_GVI) return m_pConfig -> DIR_WeeklyComposite; if (fileType == FILE_ND) return m_pConfig -> DIR_ND; if (fileType == FILE_SM) return m_pConfig -> DIR_SM; if (fileType == FILE_VH) return m_pConfig -> DIR_VH; if (fileType == FILE_CLIMAT) return m_pConfig -> DIR_CLIMAT; return nullptr; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: setInputVHfile() * FUNCTION: set parameters for Output VH file * INPUT: * pInputVHfile : object of CVhFile class for input file * fileType: output file type * year, week: output year and week * sds_names: output sds names separated by ";" * input_segments: input segments (if >1, the whole SDS will be processed by segements; * each segment covers NY/input_segments lines * OUTPUT: * paramters of pInputVHfile will be set. -----------------------------------------------------------------------------*/ int CVHsuite::setInputVHfile(CVhFile * pInputVHfile, int fileType, int year, int week, const char * sds_names, int input_segments) { char fileTag[LENGTH_SDSNAME], sFileType[LENGTH_SDSNAME]; getFolder(fileType, fileTag, sFileType); // VHP short file name , if does not exist, search last NDE file with specific type) VHP_getInputFileName( pInputVHfile -> filename, fileType, year, week); if (FileExistWarning(pInputVHfile -> filename) != 1) { printf("file '%s' does not exist\n", pInputVHfile -> filename); return (FAIL); } if (pInputVHfile -> setFileToRead(pInputVHfile -> filename, sFileType, year, week, sds_names, input_segments) == FAIL) { return FAIL; } if(showMessage(LEVEL1)) printf(" %s, %s=%s\n", __FUNCTION__, fileTag, pInputVHfile -> filename); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: setOutputVHfile() * FUNCTION: set parameters for Output VH file * INPUT: * pOutputVHfile : object of CVhFile class for output file * fileType: output file type * year, week: output year and week * sds_names: output sds names separated by ";" * inputVHfile: object of CVhFile class for input file * OUTPUT: * paramters of pOutputVHfile will be set. -----------------------------------------------------------------------------*/ int CVHsuite::setOutputVHfile(CVhFile * pOutputVHfile, int fileType, int year, int week, const char * sds_names, CVhFile * inputVHfile) { char fileTag[LENGTH_SDSNAME], sFileType[LENGTH_SDSNAME]; char * folder = getFolder(fileType, fileTag, sFileType); char * satID; int satNumber; //printf(" setOutputVHfile(): year=%d,week=%d folder=%s\n", year, week, folder); satID = getSatelliteID(fileType, year, week, & satNumber); if(m_pConfig -> UseNDEfilename ==1) { getFileName_NDE(pOutputVHfile -> filename, fileType, year, week); }else { getFilename_VHPstyle(pOutputVHfile -> filename, fileType, year, week); } int variables = pOutputVHfile -> setOutputVaraibles(sds_names); char mapName[LENGTH_FILENAME]; if (year == 0) snprintf(mapName, LENGTH_FILENAME, "%s_Week_%02d", sFileType, week); else snprintf(mapName, LENGTH_FILENAME, "%s_%04d_Week_%02d", sFileType, year, week); if(showMessage(LEVEL1)) printf(" %s(): mapName= %s folder= %s\n", __FUNCTION__, mapName, folder); pOutputVHfile -> setFileFormat(pOutputVHfile -> filename); int TotalSegments = 1; pOutputVHfile -> setSDS(variables, inputVHfile -> m_SDScolumns, inputVHfile -> m_SDSrows, inputVHfile -> m_SDSlatlonRange, TotalSegments, mapName); /* start HDF file */ if (pOutputVHfile -> open(pOutputVHfile -> filename, ACCESS_WRITE) == FAIL) { ERROR_EXIT(ERROR_VH_CREATE_FILE, pOutputVHfile -> filename); } setWeeklyFileAttributes(pOutputVHfile, satID, mapName, year, week, inputVHfile -> m_SDSlatlonRange); pOutputVHfile -> createSDS(); pOutputVHfile -> enddef(); pOutputVHfile -> close(); pOutputVHfile -> showVariableDescription(); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: setWeeklyFileAttributes() * FUNCTION: set File attributes to weekly Output VH file * INPUT: * pOutputVHfile : object of CVhFile class for output file * satID: satellite ID * productname: product name; * year, week: output year and week * sds_names: output sds names separated by ";" * inputVHfile: object of CVhFile class for input file * OUTPUT: * File attributes of pOutputVHfile will be set. -----------------------------------------------------------------------------*/ int CVHsuite::setWeeklyFileAttributes(CVhFile * pOutputVHfile, const char * satID, const char * productname, int year, int week, GeoRange latlonRange) { if(showMessage(LEVEL1)) printf(" %s()\n", __FUNCTION__); struct_FileMetadata * pMeta = & pOutputVHfile -> m_Meta; char * s = pMeta -> Ancillary_Files; stringCopy(LENGTH_SDSNAME, s, ""); if (satID != nullptr && isNoaaAvhrr(satID) == 1) { addFileName(s, "File_preLaunchCalibration", m_pConfig -> File_preLaunchCalibration); addFileName(s, "File_postLaunchCalibration", m_pConfig -> File_postLaunchCalibration); addFileName(s, "File_EDF_NDVI", m_pConfig -> File_EDF_NDVI); addFileName(s, "File_EDF_BT4", m_pConfig -> File_EDF_BT4); } addFileName(s, "File_Configure", m_pConfig -> File_Configure); addFileName(s, "File_PCF", m_pConfig -> File_PCF); addFileName(s, "File_Metadata_Regions", m_pConfig -> File_Metadata_Regions); addFileName(s, "File_IGBPLandTypes", m_pConfig -> File_IGBPLandTypes); addFileName(s, "File_AST_SurfaceTypes", m_pConfig -> File_AST_SurfaceTypes); addFileName(s, "File_AST_WaterFraction", m_pConfig -> File_AST_WaterFraction); addFileName(s, "File_Static_Metadata", m_pConfig -> File_Static_Metadata); addFileName(s, "DIR_ND", m_pConfig -> DIR_ND); addFileName(s, "DIR_SM", m_pConfig -> DIR_SM); addFileName(s, "DIR_VH", m_pConfig -> DIR_VH); // addFileName(s, "DIR_PartialDailyMap", m_pConfig -> DIR_PartialDailyMap); addFileName(s, "DIR_DailyMap", m_pConfig -> DIR_DailyMap); addFileName(s, "DIR_WeeklyComposite", m_pConfig -> DIR_WeeklyComposite); addFileName(s, "DIR_CLIMAT", m_pConfig -> DIR_CLIMAT); // addFileName(s, "DIR_Archive", m_pConfig -> DIR_Archive); // addFileName(s, "DIR_Meta", m_pConfig -> DIR_Meta); // addFileName(s, "DIR_Temp", m_pConfig -> DIR_Temp); stringCopy(LENGTH_FILENAME, pMeta -> ProductName, productname); pMeta -> year = year; pMeta -> period = week; pMeta -> DaysPerPeriod = m_pConfig -> DaysPerPeriod; pMeta -> latlonRange = latlonRange; if (satID != nullptr) stringCopy(LENGTH_SDSNAME, pMeta -> Satellite, satID); else stringCopy(LENGTH_SDSNAME, pMeta -> Satellite, "NA"); stringCopy(LENGTH_SDSNAME, pMeta -> projection, "Plate Carree"); int jday1, jday2; getJdayRangeForPeriod(week, & jday1, & jday2); setTimeString(pMeta -> Time_Begin, year, jday1, "0000000"); setTimeString(pMeta -> Time_End, year, jday2, "2359599"); pOutputVHfile -> saveMetadata(m_pConfig -> sActiveTile, m_pConfig -> File_Static_Metadata, m_pConfig -> File_Configure, m_pConfig -> File_PCF); return SUCCEED; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: create_sampled_VH_file() * FUNCTION: make sample file for GVI, ND,SM or VH file * NOTE: for version 4.3.2 or earlier, NETCDF function nc_get_vars() does not * work well (very slow) when stride is not 1 * thus, we read SDS as a whole * INPUT: * sample : sample rate * tilesX,tilesY: numbers of tiles in X and y direction * inputfilesPattern: Pattern of input files, * filenameSampled: name of output sampled file * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::create_sampled_VH_file(const char * task, int sample, int tilesX, int tilesY, char * inputfilesPattern, char * filenameSampled, char * dir_meta, int aggregate_method) { timeval time_begin; gettimeofday( & time_begin, nullptr); printf(" %s():\n input pattern=%s\n sampled file =%s\n dir_meta =%s\n", __FUNCTION__, inputfilesPattern, filenameSampled, dir_meta); int num_datasets, num_global_attrs, sds_id, rank, num_attrs, dims[8]; //copy file attibute from the fist file name //copy the SDS attributes //loop for tiles and get the sample data, set to sampled data array //Example of VHP file name //file1 ="VGVI.G04.C07.npp.P2014001.VH" int maxTiles = tilesX * tilesY; //seach input files int Nfiles; char pattern[LENGTH_FILENAME]; snprintf(pattern, LENGTH_FILENAME, "%s%s", inputfilesPattern, "*"); char ** input_filenames = getFilesWithPattern(m_pConfig -> DIR_Temp, pattern, & Nfiles); if (Nfiles == 0) { printf(" WARNING: %s() : failed to find file with pattern='%s',\n ", __FUNCTION__, pattern); printStatus(task, ERROR_SAMPLE_FILE_NO_INPUT); return FAIL; } if (Nfiles != maxTiles) { char strError[LENGTH_2000]; snprintf(strError, LENGTH_2000, " Warning: %s() :number of files found by pattern='%s' "\ "(Nfiles=%d) is not as the as expected (maxTiles=%d),\n ", __FUNCTION__, pattern, Nfiles, maxTiles); printStatus(task, strError); ;//return FAIL; } CVhFile VhOriginFile, VhSampleFile; structVAR * pVar; //from the first file, //copy the file attributes except latlonrange and filename //put file attributes for latlonrange and filename //find the number of SDS //create buffer for output //create SDS for output char * filename = input_filenames[0]; printf(" get file information from the first file: %s\n", filename); VhOriginFile.setFileFormat(filename); if (VhOriginFile.open(filename, ACCESS_READ) == FAIL) { printf(" WARNING: %s() : failed because can not open file %s,\n ", __FUNCTION__, filename); printStatus(task, ERROR_SAMPLE_FILE_NO_INPUT); VhOriginFile.close(); return FAIL; }; VhOriginFile.IO_FileInfo( & num_datasets, & num_global_attrs); printf(" num_datasets=%d, num_global_attrs=%d\n", num_datasets, num_global_attrs); if (num_datasets <= 0) { printf(" WARNING: %s() failed because num_datasets = %d ", __FUNCTION__, num_datasets); printStatus(task, ERROR_SAMPLE_FILE_NO_INPUT); VhOriginFile.close(); return FAIL; } num_datasets -= 2; //do not process the last 2 SDS, which is "latitude" and "longitude" int varid = 0; const int maxArrays = 200; int varids[maxArrays]; char * sds_names[maxArrays]; for (int i = 0; i < maxArrays; i++) sds_names[i] = (char * ) ALLOCATE(LENGTH_SDSNAME); VhOriginFile.IO_getVariableIDs(num_datasets, varids, sds_names); for (int i = 0; i < maxArrays; i++) free(sds_names[i]); int Nregions_for_metadata = 0; Region * regions_for_metadata = readRegionsList(m_pConfig -> File_Metadata_Regions, Nregions_for_metadata); VhSampleFile.setFileFormat(filenameSampled); VhSampleFile.open(filenameSampled, ACCESS_WRITE); VhSampleFile.InputFileList_initialize(LENGTH_FILENAME * maxTiles); VhSampleFile.xml_initialize(num_datasets, Nregions_for_metadata, regions_for_metadata); printf(" copy (%d) file attributes from the first file: %s\n", num_global_attrs, filename); int isNDEstyle = 0; //copying file attributes for (int i = 0; i < num_global_attrs; i++) { char attr_name[LENGTH_FILENAME]; int attr_type, count; VhOriginFile.IO_getAttributeNameByIndex(NC_GLOBAL, i, attr_name); VhOriginFile.IO_getAttributeInfoByName(NC_GLOBAL, attr_name, & attr_type, & count); char attr_type_name[LENGTH_FILENAME]; VhOriginFile.getStringOfNCtype(attr_type, attr_type_name); if (strcmp(attr_name, "FILENAME") == 0) continue; if (strcmp(attr_name, ATTRNAME_STARTLAT) == 0) continue; if (strcmp(attr_name, ATTRNAME_ENDLAT) == 0) continue; if (strcmp(attr_name, ATTRNAME_STARTLON) == 0) continue; if (strcmp(attr_name, ATTRNAME_ENDLON) == 0) continue; if (strstr(attr_name, "geospatial_") != nullptr) { isNDEstyle = 1; continue; } int32 lsize = count * 4; char * attr_buf = ALLOCATE(lsize); VhOriginFile.IO_readAttributeByName(NC_GLOBAL, attr_name, attr_buf); //below is from the first tile, not cover all tiles, if(strcmp(attr_name,"INPUT_FILES") ==0) continue; if(strcmp(attr_name,"INPUT_FILENAMES") ==0) continue; VhSampleFile.IO_put_attribute(NC_GLOBAL, attr_name, attr_type, count, attr_buf); if (showMessage(LEVEL2)) printf(" put, attribute, %d attr_name=%s value=%s\n", i, attr_name, attr_buf); free(attr_buf); } // i, for copying file attributes //add file attributes: FILENAME, and LatlongRange printf(" add file attribute: FILENAME=%s\n", filenameSampled); VhSampleFile.IO_put_attribute(NC_GLOBAL, "FILENAME", NC_CHAR, strlen(filenameSampled), filenameSampled); VhSampleFile.IO_put_attribute(NC_GLOBAL, "INPUT_FILES", NC_INT, 1, & Nfiles); GeoRange r; VhOriginFile.readSpatialRange( & r); r.minLat = r.maxLat - tilesY * (r.maxLat - r.minLat); r.maxLon = r.minLon + tilesX * (r.maxLon - r.minLon); printf(" add file attribute: latlonrange=[%f,%f,%f,%f]\n", r.minLat, r.minLon, r.maxLat, r.maxLon); printf(" get dimension from first SDS: \n"); //get size of first SDS int valid_datasets = 0; int NX = -1; int NY = -1; for (varid = 0; varid < num_datasets; varid++) { pVar = & VhOriginFile.VARs[varid]; sds_id = varid; if (VhOriginFile.m_FileFormat == HDF4) sds_id= SDselect (VhOriginFile.m_fileID , varid); VhOriginFile.IO_SDgetinfo(sds_id, pVar -> sds_name, & rank, dims, & pVar -> sds_type, & num_attrs); //printf("1 sds_name= %s sds_id=%d\n",pVar -> sds_name, sds_id); //read file again, copy the SDS information to the structure for output file pVar = & VhSampleFile.VARs[varid]; VhOriginFile.IO_SDgetinfo(sds_id, pVar -> sds_name, & rank, dims, & pVar -> sds_type, & num_attrs); if (rank <= 1) continue; if (strcmp(pVar -> sds_name, "latitude") == 0 || strcmp(pVar -> sds_name, "longitude") == 0) continue; sds_id = VhOriginFile.IO_selectSDS(pVar -> sds_name); //printf("2 sds_name= %s sds_id=%d\n",pVar -> sds_name, sds_id); VhOriginFile.readScaleParameters(sds_id, pVar); printf(" will read, SDS=%d %s [%d,%d]\n", varid, pVar -> sds_name, dims[1], dims[0]); valid_datasets++; if (valid_datasets == 1) { NX = dims[1]; NY = dims[0]; } } VhOriginFile.close(); num_datasets = valid_datasets; int dims_sampled[2]; dims_sampled[0] = (NY * tilesX) / sample; dims_sampled[1] = (NX * tilesY) / sample; printf(" dims_sampled=[%d,%d] sample=%d tilesX=%d tilesY=%d\n", dims_sampled[1], dims_sampled[0], sample, tilesX, tilesY); VhOriginFile.setSDS(num_datasets, NX, NY, r, 1, "OriginMap"); printf(" Take sample for each SDS: \n"); int Nx_sampled = dims_sampled[1]; int Ny_sampled = dims_sampled[0]; long np_sampled = Ny_sampled * Nx_sampled; int16 * pDataOrigin = (int16 * ) ALLOCATE(NX * NY * sizeof(int16)); int16 * pDataSample = (int16 * ) ALLOCATE(np_sampled * sizeof(int16)); float * pSum = (float * ) ALLOCATE(np_sampled * sizeof(float)); int16 * pNPoints = (int16 * ) ALLOCATE(np_sampled * sizeof(int16)); printf(" memory allocated: np_sampled=%ld \n", np_sampled); //add input file names to output file as file attributes struct_FileMetadata * pMeta = & VhSampleFile.m_Meta; snprintf(pMeta -> pInputFilenames, pMeta -> __maxLengthForInputFilenames, "pattern='%s'\n\nInput files:\n", pattern); for (int iTile = 0; iTile < Nfiles; iTile++) { VhSampleFile.InputFileList_add(getFileBaseName(input_filenames[iTile])); if(iTile <2 || iTile == (Nfiles-1) ) printf(" iTile=%2d file=%s \n", iTile, input_filenames[iTile]); if(iTile ==1) printf(" ...\n"); } VhSampleFile.setSDS(num_datasets, Nx_sampled, Ny_sampled, r, 1, "SampledMap"); VhSampleFile.addAttributesForLatLonRange( & r, isNDEstyle); VhSampleFile.createSDS(); VhSampleFile.enddef(); /* ; //add quality_information ; printf(" add quality_information: \n"); ; ; PixelsCounter poll; ; initializePixelsCounter(&poll); ; ; long pointsOverLand = 0; ; long missingOverLand = 0; ; //int32 sds_id; ; for(int i=0; i sds_name, "latitude") == 0 || strcmp(pVar -> sds_name, "longitude") == 0) continue; setArray(pVar -> sds_type, (VOIDP) pDataSample, np_sampled, pVar -> FillValue); setArray(pSum, np_sampled, 0.0); setArray(NC_SHORT, (VOIDP) pNPoints, np_sampled, 0); int isMax = (strstr(pVar -> sds_name, "_MAX") != nullptr) ? 1 : 0; int isMin = (strstr(pVar -> sds_name, "_MIN") != nullptr) ? 1 : 0; if(showMessage(LEVEL3)) printf("%s() process, varid= %d, sds_name= %s\n", __FUNCTION__, varid, pVar -> sds_name); for (int iTile = 0; iTile < maxTiles; iTile++) { short int row = iTile / tilesX; short int col = iTile % tilesX; if(maxTiles >1){ //input is multiple tiles filename = nullptr; char tileTag[LENGTH_30]; snprintf(tileTag, LENGTH_30,"_r%02dc%02d", row+1, col+1); // WG 20241210, allow missing files // WG 20250401, if there are multiple tiles, pick the last one for (int K=0; K sds_name); }else{ //input is 1 big array (example: sampling 4km ND ==>16km filename= input_filenames[0]; } VhOriginFile.open(filename, ACCESS_READ); //printf( " read, VhOriginFile.m_fileID=%d file=%s\n",VhOriginFile.m_fileID, filename); sds_id = varid; if (VhOriginFile.m_FileFormat == HDF4) sds_id= SDselect (VhOriginFile.m_fileID , varid); VhOriginFile.IO_SDgetinfo(sds_id, pVar -> sds_name, & rank, dims, & pVar -> sds_type, & num_attrs); if (rank <= 1) continue; timeval time_begin1; gettimeofday( & time_begin1, nullptr); int nx = dims[1]; int ny = dims[0]; //x0,x1,x2, y0,y1,y2 is the position in big array for all tiles //i1,i2,j1,j2 is the position in sample array int x0 = nx * (col); int y0 = ny * (row); int i1 = int(float(x0 - sample / 2) / sample + 0.5); int x1 = i1 * sample + sample / 2; int i2 = int(float(x0 + nx - 1 - sample / 2) / sample); int x2 = i2 * sample + sample / 2; int j1 = int(float(y0 - sample / 2) / sample + 0.5); int y1 = j1 * sample + sample / 2; int j2 = int(float(y0 + ny - 1 - sample / 2) / sample); int y2 = j2 * sample + sample / 2; int Nx1 = (i2 - i1 + 1); int Ny1 = (j2 - j1 + 1); //xx1,yy1 is the srat position in the array of this tile int xx1 = x1 - x0; int yy1 = y1 - y0; // this is the size of data piece in small array we are going to set int32 sampled_dims [2] = {Ny1,Nx1}; int32 sampled_start [2] = {yy1, xx1}; int32 sampled_stride[2] = {sample,sample}; if(showMessage(LEVEL3)) { printf("%s(), start:[%d, %d], dim:[%5d,%5d], sample=%d ==> set to big array: x:[%6d:%6d], [%6d:%6d] \n", __FUNCTION__, xx1,yy1, Nx1, Ny1,sample, i1,i2, j1,j2); } if (pVar -> sds_type == NC_SHORT || pVar -> sds_type == NC_USHORT) { if (aggregate_method == 0) { /* ; VhOriginFile.IO_readSDS(pVar -> sds_name, start, stride, dims, (VOIDP) pDataOrigin); ; ; //for sampling only: copy the points ; for (int j = 0; j < Ny1; j++) { ; int16 * p1 = pDataOrigin + (j * sample + y1 - y0) * nx + x1 - x0; ; int16 * p2 = pDataSample + (j + j1) * Nx_sampled + i1; ; for (int i = 0; i < Nx1; i++) * (p2 + i) = * (p1 + i * sample); ; } */ VhOriginFile.IO_readSDS(pVar->sds_name, sampled_start, sampled_stride, sampled_dims,(VOIDP) pDataOrigin); //for sampling only: copy the points for (int j = 0; j < Ny1; j++) { int16 * p1 = pDataOrigin + j * Nx1 ; int16 * p2 = pDataSample + (j + j1) * Nx_sampled + i1; for (int i = 0; i < Nx1; i++) * (p2 + i) = * (p1 + i ); } } else { //for averaging //x,y is the position in tile //i,j is position in sampled array VhOriginFile.IO_readSDS(pVar -> sds_name, start, stride, dims, (VOIDP) pDataOrigin); for (int y = 0; y < ny; y++) { int j = (y0 + y) / sample; long pos = j * Nx_sampled; int16 * pInt16 = (int16 * ) pDataOrigin + y * nx; uint16 * pUint16 = (uint16 * ) pDataOrigin + y * nx; for (int x = 0; x < nx; x++) { int i = (x0 + x) / sample; int value = (pVar -> sds_type == NC_SHORT) ? ( * (pInt16 + x)) : ( * (pUint16 + x)); if (value != pVar -> FillValue) { int k = pos + i; if (isMax == 1) { if ( * (pDataSample + k) < value) * (pDataSample + k) = value; } else if (isMin == 1) { if ( * (pDataSample + k) > value) * (pDataSample + k) = value; } else { *(pSum + k) += value; *(pNPoints + k) += 1; } } } //x } //y } // aggregate_method } if (pVar -> sds_type == NC_CHAR || pVar -> sds_type == NC_BYTE) { /* ; VhOriginFile.IO_readSDS(pVar -> sds_name, start, stride, dims, (VOIDP) pDataOrigin); ; for (int j = 0; j < Ny1; j++) { ; uint8 * p1 = (uint8 * ) pDataOrigin + (j * sample + y1 - y0) * nx + x1 - x0; ; uint8 * p2 = (uint8 * ) pDataSample + (j + j1) * Nx_sampled + i1; ; for (int i = 0; i < Nx1; i++) * (p2 + i) = * (p1 + i * sample); ; } */ VhOriginFile.IO_readSDS(pVar->sds_name, sampled_start, sampled_stride, sampled_dims,(VOIDP) pDataOrigin); for (int j = 0; j < Ny1; j++) { uint8 * p1 = (uint8 * ) pDataOrigin + (j ) * Nx1 ; uint8 * p2 = (uint8 * ) pDataSample + (j + j1) * Nx_sampled + i1; for (int i = 0; i < Nx1; i++) * (p2 + i) = * (p1 + i ); } } VhOriginFile.close(); if(showMessage(LEVEL3)) { printf(" iTile=%2d, row=%4d, col=%5d, j=[%5d,%5d], i=[%5d,%5d],np=[%d,%d], "\ "x=[%5d,%5d],y=[%5d,%5d],isMax=%d, isMin=%d, ", iTile, row, col, j1, j2, i1, i2, Nx1, Ny1, x1, x2, y1, y2, isMax, isMin); printTimeUsed( & time_begin1); } } if (aggregate_method == 1 && isMax == 0 && isMin == 0) { ; // calcaulte averaging value by sum and number of valid points if (pVar -> sds_type == NC_SHORT) { int16 * p = (int16 * ) pDataSample; for (long i = 0; i < np_sampled; i++) { int n = * (pNPoints + i); if (n > 0) * (p + i) = int(( * (pSum + i)) / n + 0.5); } } else if (pVar -> sds_type == NC_USHORT) { uint16 * p = (uint16 * ) pDataSample; for (long i = 0; i < np_sampled; i++) { int n = * (pNPoints + i); if (n > 0) * (p + i) = int(( * (pSum + i)) / n + 0.5); } } } // aggregate_method == 1 if (rank == 2) { printf(" write SDS %-25s start=[%d,%d],stride=[%d,%d],count=[%d,%d] \n", pVar -> sds_name, start[1], start[0], stride[1], stride[0], dims_sampled[1], dims_sampled[0]); VhSampleFile.IO_writeSDS(pVar -> sds_name, start, stride, dims_sampled, (VOIDP) pDataSample); } VhSampleFile.xml_PerformStatisticsOf1VariableForRegions(varid, pDataSample); } //varid, loop for each SDS VhOriginFile.close(); VhSampleFile.saveLatLonArray(); VhSampleFile.close(); VhSampleFile.xml_saveMetadataToXmlFile(dir_meta); //release buffer free_charArray(input_filenames, Nfiles); free1D(__FUNCTION__, (VOIDP &)pDataOrigin, "pDataOrigin"); free1D(__FUNCTION__, (VOIDP &)pDataSample, "pDataSample"); free1D(__FUNCTION__, (VOIDP &)pSum, "pSum"); free1D(__FUNCTION__, (VOIDP &)pNPoints, "pNPoints"); free1D(__FUNCTION__, (VOIDP &)regions_for_metadata, "regions_for_metadata"); printNameOfSavedFile(task, VhSampleFile.filename, & time_begin); return (SUCCEED); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: adjustPeriod() * FUNCTION: adjust year and period number, so that period falls in [1,PeriodsPerYear] * INPUT: * year,period: year and period number * OUTPUT: * year, period: the adjusted year and period number -----------------------------------------------------------------------------*/ void CVHsuite::adjustPeriod(int & year, int & period) { int N = m_pConfig -> PeriodsPerYear; if (period <= 0) { period += N; year--; } if (period > N) { period -= N; year++; } } /*----------------------------------------------------------------------------- * FUNCTION_NAME: getJdayRangeForPeriod() * FUNCTION: get range of Jday for a Period * INPUT: * year,period: year and period number * OUTPUT: * jday1,jday2: the range of jday -----------------------------------------------------------------------------*/ void CVHsuite::getJdayRangeForPeriod(int period, int * jday1, int * jday2) { * jday1 = (period - 1) * m_pConfig -> DaysPerPeriod + 1; * jday2 = period * m_pConfig -> DaysPerPeriod; } /*----------------------------------------------------------------------------- * FUNCTION_NAME: split_1km_ND_file() * FUNCTION: split 1 global ND file to 4x4 tiles * NOTE: The purpose is to provide substitute ND files for OSPO * * INPUT: * * filename: global ND file * outputFolder: output folder * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::split_1km_ND_file(char * filename,char *outputFolder ) { timeval time_begin; gettimeofday( & time_begin, nullptr); printf(" %s():\n input file=%s \n outputFolder=%s\n", __FUNCTION__,filename, outputFolder); //copy file attibute from the fist file name //copy the SDS attributes //loop for tiles and get the sample data, set to sampled data array //Example of VHP file name //filename ="data/j01/1km/VGVI.G04.C07.j01.P2014001.ND.nc" const int tilesX = 4; const int tilesY = 4; char filepath1[LENGTH_FILENAME]; char *filebase= getFileBaseName(filename); SETSTRING(LENGTH_FILENAME,filepath1, outputFolder,"/",filebase); replace(filepath1,".nc", ""); char satID[LENGTH_20]; int year, week; char *p= strstr( filebase, ".C07.") +5; sscanf(p, "%3s.P%04d%03d", satID, &year,&week); printf("satId=%s, year=%d week=%d filepath1=%s\n", satID, year,week, filepath1); int NX= 10000; int NY= 3616; int NP=NX *NY; int Global1km_NX= NX *4; int Global1km_NY= NY *4; int Global1km_NP= Global1km_NX * Global1km_NY; int16 *pGlobal1km_NDVI = (int16 *) ALLOCATE( Global1km_NP *2 * sizeof(int16)); int16 *pGlobal1km_BT = pGlobal1km_NDVI + Global1km_NP; uint8 *pGlobal1km_QA =(uint8 *) ALLOCATE( Global1km_NP ); int16 *pNDVI = (int16 *) ALLOCATE( NP*2 * sizeof(int16)); int16 *pBT = pNDVI +NP; uint8 *pQA = (uint8 *) ALLOCATE( NP ); CVhFile VhOriginFile; CVhFile ncFileOutput; CVhFile *pOutFile =& ncFileOutput; int num_datasets, num_global_attrs; if(VhOriginFile.open(filename, ACCESS_READ) == FAIL) return FAIL; VhOriginFile.IO_FileInfo(&num_datasets, &num_global_attrs); printf("num_datasets=%d, num_global_attrs=%d\n", num_datasets, num_global_attrs); GeoRange GVI_latlonRange; VhOriginFile.readSpatialRange( & GVI_latlonRange); VhOriginFile.readWholeSDS("NDVI", pGlobal1km_NDVI); VhOriginFile.readWholeSDS("BT4", pGlobal1km_BT); VhOriginFile.readWholeSDS("QA", pGlobal1km_QA); char sds_names[LENGTH_SDSNAME]= {"NDVI;BT4;QA"}; //NDE file name: //GVH1000mND_r04c01_v2r02_j01_s202408260000000_e202409012359599_c202409020343417.nc // fnameNDE2 : part 2 of NDE file name : 'data/test/GVH1000mND' char fnameNDE1[LENGTH_FILENAME]=""; SETSTRING(LENGTH_FILENAME,fnameNDE1, outputFolder,"/","GVH1000mND"); // fnameNDE2 : part 2 of NDE file name : '_v2r02_j01_s202408260000000_e202409012359599_c202409020343417.nc' char fnameNDE2[LENGTH_FILENAME]=""; snprintf(fnameNDE2, LENGTH_FILENAME, "_v2r02_%s", satID); vh_AddTimeToFilename(fnameNDE2,7, year,week); ADDSTRINGS(LENGTH_FILENAME, filename, ".nc"); const int isNDEfilename =1; float dLat = (GVI_latlonRange.maxLat -GVI_latlonRange.minLat) / Global1km_NY; for (int row = 0; row < tilesY; row++) { for (int col = 0; col < tilesX; col++) { char tileTag[LENGTH_SDSNAME]; snprintf(tileTag,LENGTH_SDSNAME,"_r%02dc%02d",row+1,col+1); GeoRange latlonRange; latlonRange.maxLat = GVI_latlonRange.maxLat - row * NY * dLat ; latlonRange.minLat = GVI_latlonRange.maxLat - (row +1) * NY * dLat ; latlonRange.minLon = GVI_latlonRange.minLon + col * NX * dLat ; latlonRange.maxLon = GVI_latlonRange.minLon + (col+1) * NX * dLat ; if(isNDEfilename ==0) { //short file name: SETSTRING( LENGTH_FILENAME, pOutFile -> filename, filepath1,tileTag, ".nc"); } else { //NDE file name SETSTRING( LENGTH_FILENAME, pOutFile -> filename, fnameNDE1,tileTag, fnameNDE2); ADDSTRINGS( LENGTH_FILENAME, pOutFile -> filename, ".nc"); } int variables = pOutFile -> setOutputVaraibles(sds_names); printf("row=%d,col=%d, %s\n latlonRange=[%f,%f,%f,%f],dLat=%f, variables=%d\n", row,col,pOutFile -> filename, latlonRange.minLat,latlonRange.minLon,latlonRange.maxLat,latlonRange.maxLon,dLat,variables); char mapName[LENGTH_SDSNAME]=""; snprintf(mapName, LENGTH_SDSNAME, "ND_%04d_Week_%02d", year, week); pOutFile->InputFileList_initialize(LENGTH_FILENAME); pOutFile->InputFileList_add(getFileBaseName(filename)); pOutFile -> setSDS(variables, NX,NY, latlonRange, 1, mapName); /* start HDF file */ pOutFile -> open(pOutFile -> filename, ACCESS_WRITE); setWeeklyFileAttributes(pOutFile, satID, mapName, year, week, latlonRange); pOutFile -> createSDS(); pOutFile -> enddef(); pOutFile -> showVariableDescription(); for (int y = 0; y < NY; y++) { int pos1= (y+row * NY) * Global1km_NX+ col * NX; int pos2= y * NX; for (int x = 0; x < NX; x++) { pNDVI[ pos2 +x]= pGlobal1km_NDVI [pos1 + x]; pBT [ pos2 +x]= pGlobal1km_BT [pos1 + x]; pQA [ pos2 +x]= pGlobal1km_QA [pos1 + x]; } } //y pOutFile->save1SDS(&pOutFile->VARs[0], pNDVI); pOutFile->save1SDS(&pOutFile->VARs[1], pBT); pOutFile->save1SDS(&pOutFile->VARs[2], pQA); pOutFile->saveLatLonArray(); pOutFile->close(); }} // col, row VhOriginFile.close(); free1D(__FUNCTION__, (VOIDP &)pGlobal1km_NDVI, "pGlobal1km_NDVI"); free1D(__FUNCTION__, (VOIDP &)pGlobal1km_QA, "pGlobal1km_QA"); free1D(__FUNCTION__, (VOIDP &)pNDVI, "pNDVI"); free1D(__FUNCTION__, (VOIDP &)pQA, "pQA"); printStatus(__FUNCTION__, COMPLETED_NORMALLY, & time_begin); return (SUCCEED); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: getFileBaseName_NDE() * FUNCTION: get NDE style file name for a VH file * INPUT: * fileType: output file type options: "ND","FinalSM","InitialSM", "FinalVH","InitialVH", "Daily"] * year, period: output year and jday * OUTPUT: * filename: output file name -----------------------------------------------------------------------------*/ void CVHsuite::getFileBaseName_NDE(char * filename, const char *prefix, const char *fileType, const char *tileTag, int year, int period) { /* GVH1KM_FINAL_WEEKLY_SM GVH1KM_INITIAL_WEEKLY_SM GVH1KM_FINAL_WEEKLY_VH GVH1KM_INITIAL_WEEKLY_VH GVH1000mFinalSM_v2r02_j01 GVH1000mInitialVH_v2r02_npp GVH500mFinalSM_r03c03_v3r1_j02 GVH500mInitialSM_r03c03_v3r1_j02 GVH1KMND GVH1000mND_r01c03_v2r02_npp GVH500mND_r03c03_v3r1_j02 VH1KM_DAILY ViirsVH1000mDaily_r01c02_v2r02_j01 ViirsVH500mDaily_r03c03_v3r1_j02 if save as 1 tile, then tile tage will be ignored GVH500mFinalSM_v3r1_j02 GVH500mInitialSM_v3r1_j02 */ snprintf(filename, LENGTH_FILENAME, "%s%s%s%s_%s_%s", prefix, m_pConfig -> sResolution_NDE, fileType,tileTag, VHP_CODEVERSION, m_pConfig -> SatelliteID); int NdaysPerPeriod = (strcmp(fileType,"Daily") ==0) ? 1: 7; // add time: _syyyymmddHHMMSSS_eyyyymmddHHMMSSS_cyyyymmddHHMMSSS vh_AddTimeToFilename(filename, NdaysPerPeriod, year,period); printf("%s(),file=%s\n", __FUNCTION__, filename); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: getFileName_NDE() * FUNCTION: get NDE style file path for a VH file * INPUT: * * fileType: output file type options:"Weekly", "ND","SM","VH"] * year, period: output year and jday * OUTPUT: * filename: output file name -----------------------------------------------------------------------------*/ void CVHsuite::getFileName_NDE(char * filename, int fileType, int year, int period) { char *folder = nullptr; char sFileType[LENGTH_20] = ""; if(fileType == FILE_GVI) {stringCopy(LENGTH_20,sFileType, "Weekly"); folder= m_pConfig -> DIR_WeeklyComposite;} if(fileType == FILE_ND) {stringCopy(LENGTH_20,sFileType, "ND"); folder= m_pConfig -> DIR_ND;} if(fileType == FILE_SM) {stringCopy(LENGTH_20,sFileType, "SM"); folder= m_pConfig -> DIR_SM;} if(fileType == FILE_VH) {stringCopy(LENGTH_20,sFileType, "VH"); folder= m_pConfig -> DIR_VH;} if(fileType == FILE_DAILY) {stringCopy(LENGTH_20,sFileType, "Dailly"); folder= m_pConfig -> DIR_DailyMap;} if(folder == nullptr) ERROR_EXIT("FileType ID is not recognized:", fileType); //set output folder --- char dir[LENGTH_FILENAME] =""; stringCopy(LENGTH_FILENAME, dir, folder); if (strcmp(m_pConfig -> sActiveTile, "") != 0) { //if the map is saved in multi-files, use yyyyddd as sub-folder char yyyyddd[LENGTH_SDSNAME]; snprintf(yyyyddd, LENGTH_SDSNAME, "/%4d%03d", year, period); ADDSTRINGS(LENGTH_FILENAME, dir, yyyyddd); } mkdir(dir); char fname_NDE[LENGTH_FILENAME]; getFileBaseName_NDE(fname_NDE, m_pConfig -> Prefix, sFileType, m_pConfig -> sActiveTile, year, period); SETSTRING(LENGTH_FILENAME, filename, dir, "/", fname_NDE); ADDSTRINGS(LENGTH_FILENAME, filename, ".nc"); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: convert_ViirsVCI_to_AvhrrSMN() * FUNCTION: convert from VIIRS SMN to AVHRR equivalent SMN * NOTE: this function implements Wenze Yangs' algorithm * * INPUT: * file_ViirsVH: VIIRS SMN file name * file_ViirsClimat: VIIRS climatology * file_AvhrrClimat: AVHRRS climatology * fileOut_AVHRR: AVHRR equivalent SMN file ( output file name) -----------------------------------------------------------------------------*/ void CVHsuite::convert_ViirsVCI_to_AvhrrSMN(char *satID,int year,int week, char *file_ViirsVH, char *file_AvhrrClimat, char *fileOut_AVHRR){ timeval time_begin; gettimeofday( & time_begin, nullptr); printf("%s():\n", __FUNCTION__); printf("file_ViirsVH = %s\n", file_ViirsVH); printf("file_AvhrrClimat= %s\n", file_AvhrrClimat); printf("fileOut_AVHRR = %s\n", fileOut_AVHRR); /* file_A36_climat = "data_yang/AVHRR_VHP/4km/climate_v2017_long/VHP.G04.C07.CLIMAT.P"+sPeriod2 +".nc" file_V36_VH = "data/" + satID + "/4km/VH/VGVI.G04.C07." + satID + ".P"+syear+sPeriod1 +".VH.nc" file_A36_SM = "data_yang/AVHRR_VHP/4km/VH/VHP.G04.C07." + satID + ".P" +syear+sPeriod1 +".nc" climate_sds_types=["MEAN", "MAX","MIN","STDDEV"] vh_sds_types=["VCI", "TCI"] sm_sds_types=["SMN", "SMT"] */ //copy file attibute from the fist file name //copy the SDS attributes //loop for tiles and get the sample data, set to sampled data array CVhFile ncFileViirsVH; CVhFile ncFileAvhrrClimat; CVhFile ncFileOutput; CVhFile *pOutFile =& ncFileOutput; pOutFile->InputFileList_initialize(LENGTH_FILENAME); pOutFile->InputFileList_add(getFileBaseName(file_ViirsVH)); pOutFile->InputFileList_add(getFileBaseName(file_AvhrrClimat)); int num_datasets, num_global_attrs; if(ncFileViirsVH.open(file_ViirsVH, ACCESS_READ) == FAIL) return ; ncFileViirsVH.IO_FileInfo(&num_datasets, &num_global_attrs); printf("num_datasets=%d, num_global_attrs=%d\n", num_datasets, num_global_attrs); int NX= 0; int NY= 0; ncFileViirsVH.getSize2D("VCI", NX,NY); int NP=NX *NY; int16 *pBuffer = (int16 *) ALLOCATE( NP *4 * sizeof(int16)); GeoRange SDS_latlonRange; ncFileViirsVH.readSpatialRange( & SDS_latlonRange); char sds_names[LENGTH_SDSNAME]= {"SMN;SMT;QA"}; char mapName[LENGTH_SDSNAME]=""; snprintf(mapName, LENGTH_SDSNAME, "SM_%04d_Week_%02d", year, week); int variables = pOutFile -> setOutputVaraibles(sds_names); pOutFile -> setSDS(variables, NX,NY, SDS_latlonRange, 1, mapName); pOutFile -> open(fileOut_AVHRR, ACCESS_WRITE); setWeeklyFileAttributes(pOutFile, satID, mapName, year, week, SDS_latlonRange); //create SDS for output pOutFile -> createSDS(); pOutFile -> enddef(); pOutFile -> showVariableDescription(); //Adjust VIIRS SMN to AVHRR SMN; //assuming the scale_factor are the same for SMN, MNN and MXN const float scale_VCI= 0.01; const float scale_TCI= 0.01; const float scale_NDVI = 1.0E-3; const float scale_BT = 0.1; const float missing = -9999; const int16 int16Missing= -999; int16 *pVCI = pBuffer +NP*3; for(int iVar=0; iVar< 3; iVar++){ if(iVar ==0){ printf("read and convert VCI ==> SMN ---\n"); int16 *pMXN = pBuffer +NP; int16 *pMNN = pBuffer +NP*2; int16 *pSMN = pBuffer ; ncFileViirsVH.readWholeSDS("VCI", pVCI); ncFileAvhrrClimat.readWholeSDS(file_AvhrrClimat, "NDVI_MAX", pMXN); ncFileAvhrrClimat.readWholeSDS(file_AvhrrClimat, "NDVI_MIN", pMNN); // const float rangeNDVI [2] ={ -0.1,1}; for (int L = 0; L < NP; L++) { float minND = (pMNN[L] <= int16Missing)? missing : pMNN[L] * scale_NDVI ; float maxND = (pMXN[L] <= int16Missing)? missing : pMXN[L] * scale_NDVI ; float VCI = (pVCI[L] < 0 )? missing : pVCI[L] * scale_VCI ; float SMN =( VCI <0 )? minND : (VCI * (maxND- minND) /100 + minND) ; pSMN[L] = getScaledValue (&pOutFile->VARs[0], SMN); //pSMN[L] =( SMN < rangeNDVI[0] || SMN> rangeNDVI[1]) ? int16Missing : int(SMN /scale_NDVI +0.5) ; if(showMessage(LEVEL8) && L == 1776) { // for debugging printf(" VCI=%f, minND=%f, maxND=%f ==> SMN=%f,pSMN[L] =%d\n", VCI,minND, maxND,SMN, pSMN[L] ); } } //L }else if(iVar ==1){ printf("read and convert TCI ==> SMT ---\n"); int16 *pMXT = pBuffer +NP; int16 *pMNT = pBuffer +NP*2; int16 *pSMT = pBuffer; int16 *pTCI = pBuffer; ncFileViirsVH.readWholeSDS("TCI", pTCI); ncFileAvhrrClimat.readWholeSDS(file_AvhrrClimat, "BT_MAX", pMXT); ncFileAvhrrClimat.readWholeSDS(file_AvhrrClimat, "BT_MIN", pMNT); //const float rangeBT [2] ={ 0,400.0}; for (int L = 0; L < NP; L++) { float minT = (pMNT[L] <= int16Missing)? missing : pMNT[L] * scale_BT ; float maxT = (pMXT[L] <= int16Missing)? missing : pMXT[L] * scale_BT ; float VCI = ( pVCI[L] < 0 )? missing : pVCI[L] * scale_VCI ; float TCI = ( pTCI[L] < 0 )? missing : pTCI[L] * scale_TCI ; float T = (maxT - TCI * (maxT- minT) /100); float SMT = ( VCI < 0) ? minT : T; pSMT[L] = getScaledValue (&pOutFile->VARs[1], SMT); //pSMT[L] = ( SMT rangeBT [1]) ? int16Missing: int(SMT/scale_BT +0.5 ); if(showMessage(LEVEL8) && L == 3285430) { // for debugging printf("L=%d VCI=%f(%d), TCI=%f(%d), T= %f, SMT=%f minT,maxT=[%f,%f] ==> pSMT[L] =%d\n", L, VCI,pVCI[L], TCI,pTCI[L], T, SMT, minT, maxT, pSMT[L] ); } } //L }else if(iVar ==2){ printf("read and save 'QA'\n"); uint8 *pQA= (uint8 *) pBuffer; ncFileViirsVH.readWholeSDS("QA", pQA); } pOutFile->save1SDS(&pOutFile->VARs[iVar], pBuffer); } pOutFile->saveLatLonArray(); pOutFile->close(); ncFileViirsVH.close(); free1D(__FUNCTION__, (VOIDP &)pBuffer, "pBuffer"); printNameOfSavedFile(__FUNCTION__, fileOut_AVHRR, & time_begin); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: searchFileName_NDE() * FUNCTION: search NDE file for a VH file * short file name: File name is as the same style as AVHRR VHP, example as * VGVI500mDaily_r01c01*.nc_ for FILE_DAILY * VGVI500mWeekly_r01c01_*.nc for FILE_GVI * VGVI500mND_r01c01_*.nc for FILE_ND * VGVIG500mSM_r01c01_*.nc for FILE_SM * VGVIG500mVH_r01c01_*.nc for FILE_VH * * INPUT: * fileType: output file type * year, week: output year and week * OUTPUT: * filename_output: output file name -----------------------------------------------------------------------------*/ void CVHsuite::searchFileName_NDE(char * filename_output, int fileType, int year, int week) { char fileTag[LENGTH_SDSNAME], sFileType[LENGTH_SDSNAME]; char * folder = getFolder(fileType, fileTag, sFileType); char NewFolder[LENGTH_FILENAME] = ""; stringCopy(LENGTH_FILENAME, NewFolder, folder); char yyyywww[LENGTH_SDSNAME]; snprintf(yyyywww, LENGTH_SDSNAME, "%4d%03d", year, week); char * satID; int satNumber; satID = getSatelliteID(fileType, year, week, & satNumber); if (strlen(m_pConfig -> sActiveTile) > 0 ){ snprintf( NewFolder, LENGTH_FILENAME,"%s/%s/%4d/%s", folder, satID, year, yyyywww); //if(m_pConfig -> bRunWithPCF==1) SETSTRING(LENGTH_FILENAME, NewFolder, folder, "/", yyyywww); } printf( "m_pConfig -> bRunWithPCF=%d\n", m_pConfig -> bRunWithPCF); printf( "NewFolder=%s\n", NewFolder); char fname1[LENGTH_SDSNAME]; char fname2[LENGTH_SDSNAME]; char pattern[LENGTH_FILENAME]; int jday1, jday2; getJdayRangeForPeriod(week, & jday1, & jday2); int month, dayInMonth; getMonth(year, jday1, & month, & dayInMonth); if(fileType == FILE_DAILY) { int jday= week; getMonth(year, jday, & month, & dayInMonth); } char sFileTypeTag[LENGTH_SDSNAME]; stringCopy(LENGTH_SDSNAME, sFileTypeTag, ""); SETSTRING(LENGTH_SDSNAME, fname1, "*", m_pConfig -> sResolution_NDE,"" ); char yyyymmdd[LENGTH_20]; snprintf(yyyymmdd, LENGTH_20, "_s%4d%02d%02d*.nc", year, month, dayInMonth); SETSTRING(LENGTH_SDSNAME, fname2, "*", m_pConfig->sActiveTile,"*"); ADDSTRINGS3(LENGTH_SDSNAME, fname2, satID, yyyymmdd); if(fileType == FILE_GVI) SETSTRING(LENGTH_FILENAME,pattern, fname1, "Weekly",fname2); if(fileType == FILE_ND) SETSTRING(LENGTH_FILENAME,pattern, fname1, "ND", fname2); if(fileType == FILE_SM) SETSTRING(LENGTH_FILENAME,pattern, fname1, "SM", fname2); if(fileType == FILE_VH) SETSTRING(LENGTH_FILENAME,pattern, fname1, "VH", fname2); if(fileType == FILE_DAILY) SETSTRING(LENGTH_FILENAME,pattern, fname1, "Daily", fname2); char fullPattern[LENGTH_FILENAME]; SETSTRING(LENGTH_FILENAME, fullPattern, NewFolder, "/", pattern); getLastFileWithPattern(m_pConfig -> DIR_Temp, fullPattern, filename_output) ; printf("%s(): filename_output=%s\n",__FUNCTION__,filename_output ); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: VHP_getInputFileName() * FUNCTION: get VHP short file name, if the file do not exist, search NDE file for a VH file * NOTE: WG:20241210, When reading daily, weekly, ND files: first try fix the file name, then try to search NDE file * INPUT: * fileType: input file type * year: input year * week: week or 'day of year' * OUTPUT: * filename_output: output file name -----------------------------------------------------------------------------*/ void CVHsuite::VHP_getInputFileName(char * filename, int fileType, int year, int week) { getFilename_VHPstyle(filename, fileType, year, week); if(FileExist(filename) == 1) { return; }else if(fileType != FILE_CLIMAT){ printf("%s(), filename=%s does not exist, will search NDE style file\n", __FUNCTION__, filename); searchFileName_NDE(filename,fileType , year, week); } printf("%s(), filename=%s\n", __FUNCTION__, filename); } /*----------------------------------------------------------------------------- * FUNCTION_NAME: create_VHI_Files() * FUNCTION: compute VHI from SM and climatology * PROCEDURE: * - search input file names (SM and climatology file), * they can be provided by PCF file; if not, search specified folder to * get file withspecified satelliteID, file type, year and week * - get the adjustment coefficients for NDVI and BT * they are provided by configure file * - set buffers for intput files to read (setFileToRead) * - set for the output VH file (setOutputFile) * - create output file, * - create SDS in the output file * - write file attributes * - initialize buffer for metadata * - intiialieze buffers for SDS arries associated to input and output files * - In order reduce RAM, use samller buffer (1/8) for input files, (segments=8) * loop for each segment, * - read input SM and climatology data * - calculate VH * - set to full size array of VCI,TCI,VHI * - save VH file * - add result file names to PSF file * INPUT: * year, week: the year and week number of the week being processed * RETURN: * status of execution -----------------------------------------------------------------------------*/ int CVHsuite::create_VHI_Files(char * task, int year, int week) { timeval time_begin; gettimeofday( & time_begin, nullptr); const int needConvertToAvhrrSMN =0; // read PCF file, and put file names in an array of strings . //read other options float NDVI_toAVHRR_slope =1; float NDVI_toAVHRR_intercept=0; float BT_toAVHRR_slope=1; float BT_toAVHRR_intercept=0; if(needConvertToAvhrrSMN ==1){ sscanf(m_pConfig -> ScalingParameters, "%f, %f, %f, %f", & NDVI_toAVHRR_slope, & NDVI_toAVHRR_intercept, & BT_toAVHRR_slope, & BT_toAVHRR_intercept); printf(" ScalingParameters: NDVI slope=%f intercept=%f, BT slope=%f intercept=%f\n", NDVI_toAVHRR_slope, NDVI_toAVHRR_intercept, BT_toAVHRR_slope, BT_toAVHRR_intercept); } //input files: CVhFile fileSM, fileSM_QA, fileCLM; int input_segments = 1; enum FILE_SDSID_SM { ID_SMN, ID_SMT, ID_SMQA }; enum FILE_SDSID_CLIMAT { ID_MXN, ID_MNN, ID_MXT, ID_MNT }; enum FILE_SDSID_VH { ID_VCI, ID_TCI, ID_VHI, ID_VHQA }; if (setInputVHfile( & fileSM_QA, FILE_SM, year, week, "QA", input_segments) == FAIL) { printStatus(task, ERROR_ARCHIVE_NO_VHP_FILE); return FAIL; } if (setInputVHfile( & fileSM, FILE_SM, year, week, "SMN;SMT;", input_segments) == FAIL) { printStatus(task, ERROR_ARCHIVE_NO_VHP_FILE); return FAIL; } if (setInputVHfile( & fileCLM, FILE_CLIMAT, year, week, "NDVI_MAX;NDVI_MIN;BT_MAX;BT_MIN;", input_segments) == FAIL) { printStatus(task, WARNING_VHC_READ_CLIMAT_FILE); return FAIL; } //output file: CVhFile fileVH; fileVH.InputFileList_initialize(LENGTH_FILENAME * 2); fileVH.InputFileList_add(fileSM.filename); fileVH.InputFileList_add(fileCLM.filename); printf(" nx=%d ny=%d, ", fileSM.m_SDScolumns, fileSM.m_SDSrows); printf(" SDS minLat= %.3f, maxLat= %.3f, minLon= %.3f, maxLon= %.3f\n", fileSM.m_SDSlatlonRange.minLat,fileSM.m_SDSlatlonRange.maxLat, fileSM.m_SDSlatlonRange.minLon,fileSM.m_SDSlatlonRange.maxLon); printf(" input SM file=%s\n", fileSM.filename); printf(" input CLIMAT file=%s\n", fileCLM.filename); if (setOutputVHfile( & fileVH, FILE_VH, year, week, "VCI;TCI;VHI;QA;", & fileSM) == FAIL) return FAIL; //calculate memory for major variables if(showMessage(LEVEL4)) VHP_print_RAM_size("VH", fileVH.m_Segment.NX, fileVH.m_Segment.NY, input_segments); fileSM.processBuffer(ALLOCATE); fileSM_QA.processBuffer(ALLOCATE); fileCLM.processBuffer(ALLOCATE); fileVH.processBuffer(ALLOCATE); fileVH.processBuffer(INITIALIZE); // uint8 * pSM_QA = (uint8 * ) fileSM.Buffers[ID_SMQA]; uint8 * pSM_QA = (uint8 * ) fileSM_QA.Buffers[0]; uint8 * pVH_QA = (uint8 * ) fileVH.Buffers[ID_VHQA]; struct_Segment * pSegment = & fileSM.m_Segment; float VCI, TCI, VHI; int test = 0; //printf(" processing input files by %d segements\n", input_segments); for (int iSegment = 0; iSegment < input_segments; iSegment++) { if (input_segments > 1) printf(" process iSegment=%d (of %d)\n", iSegment + 1, input_segments); fileSM.setSegment(iSegment); fileSM_QA.setSegment(iSegment); fileCLM.setSegment(iSegment); int NX= pSegment -> NX; int K0 = NX * pSegment -> NY * iSegment; fileSM_QA.processBuffer(OPEN_AND_READ); for (long i = 0; i < pSegment -> ValidPoints; i++) { pVH_QA[ i + K0 ] = pSM_QA[i]; } Box box; int NP_valid= getYrangeWithLandPixels(pSM_QA, NX, pSegment -> ValidLines, m_pConfig -> sActiveTile, &box); if(NP_valid == 0) { printf(" skip reading SM and Climate for this segment, \n" ); continue; } // read a small section ( part of the segment), from line y1 to y2 in this segment int32 count [2]= {box.ny, box.nx}; int32 start [2]= {box.y1 + iSegment * pSegment -> NY ,box.x1}; int32 stride[2]= {1,1}; printf(" start=[%d,%d], count=[%d,%d]\n", start[0], start[1], count[0], count[1]); fileCLM.open(fileCLM.filename, ACCESS_READ); fileCLM.IO_readSDS("NDVI_MAX", start, stride, count, fileCLM.Buffers[0]); fileCLM.IO_readSDS("NDVI_MIN", start, stride, count, fileCLM.Buffers[1]); fileCLM.IO_readSDS("BT_MAX", start, stride, count, fileCLM.Buffers[2]); fileCLM.IO_readSDS("BT_MIN", start, stride, count, fileCLM.Buffers[3]); fileCLM.close(); fileSM.open(fileSM.filename, ACCESS_READ); fileSM.IO_readSDS("SMN", start, stride, count, fileSM.Buffers[0]); fileSM.IO_readSDS("SMT", start, stride, count, fileSM.Buffers[1]); fileSM.close(); int np= box.nx * box.ny; for (long i = 0; i < np; i++) { int y = i / box.nx; int L =K0 + (box.y1 + y) * NX + box.x1 + (i % box.nx); if ((pVH_QA[L] & 12) == 0) continue; //skip water pixels float SMN = fileSM.getValue( ID_SMN, i); float MXN = fileCLM.getValue(ID_MXN, i); float MNN = fileCLM.getValue(ID_MNN, i); float SMT = fileSM.getValue( ID_SMT, i); float MXT = fileCLM.getValue(ID_MXT, i); float MNT = fileCLM.getValue(ID_MNT, i); //convert VIIRS SMN (SMT) to AVHRR equivalent NDVI (BT4) if(needConvertToAvhrrSMN == 1) { SMN = SMN * NDVI_toAVHRR_slope + NDVI_toAVHRR_intercept; SMT = SMT * BT_toAVHRR_slope + BT_toAVHRR_intercept; } VCI = calculate_VCI(SMN, MXN, MNN); TCI = calculate_TCI(SMT, MXT, MNT); VHI = calculate_VHI(VCI, TCI); if (showMessage(LEVEL8) && SMN > 0 && test < 10) { printf(" %ld, SMN=%f, MAX= %f, MIN=%f, VCI=%f\n", i, SMN, MXN, MNN, VCI); test++; } fileVH.setValue(ID_VCI, L, VCI); fileVH.setValue(ID_TCI, L, TCI); fileVH.setValue(ID_VHI, L, VHI); } } fileVH.saveVH_2NC( pVH_QA ); fileSM.processBuffer(DEALLOCATE); fileSM_QA.processBuffer(DEALLOCATE); fileCLM.processBuffer(DEALLOCATE); fileVH.processBuffer(DEALLOCATE); printStatus(task, COMPLETED_NORMALLY, & time_begin); return (SUCCEED); }