diff --git a/docs/Release_Notes/README.md b/docs/Release_Notes/README.md index b290b7d..0c71227 100755 --- a/docs/Release_Notes/README.md +++ b/docs/Release_Notes/README.md @@ -2,12 +2,16 @@ Institute for the Environment, UNC-Chapel Hill - - The Spatial Allocator is now being distributed through GitHub (https://github.com/CMASCenter/Spatial-Allocator). * Releases available on github: 4.3, 4.3.1, 4.3.2, and 4.4. * Releases are also available from the [CMAS Center](http://www.cmascenter.org). +### 02/14/2020 (last update at EPA by Limei Ran) +* Updated land use tools for MODIS IGBP land cover name change in MODIS collection 6, +* WRF/CMAQ to EPIC tool update for CMAQ5.3 - added species and get rid of orgnaic Ndep output +* EPIC to CMAQ tool update to handle leap year day +* Utility tool updates + ### 06/2019 Release Notes for Spatial Allocator v4.4 * Includes the new Postgres Surrogate Tool, for generating surrogates using a PostgreSQL database with the PostGIS extension diff --git a/src/raster/README b/src/raster/README index 51069a6..fbb8edb 100755 --- a/src/raster/README +++ b/src/raster/README @@ -1,46 +1,163 @@ -10/2018 -SA 4.3.2 release -updated Raster tools for EPIC +Spatial Allocator Tools on a Linux Server +L.Ran and D. Yang, EPA and IE UNC-Chapel Hill - extractEPIC2CMAQ.cpp - added additional variable FBARE: Bare Land Fraction for Wind Erosion - computeGridGOES.cpp - removed hardcoded dimenstions so that program can be used to read new or old versions of GEOS Imager data. The coverage was extended in the new retrievals (from 1500x800 to 1552x1300). - computeSiteDailyWeather.cpp - increase the number of wet and dry deposition variables available from CMAQ - compute_EPICSiteData.cpp - Changed how a crop_acre is defined: crop_acre > cropArea_threshold && waterP < 50. +02/14/2020 (last update at EPA by Limei Ran) +Updated: + Raster Tools: + (1) Updated land use tools for MODIS IGBP land cover name change in MODIS collection 6, + (2) WRF/CMAQ to EPIC tool update for CMAQ5.3 - added species and get rid of orgnaic Ndep output + (3) EPIC to CMAQ tool update to handle leap year day + (4) Utility tool updates -09/2015 -SA 4.2 update release - 1. Modified the extraction tools for FEST-C. +06/2019 Release Notes for Spatial Allocator v4.4 -05/2014 -SA 4.2 release - 1. Updated computeGridLandUse_beld4.cpp for 2001 BELD4 generation - 2. Updated extractEPICYearlyAverage2CMAQ.cpp for crop area weighted variables + (1) Includes the new Postgres Surrogate Tool, for generating surrogates using a PostgreSQL database with the PostGIS extension +10/2018 Update Release Notes for SA 4.3.2 -05/2012 -1. preProcessNLCD.exe - (1) Images have to be GByte format. - (2) Created no-overlapping images in ESRI BIL format - EHdr -- ESRI .hdr Labelled inlimited size. + (1) Updated FEST-C tools for the new release of FEST-C v1.4 + (2) Updated GOES data processing tool for input directory + (3) Updated all geoprocessing classes based on the change of GDAL library + (4) Added new land use processing tool with MODIS vegetation product + (5) Updated Raster tools for EPIC + extractEPIC2CMAQ.cpp - added additional variable FBARE: Bare Land Fraction for Wind Erosion + computeGridGOES.cpp - removed hardcoded dimenstions so that program can be used to read new or old versions of GEOS Imager data. The coverage was extended in the new retrievals (from 1500x800 to 1552x1300). + computeSiteDailyWeather.cpp - increase the number of wet and dry deposition variables available from CMAQ + compute_EPICSiteData.cpp - Changed how a crop_acre is defined: crop_acre > cropArea_threshold && waterP < 50. -2. gl_sin -- global MODIS data. If you convert it to NLCD Albers projection. It shift to NW around 25km. - Extract out NA area from gl_sin, then project it to Laea and then to NLCD aea. It got rid of shift. - Use projected MODIS data in NLCD for North American and check the shifting. - Defined the data projection by (due to error in reading projection ) - ../src/libs/gdal-1.5.2/local/bin/gdal_translate -of EHdr -a_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" na_modis.bil na_modis_nlcd.bil +06/2017 Update Release Notes for SA 4.3.1 -3 Processed NOAA wa_01 image by ../src/libs/gdal-1.5.2/local/bin/gdal_translate -a_ullr -2144475 3177765 -1736265 2777085 -outsize 13607 13356 wa_01.img wa_01_regrided.img + Raster Tools: + (1) Added input data files for the FESTC release, includes crop fraction data for US and Canada for 42 different crops for the year 2011 which is used for 2011 BELD4 data generation within FEST-C. + (2) Updated executables and source code for extractEPIC2CMAQ and extractEPICYearlyAverage2CMAQ + (3) (Minor modification of the source code/executables for the averaging and the extraction of EPIC output for CMAQ) -4. toNLCDRaster program needs at least one image file in DATADIR to obtain the standard NLCD projection info. +02/2017 Update Release Notes for SA 4.3 -5. src/libs/gdal-1.5.2/local/bin/gdal_rasterize -a GRIDID -l temp_grdshape_nlcd.shp ../output/wrf12km_nc_30m.bil can not take input - shapefile in other directory (ONLY CURRENT). + Vector Tools: + (1) Update Earth Radius to 6370000 + (2) Correct bug in Polar Stereographic grid description format + (3) Added Beld4SMK tool to the release + (4) Added 64-bit executables + (5) Update reference data -6. test test_geos.cpp -g++ -I/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/gdal-1.9.1/local/include -I/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/geos-3.3.6/local/include test_geos.cpp -o test_geos -L/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/gdal-1.9.1/local/lib -lgdal -L/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/geos-3.3.6/local/lib -lgeos_c + Documentation + (1) All documentation ported to GitHub Markdown +### 05/31/2016 Update notes for SA 4.2: -g++ -I/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/gdal-1.9.1/local/include -I/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/geos-3.3.6/local/include -I/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/proj-4.8.0/local/include test_geos.cpp -o test_geos -static -pthread -L/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/gdal-1.9.1/local/lib -lgdal -L/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/geos-3.3.6/local/lib -lgeos_c -lgeos -L/nas01/depts/ie/cempd/SA/sa_06_2009/src/libs/proj-4.8.0/local/lib -lproj -ldl + Raster Tools: + (1) EPIC site info tool update for elevation output + Vector Tools: + (1) Update projection radius to WRF 6370000 + (2) Bug fix for makefile +02/05/2016 +Update Update Release Notes for SA 4.2: -http://trac.osgeo.org/geos/browser/trunk/tests/unit/geom/ + Raster Tools: + (1) Updated elevation output format for EPIC site information + (2) Modified EPIC extraction tools for added EPIC output variables + + Vector Tools: + (1) Fixed an error in polygon surrogate computation + + +09/30/2015 + +Update Release Notes for SA 4.2: + This update release has the following updates: + + Vector Tools: + (1) Fixed an error in re-gridding the BELD3 data for CMAQ. + (2) Fixed an error in surroagte computation related to holes in polyogns. + + Raster Tools: + (1) Updated extraction tools for FEST-C output because of changed format and + variable name definitions. + + Surrogate Tools: + (1) Fixed a bug related with polygon surrogate normaliztion when st+county code is 4 digit. + (2) Fixed a bug in surrogate normalization in reading surrogate data lines with leading spaces and no "!" sign. + (3) Added QA threshold for reporting. + (4) Updated the surrogate Merging Tool to output QA information on the source surrogates. + (5) Added leading "0" for one digit state codes in both merging and normalization tools. + + + +1. Set the system: + (1) Modify ./bin/sa_setup.csh + + (2) Include sa_setup.csh in your .cshrc + + (3) Troubleshooting for library errors in the SA Raster Tools: + >Recompile all libraries under src/libs following instructions in the src/libs/README file. + >Modiy src/raster/Makefile for correct paths. + >Type "make clean" to clean previous compiled programs. + >Type "make" to compile the tools. + >Type "make -B install" or "make install" to install compiled tools. + + (4) Troubleshooting for library errors in the SA Vector Tools: + >Modiy src/vector/Makefile for correct library paths. Users can use src/vector/libs_32bits + >Type "make clean" to clean previous compiled programs. + >Type "make" to compile the tools. + >Type "make install" to install compiled tools. + + +2. Vector tools are stored in: bin/32bits and sample vector script files are in sa_052014/scripts. + Users normally do not need to recompile the Vector Tools as they are statically compiled. + The Raster Tools often need to be re-compiled as they are not statically built. + + +3. Important Notes: + tmp*.* files created under ./raster_script can be deleted after the completion of the run + + +4. Sample land use and satellite data can be stored in ./data/sat. See ./data/sat/README for data download. + Sample Raster Tools running script files in ./raster_scripts: + + +5. SA Raster Tools for Satellite data and land use data processing: + +(1) allocateGOES2WRFGrids.csh -- GOES data processing tool. Modify it based on your inputs. + Download GOES data and put all data under one directory which is defined in the script. + Under this directory, each day is a subdirectory as YYYYMMDD + which contains GOES files for that day. + Contact: Dr. Arastoo Pour Biazar (biazar@nsstc.uah.edu) at + UAH's Earth System Science Center (ESSC) for the GOES data which can be processed by the tool. + + +(2) allocateMODISL2CloudVars2Grids.csh -- MODIS L2 cloud and aerosol product processing tool. Modify it for your inputs. + Down load data from the NASA LAADS web site: http://ladsweb.nascom.nasa.gov/data/search.html + + MODIS Cloud Products: + Cloud product variables contain 5km and 1km data. To use this re-gridding tool, + users have to download two sets data: MOD06_L2 and MOD03 (Level 1 Geolocation 1km) into the input directory. + for Terra or MYD06_L2 and MYD03 for Aqua in input data directory. + + MODIS Aerosol Products: + Aerosol products contain variable data at 10km resolution (nadir). Users have + to download MOD04 for Terra or MYD04 for Aqua in the input data directory. + + +(3) allocateOMIL2vars2Grids.csh - OMI L2 aerosol and NO2 product processing tool. Modify it for your inputs. + Download data from the NASA mirador site: + http://mirador.gsfc.nasa.gov/cgi-bin/mirador/presentNavigation.pl?tree=project&project=OMI + Put downloaded data in one directory which is defined in the script file. + + +(4) allocateOMIvar2Grids.csh - OMI L2G and L3 aerosol and NO3 processing tool. Modify it for your inputs. + Download data from the NASA Geonanni web site: + http://gdata1.sci.gsfc.nasa.gov/daac-bin/G3/gui.cgi?instance_id=omi + The tool is also may also be used process MODIS L3 products (not tested well) from NASA MODIS web site. + + +(5) NLCD_MODIS_processor.csh - generate WRF grid NLCD and MODIS landuse data. + Download NLCD data from: http://www.mrlc.gov/nlcd06_data.php + Download MODIS tiled data (MCD12Q1) from: https://lpdaac.usgs.gov/products/modis_products_table + +(6) landuseTool_WRFCMAQ_BELD4.csh - generate BELD4 data from 2001 or 2006 NLCD/MODIS and crop tables. + Processed crop and tree tables and shapefiles are stored under "data" directory. + Download NLCD data from: http://www.mrlc.gov/nlcd06_data.php + Download MODIS tiled data (MCD12Q1) from: https://lpdaac.usgs.gov/products/modis_products_table diff --git a/src/raster/computeGridGOES.exe b/src/raster/computeGridGOES.exe index 35321ba..fe414a0 100755 Binary files a/src/raster/computeGridGOES.exe and b/src/raster/computeGridGOES.exe differ diff --git a/src/raster/computeGridLandUse.cpp b/src/raster/computeGridLandUse.cpp index 283b1e5..de126ae 100755 --- a/src/raster/computeGridLandUse.cpp +++ b/src/raster/computeGridLandUse.cpp @@ -93,7 +93,9 @@ GByte MODIS_NODATAVALUE = 255; //245: unclassified water and 0,17: water //MODIS tile land cover class variable: may change depending on NASA MODIS LC products -string modisVarName = string ( "Land_Cover_Type_1" ); +string modisVarName = string ( "Land_Cover_Type_1" ); //IGBP NASA type +string modisVarName_v6 = string ( "LC_Type1" ); + //match NLCD classes by array index: 0 or 127 is background @@ -464,6 +466,12 @@ int main( int nArgc, char* papszArgv[] ) //Get MODIS land cover file info imageFile = modisFiles.at(0); //use the first MODIS image files + + if ( imageFile.find( "MCD12Q1" ) != string::npos && imageFile.find( ".006." ) != string::npos ) + { + modisVarName = modisVarName_v6; + } + imageInfo = getHDF4VarInfo (imageFile, modisVarName); } else diff --git a/src/raster/computeGridLandUse.exe b/src/raster/computeGridLandUse.exe index c5bd4f9..3752724 100755 Binary files a/src/raster/computeGridLandUse.exe and b/src/raster/computeGridLandUse.exe differ diff --git a/src/raster/computeGridLandUse_LAI_MODIS.cpp b/src/raster/computeGridLandUse_LAI_MODIS.cpp index 77a3676..4da3103 100644 --- a/src/raster/computeGridLandUse_LAI_MODIS.cpp +++ b/src/raster/computeGridLandUse_LAI_MODIS.cpp @@ -98,6 +98,7 @@ GByte MODISLAI_NODATAVALUE = 249; //249-255 not LAI and PAR values //MODIS tile land cover class variable: may change depending on NASA MODIS LC products string modisVarName = string ( "Land_Cover_Type_1" ); +string modisVarName_v6 = string ( "LC_Type1" ); //match NLCD classes by array index: 0 or 127 is background @@ -528,6 +529,12 @@ int main( int nArgc, char* papszArgv[] ) //Get MODIS land cover file info imageFile = modisFiles.at(0); //use the first MODIS image files + + if ( imageFile.find( "MCD12Q1" ) != string::npos && imageFile.find( ".006." ) != string::npos ) + { + modisVarName = modisVarName_v6; + } + imageInfo = getHDF4VarInfo (imageFile, modisVarName); } else diff --git a/src/raster/computeGridLandUse_LAI_MODIS.exe b/src/raster/computeGridLandUse_LAI_MODIS.exe index 04a5c8f..a7de837 100755 Binary files a/src/raster/computeGridLandUse_LAI_MODIS.exe and b/src/raster/computeGridLandUse_LAI_MODIS.exe differ diff --git a/src/raster/computeGridLandUse_beld4.cpp b/src/raster/computeGridLandUse_beld4.cpp index 786e21b..6bf2fca 100755 --- a/src/raster/computeGridLandUse_beld4.cpp +++ b/src/raster/computeGridLandUse_beld4.cpp @@ -115,6 +115,7 @@ GByte MODIS_NODATAVALUE = 255; //245: unclassified water and 0,17: water //MODIS tile land cover class variable: may change depending on NASA MODIS LC products string modisVarName = string ( "Land_Cover_Type_1" ); //IGBP NASA type +string modisVarName_v6 = string ( "LC_Type1" ); //match NLCD classes by array index: 0 or 127 is background @@ -579,6 +580,12 @@ int main( int nArgc, char* papszArgv[] ) //Get MODIS land cover file info imageFile = modisFiles.at(0); //use the first MODIS image files + + if ( imageFile.find( "MCD12Q1" ) != string::npos && imageFile.find( ".006." ) != string::npos ) + { + modisVarName = modisVarName_v6; + } + imageInfo = getHDF4VarInfo (imageFile, modisVarName); } else diff --git a/src/raster/computeGridLandUse_beld4.exe b/src/raster/computeGridLandUse_beld4.exe index d48bf99..081ba9e 100755 Binary files a/src/raster/computeGridLandUse_beld4.exe and b/src/raster/computeGridLandUse_beld4.exe differ diff --git a/src/raster/computeGridMODISL2Clouds.exe b/src/raster/computeGridMODISL2Clouds.exe index 907eea4..8c89d70 100755 Binary files a/src/raster/computeGridMODISL2Clouds.exe and b/src/raster/computeGridMODISL2Clouds.exe differ diff --git a/src/raster/computeGridOMI.exe b/src/raster/computeGridOMI.exe index 32562d0..80d4d51 100755 Binary files a/src/raster/computeGridOMI.exe and b/src/raster/computeGridOMI.exe differ diff --git a/src/raster/computeGridOMIL2.exe b/src/raster/computeGridOMIL2.exe index 3fdc3f7..6292312 100755 Binary files a/src/raster/computeGridOMIL2.exe and b/src/raster/computeGridOMIL2.exe differ diff --git a/src/raster/computeSiteDailyWeather.cpp_cmaq52 b/src/raster/computeSiteDailyWeather.cpp_cmaq52 index a4fd186..109098a 100755 --- a/src/raster/computeSiteDailyWeather.cpp_cmaq52 +++ b/src/raster/computeSiteDailyWeather.cpp_cmaq52 @@ -131,14 +131,14 @@ const char *dryNDName = "DNDEP"; const char *wetNDName = "WNDEP"; //CMAQ dry deposition variables -const int numDryDVars = 17; //15 CMAQ variables related to dry N -std::string cmaqDryDVarNames[] = {"NO2","NO","HNO3","ANO3I","ANO3J","ANO3K","PAN","PANX","NTR1","NTR2","INTR","N2O5","HONO","ANH4I","ANH4J","ANH4K","NH3"}; -double cmaqDryDVarFactors[] = {0.30435, 0.46667, 0.22222, 0.22581, 0.22581, 0.22581, 0.11570, 0.11570, 0.11696, 0.10366, 0.09519, 0.25926, 0.29787,1.00000, 1.00000, 1.00000, 1.05900}; +const int numDryDVars = 18; //18 CMAQ variables related to dry N +std::string cmaqDryDVarNames[] = {"NO2","NO","HNO3","ANO3I","ANO3J","ANO3K","PAN","PANX","OPAN","NTR1","NTR2","INTR","N2O5","HONO","ANH4I","ANH4J","ANH4K","NH3_DDEP"}; +double cmaqDryDVarFactors[] = {0.30435, 0.46667, 0.22222, 0.22581, 0.22581, 0.22581, 0.11570, 0.11570, 0.11570, 0.11696, 0.10366, 0.09519, 0.25926, 0.29787,1.00000, 1.00000, 1.00000, 1.05900}; //CMAQ wet deposition variables -const int numWetDVars = 18; //16 CMAQ variables related to wet N -std::string cmaqWetDVarNames[] = {"NO2","NO","ANO3I","ANO3J","ANO3K","HNO3","PAN","PANX","NTR1","NTR2","INTR","N2O5","HONO","PNA","ANH4I","ANH4J","ANH4K","NH3"}; -double cmaqWetDVarFactors[] = {0.30435, 0.46667, 1.00000, 1.00000, 1.00000, 0.98400, 0.11570, 0.11570, 0.11696, 0.10366, 0.09519, 0.25926, 0.29787, 0.177720, 1.00000, 1.00000, 1.00000, 1.05900}; +const int numWetDVars = 19; //19 CMAQ variables related to wet N +std::string cmaqWetDVarNames[] = {"NO2","NO","ANO3I","ANO3J","ANO3K","HNO3","PAN","PANX","OPAN", "NTR1","NTR2","INTR","N2O5","HONO","PNA","ANH4I","ANH4J","ANH4K","NH3"}; +double cmaqWetDVarFactors[] = {0.30435, 0.46667, 1.00000, 1.00000, 1.00000, 0.98400, 0.11570, 0.11570, 0.11570, 0.11696, 0.10366, 0.09519, 0.25926, 0.29787, 0.177720, 1.00000, 1.00000, 1.00000, 1.05900}; const double sigmaLevels[] = {0.50, 1.0}; @@ -1316,7 +1316,7 @@ void readCMAQDryDFile ( gridInfo grid, string cmaqDryDFile, int julianDateInt, s int ncStatus; ncStatus = nc_inq_varid (ncid, cmaqDryDVarNames[i].c_str(), &var_id); - //handle last variable: NH3_Dep or NH3 + //handle last variable: NH3_DDEP or NH3_Dep, NH3 if ( ncStatus != NC_NOERR && i != numDryDVars-1 ) { printf ("\tError: getting netCDF variable ID for variable: %s.\n", cmaqDryDVarNames[i].c_str() ); @@ -1325,15 +1325,23 @@ void readCMAQDryDFile ( gridInfo grid, string cmaqDryDFile, int julianDateInt, s if ( ncStatus != NC_NOERR && i == numDryDVars-1 ) { - //try NH3 - string diffNH3Str = string ( "NH3" ); + //try NH3_Dep + string diffNH3Str = string ( "NH3_Dep" ); cmaqDryDVarNames[i] = diffNH3Str; ncStatus = nc_inq_varid (ncid, cmaqDryDVarNames[i].c_str(), &var_id); if ( ncStatus != NC_NOERR ) { - printf ("\tError: getting netCDF variable ID for variable: %s.\n", cmaqDryDVarNames[i].c_str() ); - exit ( 1 ); + //try NH3 + diffNH3Str = string ( "NH3" ); + cmaqDryDVarNames[i] = diffNH3Str; + + ncStatus = nc_inq_varid (ncid, cmaqDryDVarNames[i].c_str(), &var_id); + if ( ncStatus != NC_NOERR ) + { + printf ("\tError: getting netCDF variable ID for variable: %s.\n", cmaqDryDVarNames[i].c_str() ); + exit ( 1 ); + } } } @@ -1672,12 +1680,14 @@ double cmaqDryDVarFactors[] = {0.30435, 0.46667, 0.22222, 0.22581, 0.22581, 0. varItemPos = 7; double ddepNHX = 0.0; - for ( k=0; k<=12; k++) + for ( k=0; k<=13; k++) { + if ( mcipVars[k][index] < 0.0 ) mcipVars[k][index] = 0.0; oNDep += mcipVars[k][index] * cmaqDryDVarFactors[k]; } - for (k=13; k<=16; k++) + for (k=14; k<=17; k++) { + if ( mcipVars[k][index] < 0.0 ) mcipVars[k][index] = 0.0; ddepNHX += mcipVars[k][index] * cmaqDryDVarFactors[k]; } rNDep += ddepNHX * 0.77778; @@ -1708,12 +1718,12 @@ double cmaqWetDVarFactors[] = {0.30435, 0.46667, 1.00000, 1.00000, 1.00000, 0. } oNDep += wdepTNO3 * 0.22581; - for ( k=6; k<=13; k++) + for ( k=6; k<=14; k++) { oNDep += mcipVars[k][index] * cmaqWetDVarFactors[k]; } - for ( k=14; k<=17; k++) + for ( k=15; k<=18; k++) { wdepNHX += mcipVars[k][index] * cmaqWetDVarFactors[k]; } @@ -1757,7 +1767,8 @@ double cmaqWetDVarFactors[] = {0.30435, 0.46667, 1.00000, 1.00000, 1.00000, 0. maxDep = max (maxDep, epicValue); if ( strcmp (mcipVarName, wetNDName) == 0 ) { - gNDep = 0.33 * epicValue; + //gNDep = 0.33 * epicValue; //get grid it for CMAQ5.3 + gNDep = 0.0; } } diff --git a/src/raster/computeSiteDailyWeather.exe b/src/raster/computeSiteDailyWeather.exe index 7f7694b..2ba08b9 100755 Binary files a/src/raster/computeSiteDailyWeather.exe and b/src/raster/computeSiteDailyWeather.exe differ diff --git a/src/raster/compute_EPICSiteData.exe b/src/raster/compute_EPICSiteData.exe index 271ba17..c190f05 100755 Binary files a/src/raster/compute_EPICSiteData.exe and b/src/raster/compute_EPICSiteData.exe differ diff --git a/src/raster/create_gridPolygon.exe b/src/raster/create_gridPolygon.exe index 6538831..1663a81 100755 Binary files a/src/raster/create_gridPolygon.exe and b/src/raster/create_gridPolygon.exe differ diff --git a/src/raster/extractEPIC2CMAQ.cpp b/src/raster/extractEPIC2CMAQ.cpp index 104e758..141e5be 100644 --- a/src/raster/extractEPIC2CMAQ.cpp +++ b/src/raster/extractEPIC2CMAQ.cpp @@ -612,6 +612,7 @@ int main(int nArgc, char *argv[]) //write output by day for each file string currentDay = stepTimeStr; + for (i=0; i #include #include +#include +#include +#include + #include "geotools.h" #include "commontools.h" @@ -95,7 +103,7 @@ int MODIS_ALB_PARAMS = 3; //3 albedo MCD43A1 variable parameters /****************************************/ string createGridShapes ( gridInfo grid ) { - string shapeName; + string shapeName, polyName; GDALDriver *poDriver = NULL; GDALDataset *poDS = NULL; @@ -108,7 +116,7 @@ string createGridShapes ( gridInfo grid ) double xmin,xmax; double ymin,ymax; - int i,j; + int i,j,k,m; int gid = 0; char *pszProj4; @@ -116,7 +124,7 @@ string createGridShapes ( gridInfo grid ) printf("\nGenerating grid shapefile for a given domain...\n"); //get Shapefile Name - if ( grid.name.empty() ) + if ( grid.name.empty() || grid.name.find(".shp") == string::npos ) { //create temp Shapefile name @@ -134,6 +142,23 @@ string createGridShapes ( gridInfo grid ) } + // for polygon text file + ifstream imageStream; //stream to read input polygon text file + if ( (! grid.polyID.empty()) && grid.polyID.find("GRIDID") == string::npos ) + { + polyName = grid.polyID; + + imageStream.open( polyName.c_str() ); + + if (! imageStream.good() ) + { + printf( "\tError: opening text file stream - %s\n", polyName.c_str() ); + exit ( 1 ); + } + } + + + /* -------------------------------------------------------------------- */ /* Register format(s). */ /* -------------------------------------------------------------------- */ @@ -188,7 +213,7 @@ string createGridShapes ( gridInfo grid ) /* Create Polygon Attribute Name. */ /* -------------------------------------------------------------------- */ - OGRFieldDefn oField( grid.polyID.c_str(), OFTInteger); + OGRFieldDefn oField( "GRIDID", OFTInteger); printf ("\tCreated item.\n"); if( poLayer->CreateField( &oField ) != OGRERR_NONE ) @@ -201,52 +226,345 @@ string createGridShapes ( gridInfo grid ) /* Assign each polygon attribute and geometry. */ /* -------------------------------------------------------------------- */ - for ( i = 0; i < grid.rows; i++ ) + if ( grid.polyID.empty() || grid.polyID.compare("GRIDID") == 0) { - ymin = grid.ymin + i * grid.yCellSize; - ymax = ymin + grid.yCellSize; + for ( i = 0; i < grid.rows; i++ ) + { + ymin = grid.ymin + i * grid.yCellSize; + ymax = ymin + grid.yCellSize; - for( j = 0; j < grid.cols; j++ ) - { - //printf ("\trow=%d col=%d\n",i,j); + for( j = 0; j < grid.cols; j++ ) + { + //printf ("\trow=%d col=%d\n",i,j); - gid ++; + gid ++; - poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() ); - poFeature->SetField( grid.polyID.c_str(), gid ); + poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() ); + poFeature->SetField( "GRIDID", gid ); - xmin = grid.xmin + j * grid.xCellSize; - xmax = xmin + grid.xCellSize; + xmin = grid.xmin + j * grid.xCellSize; + xmax = xmin + grid.xCellSize; - oRing.addPoint( xmin, ymin ); - oRing.addPoint( xmin, ymax ); - oRing.addPoint( xmax, ymax ); - oRing.addPoint( xmax, ymin ); - oRing.addPoint( xmin, ymin ); - //printf("ID = %d x: %lf %lf, y: %lf %lf\n",gid,xmin,xmax,ymin,ymax); + oRing.addPoint( xmin, ymin ); + oRing.addPoint( xmin, ymax ); + oRing.addPoint( xmax, ymax ); + oRing.addPoint( xmax, ymin ); + oRing.addPoint( xmin, ymin ); + //printf("ID = %d x: %lf %lf, y: %lf %lf\n",gid,xmin,xmax,ymin,ymax); - oPoly.addRing( &oRing ); + oPoly.addRing( &oRing ); - if ( poFeature->SetGeometry( &oPoly ) != OGRERR_NONE ) - { - printf( "\tSetting a polygon feature failed.\n" ); - exit( 1 ); - } + if ( poFeature->SetGeometry( &oPoly ) != OGRERR_NONE ) + { + printf( "\tSetting a polygon feature failed.\n" ); + exit( 1 ); + } - if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE ) - { - printf( "\tFailed to create feature in shapefile.\n" ); - exit( 1 ); - } + if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE ) + { + printf( "\tFailed to create feature in shapefile.\n" ); + exit( 1 ); + } - oRing.empty(); - oPoly.empty(); - OGRFeature::DestroyFeature( poFeature ); - } + oRing.empty(); + oPoly.empty(); + OGRFeature::DestroyFeature( poFeature ); + } //j + } //i } + else + { + // for arcGIS polygon shapefile text format + + bool newPoly = true; + string lineStr, dataStr; + i = 0.0; + vector vecString; + double x, y, x1, y1, x2, y2; + + + //get first polyID line + getline( imageStream, lineStr); + + while ( newPoly ) + { + lineStr = trim (lineStr); //get rid of spaces at edges + i++; //count the line number + + vecString = string2stringVector (lineStr, ","); + + if ( vecString.size() != 3 ) + { + printf( "\tError: new poly has ID, Longitude, Latitude - line %d %s\n",i, lineStr.c_str() ); + exit( 1 ); + } + + dataStr = vecString[0]; + gid = atoi ( dataStr.c_str() ); + + dataStr = vecString[1]; + x = atof ( dataStr.c_str() ); + + dataStr = vecString[2]; + y = atof ( dataStr.c_str() ); + + if ( x > 180 ) x = x - 360.0; + //printf ("%s\n", lineStr.c_str()); + printf ("%d,%.13lf,%.13lf,%d\n", gid, x, y,gid ); //print out lable point for arcGIS + //printf ("\tPoly: gid=%d line %d:%s\n", gid, i, lineStr.c_str()); + + + int cross180 = 0; //to check whether the polygon crosses 180 long. line + + //get all points + vector vecX, vecY; + int j1, j2; + j = 0; + getline( imageStream, lineStr); + + do { + lineStr = trim (lineStr); //get rid of spaces at edges + i++; //count the line number + j++; //total points in each polygon + + vecString = string2stringVector (lineStr, ","); + + if ( vecString.size() != 2 ) + { + printf( "\tError: Poly point has Longitude, Latitude - line %d %s\n",i, lineStr.c_str() ); + exit( 1 ); + } + + dataStr = vecString[0]; + x = atof ( dataStr.c_str() ); + + dataStr = vecString[1]; + y = atof ( dataStr.c_str() ); + + //point behind from the second point + j1 = j - 2; + //checking crossing 180 from the second point + if ( j > 1 && ( (vecX[j1] < 180.0 && vecX[j1] > 20.0 && x > 180.0 && x < 260.0) + || (vecX[j1] > 180.0 && vecX[j1] < 260.0 && x < 180.0 && x > 100.0) ) ) + { + //crossing 180 and compute crossing point + cross180 ++; + double m; + m = (y - vecY[j1]) / (x - vecX[j1]); + double b = y - m*x; + y1 = m * 180.0 + b; + + //add the 180 long. intersection point + vecX.push_back ( 180.0 ); + vecY.push_back ( y1 ); + j++; + } + + + vecX.push_back ( x ); + vecY.push_back ( y ); + + // printf ("\tPoint i=%d x=%lf y=%lf line: %s\n", i, x, y, lineStr.c_str() ); + + getline( imageStream, lineStr); + + } while ( lineStr.find ( "end" ) == string::npos ); + + + //check last point again the first point + x = vecX[0]; + y = vecY[0]; + j1 = vecX.size() - 1; + if ( ( (vecX[j1] < 180.0 && x > 180.0 ) || (vecX[j1] > 180.0 && x < 180.0 ) ) && cross180 == 1 ) + { + cross180 ++; + j++; + + double m; + m = (y - vecY[j1]) / (x - vecX[j1]); + double b = y - m*x; + y1 = m * 180.0 + b; + + //add the 180 long. intersection point + vecX.push_back ( 180.0 ); + vecY.push_back ( y1 ); + } + + + //printf ("\t\t: cross180=%d\n", cross180); + if ( cross180 > 2 ) + { + exit (1); + } + if ( cross180 ==2 ) + { + //cross180 = 0; + } + else if ( cross180 == 1 ) + { + printf ("\t\t: cross180=%d\n", cross180); + exit (1); + } + else if ( cross180 == 0 ) + { + cross180 = 1; //loop through one poly + } + + + if ( j != vecX.size() ) + { + printf ("\t\t: j=%d not equal vecX size=%d\n", j, vecX.size()); + exit (1); + } + + + int k=0; //beginning index in m polygon + int k2 = vecX.size() - 1; // end index in m polygon + + int begin_j = 0; + int end_j = 0; + int jj = 0; + double x0, y0, x3, y3; //keep beginning point + + //write out 1 or 2 polys + for (m=0; mGetLayerDefn() ); + poFeature->SetField( "GRIDID", gid ); + + for (j=k; j<=k2; j++) + { + x = vecX[j]; + y = vecY[j]; + + if ( cross180 == 2 ) + { + //printf ("\tPoly: gid=%d m=%d pt=%d: x=%lf y=%lf\n", gid,m,j,x,y); + } + + + if (x == 180.0 && m == 0 ) + { + //point behind and front + j1 = j - 1; + j2 = j + 1; + if ( j == 0 ) j1 = vecX.size() - 1; + if ( j == (vecX.size() - 1) ) j2 = 0; + + x1 = vecX[j1]; + x2 = vecX[j2]; + + if ( x1 > 180.0 ) x = -180.0; //on the left + if ( x1 < 180.0 ) x = 180.0; //on the right + + if ( cross180 == 2 ) + { + //printf ("\t\tIn1: gid=%d m=%d pt=%d: x=%lf y=%lf\n", gid,m,j,x,y); + } + + oRing.addPoint( x, y ); + + //crossing point + if ( (x1 > 180.0 && x2 < 180.0) || (x1 < 180.0 && x2 > 180.0) ) + { + begin_j = j; //beginning point for next splitted polygon + + //find next 180 point in the same 180 side + for (jj=j+1; jj 180 ) x = x - 360.0; //get nagative -180 to 180 for longitude + + if ( cross180 == 2 ) + { + //printf ("\t\tIn4: gid=%d m=%d pt=%d: x=%lf y=%lf\n", gid,m,j,x,y); + } + + oRing.addPoint( x, y); + } + + //keep first point + if (j == k ) + { + x0 = x; + y0 = y; + } + + } //j + + //close poly + oRing.addPoint( x0, y0 ); + + oPoly.addRing( &oRing ); + + if ( poFeature->SetGeometry( &oPoly ) != OGRERR_NONE ) + { + printf( "\tSetting a polygon feature failed.\n" ); + exit( 1 ); + } + + if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE ) + { + printf( "\tFailed to create feature in shapefile.\n" ); + exit( 1 ); + } + + oRing.empty(); + oPoly.empty(); + OGRFeature::DestroyFeature( poFeature ); + + //set index range for next poly = two splited ploys + k = begin_j; + k2 = end_j; + + } // end of m + + + getline( imageStream, lineStr); + + if ( lineStr.find ( "end" ) != string::npos ) + { + newPoly = false; //double end to end of the reading + } + } // while new poly + } //reading poly text file + GDALClose( poDS ); - printf ("\tCompleted in creating the grid shapefile.\n\n"); + + printf ("\tCompleted in creating the grid shapefile: %s.\n\n", shapeName.c_str() ); return shapeName; @@ -286,6 +604,8 @@ gridInfo getImageInfo (string imageFile) { // GDAL image, not HDF4 image isHDF4 = false; + + printf ("\ttest1: not hdf4\n"); } else if ( strValues.size() == 2 ) { @@ -306,9 +626,17 @@ gridInfo getImageInfo (string imageFile) /* -------------------------------------------------------------------- */ /* Register format(s). */ /* -------------------------------------------------------------------- */ + + printf ("\ttest2: open\n"); + GDALAllRegister(); + printf ("\ttest3: open\n"); + poRDataset = (GDALDataset *) GDALOpen( imageFile.c_str(), GA_ReadOnly ); + + printf ("\ttest4 end of open\n"); + if( poRDataset == NULL ) { printf( "\tOpen GDAL image file failed: %s.\n", imageFile.c_str() ); @@ -1155,9 +1483,14 @@ void computeGridCornerCenterLatLong ( gridInfo *grid ) exit ( 1 ); } - //get the maximum x and y for the domain - grid->xmax = grid->xmin + grid->xCellSize*grid->cols; - grid->ymax = grid->ymin + grid->yCellSize*grid->rows; + //get the maximum x and y for the domain + + // for grid domain + if ( grid->name.find( ".shp" ) == string::npos ) + { + grid->xmax = grid->xmin + grid->xCellSize*grid->cols; + grid->ymax = grid->ymin + grid->yCellSize*grid->rows; + } if ( pj_is_latlong(proj4DF) ) { @@ -1556,7 +1889,8 @@ void writeWRFGlobalAtributes ( int ncid, string title, gridInfo grid, string sta //for WRF simulation - 40 NLCD data - string gmtS = startDateTime.substr(8, 4); + //string gmtS = startDateTime.substr(8, 4); + string gmtS = string ("0.0"); string yearS = startDateTime.substr(0, 4); tmp_str = startDateTime.substr(0, 8); @@ -2416,7 +2750,6 @@ gridInfo getHDF4VarInfo ( string imageFile, string varName ) //if ( imageFile.find( "MCD43A3" ) != string::npos ) exit ( 1 ); - //data type 6 sprintf (buffer, "%d\0",num_type); imageInfo.attsStr.push_back ( string ( buffer ) ); @@ -2425,17 +2758,19 @@ gridInfo getHDF4VarInfo ( string imageFile, string varName ) //add_offset 7 valueStr = getHDF4VarAttrInfo (sds_id, "long_name"); + if ( valueStr.size() != 0) { imageInfo.attsStr.push_back ( valueStr ); //printf("\tlong_name=%s\n", valueStr.c_str() ); } + anyHDF4Errors ( SDendaccess(sds_id) ); if ( imageFile.find( "MCD12Q1" ) != string::npos || imageFile.find( "MOD12Q1" ) != string::npos || + imageFile.find( "MCD15A3H" ) != string::npos || imageFile.find( "MCD15A2H" ) != string::npos || imageFile.find( "MOD15A2GFS" ) != string::npos || imageFile.find( "MOD15A2" ) != string::npos || - imageFile.find( "MCD15A2H" ) != string::npos || imageFile.find( "MCD43A3" ) != string::npos || imageFile.find( "MCD43A1" ) != string::npos || imageFile.find( "MCD43A2" ) != string::npos ) { @@ -2468,10 +2803,15 @@ gridInfo getHDF4VarInfo ( string imageFile, string varName ) { imageInfo.name = string ( "MOD15A2" ); } + else if ( imageFile.find( "MCD15A3H" ) != string::npos ) + { + imageInfo.name = string ( "MCD15A3H" ); + } else if ( imageFile.find( "MCD15A2H" ) != string::npos ) { imageInfo.name = string ( "MCD15A2H" ); } + else if ( imageFile.find( "MCD43A3" ) != string::npos ) { imageInfo.name = string ( "MCD43A3" ); @@ -2488,7 +2828,7 @@ gridInfo getHDF4VarInfo ( string imageFile, string varName ) anyHDF4Errors ( SDattrinfo(sd_id, i, attr_name, &data_type, &count) ); - //printf ( "\tAttribute Index: %d Attribute name: %s Attribute Type: %d Attribute count: %d\n", i,attr_name,data_type, count); + printf ( "\tAttribute Index: %d Attribute name: %s Attribute Type: %d Attribute count: %d\n", i,attr_name,data_type, count); if ( ( attBuff = (char *) calloc (count, sizeof(char)) ) == NULL) @@ -2497,11 +2837,12 @@ gridInfo getHDF4VarInfo ( string imageFile, string varName ) exit ( 1 ); } + anyHDF4Errors ( SDreadattr(sd_id, i, attBuff) ); string attStr = string ( attBuff ); //printf ("\t\tAttribute: %s:\n",attBuff); - + metaStr = string2stringVector (attStr, "\n"); @@ -2586,10 +2927,10 @@ gridInfo getHDF4VarInfo ( string imageFile, string varName ) } } //end of lines + if ( ( pos = attStr.find ( "CHARACTERISTICBINSIZE" ) ) != string::npos ) { - - if ( attStr.find ( "CHARACTERISTICBINSIZE500M" ) != string::npos ) + if ( attStr.find ( "CHARACTERISTICBINSIZE500M" ) != string::npos ) { valueStr = attStr.substr( pos+82, 11 ); } @@ -2688,8 +3029,11 @@ string getHDF4VarAttrInfo (int32 sds_id, const char *attName ) } anyHDF4Errors ( SDreadattr(sds_id, attr_index, attBuff) ); + attStr = string ( attBuff ); - attStr.erase (count); + + //attStr.erase (count); + } else if ( data_type == DFNT_FLOAT64 ) //data type = 6 { @@ -2803,21 +3147,32 @@ gridInfo getHDF5VarInfo ( string imageFile, string varName ) int numDims, numMems; hsize_t dims[6]; //set max dimensions to 6 string valueStr; + string fileType; int j; //open the file + anyHDF5Errors ( file_id = H5Fopen( imageFile.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT ) ); + if ( imageFile.find ( ".he5" ) != string::npos ) { groupPath = getHDF5VarGroups ( file_id, varName ); + fileType = string ( "omi"); + } + else if ( imageFile.find ( "S5P_OFFL_L2" ) != string::npos ) + { + groupPath =string ("/PRODUCT/"); + fileType = string ( "tropomi"); } else { groupPath =string ("/"); //no path + fileType = string ( "other"); } + dataSet = groupPath + varName; @@ -2836,9 +3191,14 @@ gridInfo getHDF5VarInfo ( string imageFile, string varName ) for (j=0; j= 2 ) { imageInfo.rows = dims[0]; @@ -2863,6 +3223,17 @@ gridInfo getHDF5VarInfo ( string imageFile, string varName ) imageInfo.xCellSize = 0.05; imageInfo.yCellSize = 0.05; } + + //this is TROPOMI L2 data + if (fileType.compare ("tropomi") == 0 ) + { + imageInfo.rows = dims[1]; + imageInfo.cols = dims[2]; + + // tropomi resolution is iaround 3.5x7 km, 3/111.1 degree resolution + imageInfo.xCellSize = 3.5/111.1; + imageInfo.yCellSize = 3.5/111.1; + } //Processing data set attributes @@ -2871,15 +3242,36 @@ gridInfo getHDF5VarInfo ( string imageFile, string varName ) imageInfo.attsStr.push_back ( varName ); //Units 2 - valueStr = getHDF5VarAttrInfo (dataset_id, "Units"); + if (fileType.compare ("tropomi") == 0 ) + { + valueStr = getHDF5VarAttrInfo (dataset_id, "units"); + } + else + { + valueStr = getHDF5VarAttrInfo (dataset_id, "Units"); + } imageInfo.attsStr.push_back ( valueStr ); //scale factor 3 - valueStr = getHDF5VarAttrInfo (dataset_id, "ScaleFactor"); + if (fileType.compare ("tropomi") == 0 ) + { + valueStr = string ("1.0"); + } + else + { + valueStr = getHDF5VarAttrInfo (dataset_id, "ScaleFactor"); + } imageInfo.attsStr.push_back ( valueStr ); //add_offset 4 - valueStr = getHDF5VarAttrInfo (dataset_id, "Offset"); + if (fileType.compare ("tropomi") == 0 ) + { + valueStr = string ("0.0"); + } + else + { + valueStr = getHDF5VarAttrInfo (dataset_id, "Offset"); + } imageInfo.attsStr.push_back ( valueStr ); //_FillValue 5 @@ -2910,7 +3302,14 @@ gridInfo getHDF5VarInfo ( string imageFile, string varName ) //add LONGNAME 7 - valueStr = getHDF5VarAttrInfo (dataset_id, "Title"); + if (fileType.compare ("tropomi") == 0 ) + { + valueStr = getHDF5VarAttrInfo (dataset_id, "long_name"); + } + else + { + valueStr = getHDF5VarAttrInfo (dataset_id, "Title"); + } imageInfo.attsStr.push_back ( valueStr ); @@ -3037,14 +3436,22 @@ void readHDF5SatVarData (string imageFile, string varName, double *poImage) hsize_t dims[6]; //set max dimensions to 6 string valueStr; + string fileType; + //open the file anyHDF5Errors ( file_id = H5Fopen( imageFile.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT ) ); /* obtain group path under "/HDFEOS/SWATHS" */ - if ( imageFile.find ( ".he5" ) != string::npos ) + if ( imageFile.find ( ".he5" ) != string::npos) { groupPath = getHDF5VarGroups ( file_id, varName ); + fileType = string ( "omi"); + } + else if ( imageFile.find ( "S5P_OFFL_L2" ) != string::npos ) + { + groupPath =string ("/PRODUCT/"); + fileType = string ( "tropomi"); } else { @@ -3078,6 +3485,13 @@ void readHDF5SatVarData (string imageFile, string varName, double *poImage) int rows = dims[0]; int cols = dims[1]; + if ( imageFile.find ( "S5P_OFFL_L2" ) != string::npos ) + { + //1d is 1 for tropomi + rows = dims[1]; + cols = dims[2]; + } + //Define the memory space to read dataset. mspace_id = H5Screate_simple(numDims,dims,NULL); @@ -3128,6 +3542,7 @@ string getHDF5VarGroups (hid_t file_id, string varName) size_t size; + string startPath = string ( "/HDFEOS/SWATHS" ); printf("\tGet variable HDF5 group path: %s\n", varName.c_str() ); @@ -3377,12 +3792,22 @@ int defineWRFNetCDFSatVar ( int ncid, string satVarName, int numDims, int *dimIn anyErrors( nc_put_att_text(ncid, satV_id, "description", imageInfo.attsStr[6].size(), imageInfo.attsStr[6].c_str() )); - anyErrors( nc_put_att_text(ncid, satV_id, "units", imageInfo.attsStr[1].size(), imageInfo.attsStr[1].c_str() ) ); + + //tropomi hour variable + if ( satVarName.find ("scanline_average_hours") != string::npos ) + { + string unitInfo = string ("hours of the day"); + anyErrors( nc_put_att_text(ncid, satV_id, "units", unitInfo.size(), unitInfo.c_str() ) ); + } + else + { + anyErrors( nc_put_att_text(ncid, satV_id, "units", imageInfo.attsStr[1].size(), imageInfo.attsStr[1].c_str() ) ); + } if ( satVarName.find ("_LAI") != string::npos || satVarName.find ("_FPAR") != string::npos || satVarName.find ("Albedo") != string::npos || - satVarName.find ("ALBEDO") != string::npos || satVarName.find ("Lai_1km") != string::npos || satVarName.find ("MODIS_ALB_") != string::npos || satVarName.find ("Lai_500m") != string::npos || satVarName.find ("Fpar_500m") != string::npos ) + satVarName.find ("ALBEDO") != string::npos || satVarName.find ("Lai_1km") != string::npos || satVarName.find ("MODIS_ALB_") != string::npos || satVarName.find ("Lai_500m") != string::npos || satVarName.find ("Fpar_500m") != string::npos ) { attValue[0] = 1.0; anyErrors( nc_put_att_float(ncid, satV_id, "scale_factor", NC_FLOAT, 1, attValue) ); @@ -3442,6 +3867,7 @@ void defineWRFNCTimeVars ( int ncid, int time_dim, int dateStr_dim, string start { int dimIndex[NC_MAX_DIMS]; //dimension index array for write out array + const int timeStrLength = 19; //time string length in WRF Netcdf output printf( "Define time variables in output netcdf file...\n" ); @@ -3461,7 +3887,8 @@ void defineWRFNCTimeVars ( int ncid, int time_dim, int dateStr_dim, string start dimIndex[1] = dateStr_dim; anyErrors( nc_def_var(ncid, "Times", NC_CHAR, 2, dimIndex, timeStr_id) ); anyErrors( nc_put_att_text(ncid, *timeStr_id, "description", 21, "Time string (19 char)" ) ); - anyErrors( nc_put_att_text(ncid, *timeStr_id, "units", 15, "yyyymmddhhmm:ss") ); + //anyErrors( nc_put_att_text(ncid, *timeStr_id, "units", 15, "yyyymmddhhmm:ss") ); + anyErrors( nc_put_att_text(ncid, *timeStr_id, "units", timeStrLength, "yyyy-mm-dd_hh:mm:ss") ); } @@ -3683,6 +4110,10 @@ void computeGridSatValues ( GUInt32 *poImage_grd, int *grdIndex, float *satV, d int i,j,k,m; + //for tropomi + int *poImage_qa = NULL; + int *deltaTime_2010010 = NULL; //millisecond + /******************************************* * get info from the model domain grid * @@ -3694,11 +4125,12 @@ void computeGridSatValues ( GUInt32 *poImage_grd, int *grdIndex, float *satV, d int totalDomainSize = gridCols*gridRows; //for output dimensions + for ( i=2; i modisFiles, string varName if ( i == 0 ) { + + allInfo = getHDF4VarInfo (imageFile, varName); - //printGridInfo ( allInfo ); } else { gridInfo modisInfo = getHDF4VarInfo (imageFile, varName); - //printGridInfo ( modisInfo ); allInfo.xmin = min ( allInfo.xmin, modisInfo.xmin ); @@ -5931,6 +6398,8 @@ string extractDomainLAIData ( std::vector modisFiles, std::vector modisFiles, std::vector modisFiles, std::vectorname; + printf( "\nShapefile extent: get shapefile extent - %s.\n", shapeFile.c_str() ); +/* -------------------------------------------------------------------- */ +/* Open shapefile */ +/* -------------------------------------------------------------------- */ + poDS = (GDALDataset*) GDALOpenEx( shapeFile.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL ); + if( poDS == NULL ) + { + printf( "\tError: Opening Shapefile file failed - %s.\n", shapeFile.c_str() ); + exit( 1 ); + } + +/* -------------------------------------------------------------------- */ +/* Get layer extent to make new raster domain. */ +/* -------------------------------------------------------------------- */ + poLayer = poDS->GetLayer( 0 ); + if( poLayer == NULL ) + { + printf( "\tError: Couldn't fetch layer %s!\n", shapeFile.c_str() ); + exit( 1 ); + } + + if (poLayer->GetExtent(&oExt, TRUE) == OGRERR_FAILURE) + { + printf("\tNo extent from the shapefile\n", shapeFile.c_str() ); + exit ( 1 ); + } + + printf("\tShapefile extent min_xy: (%lf, %lf) max_xy: (%lf, %lf)\n", + oExt.MinX, oExt.MinY, oExt.MaxX, oExt.MaxY); + + //compute new extent to match resolution + xMin = oExt.MinX; + yMin = oExt.MinY; + + xMax = oExt.MaxX; + yMax = oExt.MaxY; + + GDALClose ( poDS ); + + //assign info to grid + grid->xmin = xMin; + grid->ymin = yMin; + + grid->xmax = xMax; + grid->ymax = yMax; + + grid->xCellSize = 0.0; + grid->yCellSize = 0.0; +} + +/********************************************* +* 67. readGridCoordinates * +*********************************************/ +std::map readGridCoordinates ( gridInfo *grid, string polyCenterFile) +{ + + ifstream imageStream; //stream to read input EPIC site file + double x,y; + int row,col; //starting from 0 + string lineStr,dataStr; //for reading in site line + vector vecString; + + vector centerX, centerY; //center point long. and lat. + + int i,j; + std::map mpasIDs; //gridID, mpasID + std::map gridIDs; //reading in sequence, gridID + int mpasID, gridID, gridN; + + + printf( "\nReading polygon center point locations from file - %s\n",polyCenterFile.c_str() ); + + i = 0; //lines read + + imageStream.open( polyCenterFile.c_str() ); + if (imageStream.good() ) + { + getline( imageStream, lineStr, '\n'); + while ( !imageStream.eof() ) + { + lineStr = trim (lineStr); //get rid of spaces at edges + i++; //count the line number + + //put data in string vector + vecString = string2stringVector (lineStr, ","); + if ( vecString.size() < 4 ) //mpasID, long, lat, gridID + { + printf( "\tError: should contain at least 4 items: MPAS_ID,Longitude,Latitude,GRIDID. %d - %s\n", i,lineStr.c_str() ); + exit ( 1 ); + } + + if ( i == 1 ) + { + printf( "\tFile header: %s\n", lineStr.c_str() ); + } + else + { + //get mpas ID + dataStr = vecString[0]; + mpasID = atoi ( dataStr.c_str() ); + + //get longitude + dataStr = vecString[1]; + x = atof ( dataStr.c_str() ); + + //get latitude + dataStr = vecString[2]; + y = atof ( dataStr.c_str() ); + + //get grid ID + dataStr = vecString[3]; + gridID = atoi ( dataStr.c_str() ); + + printf ( "\tLine: %s\n", lineStr.c_str() ); + printf ( "\tLine: %d: %d %lf %lf %d\n", i,mpasID,x,y,gridID ); + + + //store site info in vectors + mpasIDs[gridID]=mpasID; //gridID from 1 + gridIDs[i-2]= gridID; //keep a record on the sequence with gridID from 0 + centerX.push_back ( x ); + centerY.push_back ( y ); + } //site lines + + getline( imageStream, lineStr); + } // while loop + + printf ("\tFinished reading: %d lines\n\n", i); + imageStream.close(); + + } // if good + else + { + printf ("\tError: Open file - %s\n", polyCenterFile.c_str() ); + exit ( 1 ); + } + + + gridN = centerX.size(); + printf ("\tNumber of polygons: %d\n\n", gridN); + + //set the size for 1-D netCDF output + grid->cols = gridN; + grid->rows = 1; + + + //allocate memory + //x and y array 1d + if ( (grid->x = (float*) calloc (grid->cols, sizeof(float)) ) == NULL) + { + printf( "Calloc x failed.\n"); + exit ( 1 ); + } + if ( (grid->y = (float*) calloc (grid->rows, sizeof(float)) ) == NULL) + { + printf( "Calloc y failed.\n"); + exit ( 1 ); + } + + //lat and long array 2d + if ( (grid->lat = (float*) calloc (grid->rows*grid->cols, sizeof(float)) ) == NULL) + { + printf( "Calloc lat failed.\n"); + exit ( 1 ); + } + if ( (grid->lon = (float*) calloc (grid->rows*grid->cols, sizeof(float)) ) == NULL) + { + printf( "Calloc lon failed.\n"); + exit ( 1 ); + } + + //x and y array for grid is dummy values for polygon shapefile output + //mass x + for (i=0; icols; i++) + { + gridID = gridIDs[i]; + grid->x[gridID] = centerX[i]; + } + + //mass y + for (j=0; jrows; j++) + { + gridID = gridIDs[j]; + grid->y[gridID] = centerY[j]; + } + + for(j=0; jrows; j++) + { + for (i=0; icols; i++) + { + int index = j * grid->cols + i; + gridID = gridIDs[index]; + + grid->lon[gridID] = centerX[index]; + grid->lat[gridID] = centerY[index]; + } + } + + centerX.clear(); + centerY.clear(); + gridIDs.clear(); + + return mpasIDs; +} + + +/********************************************* +* 68. readHDF5SatVarDataInt * +* read a HDF5 file variable data * +*********************************************/ +void readHDF5SatVarDataInt (string imageFile, string varName, int *poImage) +{ + + string groupPath, dataSet; + int totalSize = 1; + int i,j; + + /********************************* + * read HDF5 variable info * + *********************************/ + hid_t file_id, dataset_id; + hid_t fspace_id, mspace_id; + hid_t dtype_id, ndtype_id; /* identifiers */ + H5T_class_t dclass_id; + int numDims, numMems; + hsize_t dims[6]; //set max dimensions to 6 + string valueStr; + + string fileType; + + + //open the file + anyHDF5Errors ( file_id = H5Fopen( imageFile.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT ) ); + + /* obtain group path under "/HDFEOS/SWATHS" */ + if ( imageFile.find ( ".he5" ) != string::npos) + { + groupPath = getHDF5VarGroups ( file_id, varName ); + fileType = string ( "omi"); + } + else if ( imageFile.find ( "S5P_OFFL_L2" ) != string::npos ) + { + groupPath =string ("/PRODUCT/"); + fileType = string ( "tropomi"); + } + else + { + groupPath =string ("/"); //no path for hdf5 file like OMI NO2 L3 0.05 grids + } + + dataSet = groupPath + varName; + + printf( "\tRead data set: %s\n", dataSet.c_str() ); + + //open data set + anyHDF5Errors ( dataset_id = H5Dopen (file_id, dataSet.c_str(), H5P_DEFAULT ) ); + + //get dimension size of HDF5 variable + anyHDF5Errors ( fspace_id = H5Dget_space(dataset_id) ); /* Get filespace handle first. */ + anyHDF5Errors ( numDims = H5Sget_simple_extent_dims(fspace_id, dims, NULL) ); + printf ( "\tNumber of Dimensions =%d \n", numDims ); + if ( numDims < 2 ) + { + printf ( "\tError: Variable dimension should be >= 2. numDims = %d\n", numDims ); + exit ( 1 ); + } + + totalSize = 1; + for (i=0; i 2 ) + { + //1d is 1 for tropomi + rows = dims[1]; + cols = dims[2]; + } + + //Define the memory space to read dataset. + mspace_id = H5Screate_simple(numDims,dims,NULL); + + /* + * Read dataset back and display. + */ + anyHDF5Errors ( H5Dread(dataset_id, H5T_NATIVE_INT, mspace_id, fspace_id, H5P_DEFAULT, poImage) ); + + /*checking + double tempD[rows][cols]; + anyHDF5Errors ( H5Dread(dataset_id, H5T_NATIVE_DOUBLE, mspace_id, fspace_id, H5P_DEFAULT, tempD) ); + printf("Dataset: \n"); + for (j = 0; j < dims[0]; j++) + { + for (i = 0; i < dims[1]; i++) + { + int index = j*cols + i; + printf("%.4lf %.4lf", poImage[index], tempD[j][i]); + if ( poImage[index] != tempD[j][i] ) + { + printf("\t\tDo not match\n"); + exit ( 1 ); + } + printf("\n"); + } + } + */ + + anyHDF5Errors ( H5Sclose ( mspace_id ) ); + anyHDF5Errors ( H5Sclose ( fspace_id ) ); + anyHDF5Errors ( H5Dclose ( dataset_id ) ); + anyHDF5Errors ( H5Fclose ( file_id ) ); + + //printf ( "\tRead data: %s | %s\n", imageFile.c_str(), varName.c_str() ); + +} + + + +/**********************************************************/ +/* 69. compute domain grid satellite image value */ +/**********************************************************/ +void computeGridSatValues_hour ( GUInt32 *poImage_grd, int *grdIndex, float *satV, float *hourV, double *poImage, + int *poImage_qa, int *poImage_time, gridInfo imageInfo, gridInfo newRasterInfo, + gridInfo grid, int qa_value_Min ) +{ + double halfGridCellArea; //half of domain cell area + int *gridIDs=NULL; //array to store number of pixels in each modeling grids + double *gridValues = NULL; //array to store total values in each modeling grids + double *gridTimes = NULL; //array to store total Times in each modeling cells + + int i,j,k,m; + + + /******************************************* + * get info from the model domain grid * + *******************************************/ + //total dimension size + + int gridRows = grid.rows; + int gridCols = grid.cols; + + int totalDomainSize = gridCols*gridRows; //for output dimensions + + + for ( i=2; i 0 && gridID <= gridPixels) + { + //get row and col indices in the model domain from LL corner + int dmnRow = (int) ( floor ( ( gridID - 1 ) / gridCols) ); //start from 0 + int dmnCol = ( gridID - 1 ) % gridCols; //start from 0 + //printf ("\t\tgridID=%d dmnRow=%d dmnCol=%d\n",gridID,dmnRow,dmnCol); + + geoIndex = grdIndex[pixelIndex]; //get the index to get the value image rowsXcols starting from 0 + + if ( geoIndex != -999 ) + { + //get sat image geo indices for (1st and 2nd dimensions) based on grdIndex result + //compute geo-site row, col + int satRow = (int) ( floor (geoIndex / xCells) ); //start from 0 + int satCol = geoIndex % xCells; //start from 0 + //printf ("\tgeoIndex=%d row=%d col=%d\n",geoIndex,satRow,satCol); + + if ( imageInfo.dims.size() == 2 ) + { + satIndex = geoIndex; + dmnIndex = gridID - 1; + + value = poImage[satIndex]; + qaValue = poImage_qa[satIndex]; + timeD = poImage_time[satRow]; + timeD = timeD / 3600000.0; //milliseconds to hours + + if ( value != SAT_MISSIING_VALUE && qaValue >= qa_value_Min ) + { + gridIDs[dmnIndex] += 1; //count Satellite cells + gridValues[dmnIndex] += value; //sum the Satellite value for a domain grid + gridTimes[dmnIndex] += timeD; + } + } //2dims + else if ( imageInfo.dims.size() == 3 ) + { + for ( k=0; k= halfGridCellArea ) + { + value = gridValues[dmnIndex] / gridIDs[dmnIndex]; //take average for all satellite variables + satV[dmnIndex] = value; + if ( hourV != NULL ) hourV[dmnIndex] = gridTimes[dmnIndex] / gridIDs[dmnIndex]; + } + else + { + satV[dmnIndex] = MISSIING_VALUE; + if ( hourV != NULL ) hourV[dmnIndex] = MISSIING_VALUE; + } + } //2dims + + else if ( imageInfo.dims.size() == 3 ) + { + for ( k=0; k= halfGridCellArea ) + { + value = gridValues[dmnIndex] / gridIDs[dmnIndex]; //take average for all satellite variables + satV[dmnIndex] = value; + if ( hourV != NULL ) hourV[dmnIndex] = gridTimes[dmnIndex] / gridIDs[dmnIndex]; + } + else + { + satV[dmnIndex] = MISSIING_VALUE; + if ( hourV != NULL ) hourV[dmnIndex] = MISSIING_VALUE; + } + } //k + } //3dims + + else if ( imageInfo.dims.size() == 4 ) + { + for ( k=0; k= halfGridCellArea ) + { + value = gridValues[dmnIndex] / gridIDs[dmnIndex]; //take average for all satellite variables + satV[dmnIndex] = value; + if ( hourV != NULL ) hourV[dmnIndex] = gridTimes[dmnIndex] / gridIDs[dmnIndex]; + } + else + { + satV[dmnIndex] = MISSIING_VALUE; + if ( hourV != NULL ) hourV[dmnIndex] = MISSIING_VALUE; + } + } //m + } //k + } //4 dims + + } //j + } //i + + CPLFree (gridIDs); + CPLFree (gridValues); + CPLFree (gridTimes); +} diff --git a/src/raster/geotools.h b/src/raster/geotools.h index dd9f149..47f0034 100755 --- a/src/raster/geotools.h +++ b/src/raster/geotools.h @@ -141,3 +141,8 @@ void getPointsInRasterValues (gridInfo grid, std::vector epicSiteIDs, std: int getNetCDFDim ( char * dimName, string fileName ); string extractDomainLAIData ( std::vector modisFiles, std::vector satVars, gridInfo grid); string extractDomainALBData ( std::vector modisFiles, std::vector satVars, gridInfo grid); +void getShapeFile_Extent ( gridInfo *grid ); +std::map readGridCoordinates ( gridInfo *grid, string polyCenterFile); +void readHDF5SatVarDataInt (string imageFile, string varName, int *poImage); +void computeGridSatValues_hour ( GUInt32 *poImage_grd, int *grdIndex, float *satV, float *hourV, double *poImage, int *poImage_qa, int *poImage_time, gridInfo imageInfo, gridInfo newRasterInfo, gridInfo grid, int qa_value_Min ); + diff --git a/src/raster/preProcessNLCD.exe b/src/raster/preProcessNLCD.exe index 1f6b39b..1df0eb1 100755 Binary files a/src/raster/preProcessNLCD.exe and b/src/raster/preProcessNLCD.exe differ diff --git a/src/raster/rasterWtoPolygons.exe b/src/raster/rasterWtoPolygons.exe index 2b4ce5e..224d11a 100755 Binary files a/src/raster/rasterWtoPolygons.exe and b/src/raster/rasterWtoPolygons.exe differ diff --git a/src/raster/toDataAssimilationFMT.exe b/src/raster/toDataAssimilationFMT.exe index 3162e22..61ad1f7 100755 Binary files a/src/raster/toDataAssimilationFMT.exe and b/src/raster/toDataAssimilationFMT.exe differ diff --git a/src/raster/toNLCDRaster.exe b/src/raster/toNLCDRaster.exe index 86ea074..322df01 100755 Binary files a/src/raster/toNLCDRaster.exe and b/src/raster/toNLCDRaster.exe differ diff --git a/src/raster/txt2ncf.exe b/src/raster/txt2ncf.exe index 63dd35e..7c12a18 100755 Binary files a/src/raster/txt2ncf.exe and b/src/raster/txt2ncf.exe differ diff --git a/src/raster/utilities.cpp b/src/raster/utilities.cpp index 2a31768..88747f9 100755 --- a/src/raster/utilities.cpp +++ b/src/raster/utilities.cpp @@ -51,6 +51,7 @@ * 49. stringVector2string - convert string vector to string separated by "," * 50. fillFloatArrayMissingValueVar - fill missing value variable * 51. dayofweek - find the day of week (1-monday ... 7-sunday) + * 52. fillTimeArrays_Day - fill day time array * * Written by the Institute for the Environment at UNC, Chapel Hill * in support of the EPA NOAA CMAS Modeling, 2007-2008. @@ -1141,7 +1142,7 @@ vector obtainSatFileNames (string dataDir, string startDateTime, string timeEndStr = endDateTime.substr(8, 4); dayEndStr = dayEndStr_julian + timeEndStr; } - else if ( satType.compare ("OMI") == 0 ) + else if ( satType.compare ("OMI") == 0 || satType.compare ("TROPOMI") == 0 ) { dayStartStr = startDateTime; dayEndStr = endDateTime; @@ -1163,7 +1164,8 @@ vector obtainSatFileNames (string dataDir, string startDateTime, string //MODIS LAIFPAR: MOD15A2GFS 1km: MOD15A2GFS.A2006001.h07v03.005.2009023132100.hdf //MODIS LAIFPAR: MOD15A2 1km: MOD15A2.A2006057.h10v06.005.2008077144707.hdf //MODIS LAIFPAR: MCD15A2 1km: MCD15A2.A2011089.h09v03.005.2011103004948.hdf - //MODIS LAIFPAR: MCD15A2H 500m 8-day: MCD15A2H.A2016241.h15v03.006.2016250072611.hdf + //MODIS LAIFPAR: MCD15A2H 500m: MCD15A2H.A2016241.h15v03.006.2016250072611.hdf + //MODIS LAIFPAR: MCD15A3H 500m: MCD15A3H.A2016073.h11v04.006.2016110203715.hdf //MODIS albedo: MCD43A3 500km: MCD43A3.A2006177.h14v03.005.2008134144426.hdf //MODIS albedo: MCD43A1 500km: MCD43A1.A2006121.h08v04.005.2008117094839.hdf //MODIS albedo: MCD43A2 500km: MCD43A2.A2006121.h08v04.005.2008117094839.hdf @@ -1181,16 +1183,18 @@ vector obtainSatFileNames (string dataDir, string startDateTime, string tmp_str = string(dirpFiles->d_name); tmp_str = trim(tmp_str); + if ( ( ( tmp_str.find( "D06_L2.A" ) != string::npos || tmp_str.find( "D04_L2.A" ) != string::npos || tmp_str.find( "MISR_AS_AEROSOL_F12_0022" ) != string::npos || tmp_str.find( "MCD12Q1.A" ) != string::npos || tmp_str.find( "MOD12Q1.A" ) != string::npos || tmp_str.find( "MOD15A2GFS.A" ) != string::npos || tmp_str.find( "MOD15A2.A" ) != string::npos || - tmp_str.find( "MCD15A2.A" ) != string::npos || tmp_str.find( "MCD15A2H.A" ) != string::npos || + tmp_str.find( "MCD15A2.A" ) != string::npos || tmp_str.find( "MCD15A3H.A" ) != string::npos || + tmp_str.find( "MCD15A2H.A" ) != string::npos || tmp_str.find( "MCD43A3.A" ) != string::npos || tmp_str.find( "MCD43A1.A" ) != string::npos || tmp_str.find( "MCD43A2.A" ) != string::npos ) && tmp_str.find( ".hdf" ) != string::npos ) || ( tmp_str.find( "OMI-Aura_L2-OM" ) != string::npos && tmp_str.find( ".he5" ) != string::npos && - tmp_str.find( ".he5.xml" ) == string::npos) ) + tmp_str.find( ".he5.xml" ) == string::npos) || (tmp_str.find( "S5P_OFFL_L2" ) != string::npos) ) { imageDate = getDayTimeStrFromSatFileName ( tmp_str, satType ); dayTimeMid = atol ( imageDate.c_str() ); @@ -1306,6 +1310,9 @@ string getDayTimeStrFromSatFileName ( string imageFileName, string satType ) string dayStr, timeStr; + char numChars[10]; + + i = imageFileName.rfind("/", imageFileName.length() ); if (i != string::npos) { @@ -1323,14 +1330,14 @@ string getDayTimeStrFromSatFileName ( string imageFileName, string satType ) } else if ( imageFileName.find( "MCD12Q1.A" ) != string::npos || imageFileName.find( "MOD12Q1.A" ) != string::npos || imageFileName.find( "MOD15A2.A" ) != string::npos || - imageFileName.find( "MCD15A2.A" ) != string::npos || + imageFileName.find( "MCD15A2.A" ) != string::npos || imageFileName.find( "MCD43A3.A" ) != string::npos || imageFileName.find( "MCD43A1.A" ) != string::npos || imageFileName.find( "MCD43A2.A" ) != string::npos ) //MODIS land cover products 500m and 1km { dayStr = imageFileName.substr(9, 7); timeStr = string ( "0000" ); //yearly products } - else if ( imageFileName.find( "MCD15A2H.A" ) != string::npos ) + else if ( imageFileName.find( "MCD15A3H.A" ) != string::npos || imageFileName.find( "MCD15A2H.A" ) != string::npos ) { dayStr = imageFileName.substr(10, 7); timeStr = string ( "0000" ); //yearly products @@ -1362,6 +1369,30 @@ string getDayTimeStrFromSatFileName ( string imageFileName, string satType ) dayStr.erase (4,1); //get rid of m in time: 2006m0422 timeStr = imageFileName.substr(28, 4); } + else if ( imageFileName.find( "S5P_OFFL_L2" ) != string::npos ) + { + dayStr = imageFileName.substr(20, 8); + + //get mid of HH MM + string tempH1 = imageFileName.substr(29, 2); + string tempM1 = imageFileName.substr(31, 2); + + string tempH2 = imageFileName.substr(45, 2); + string tempM2 = imageFileName.substr(47, 2); + + int tempH1_V = atoi ( tempH1.c_str() ); + int tempM1_V = atoi ( tempM1.c_str() ); + int tempH2_V = atoi ( tempH2.c_str() ); + int tempM2_V = atoi ( tempM2.c_str() ); + + float tempMV_f = (tempH1_V + tempM1_V/60.0 + tempH2_V + tempM2_V/60.0) / 2; + int tempHM_V = floor ( tempMV_f ); + int tempMM_V = (tempMV_f - tempHM_V) * 60; + + timeStr = convert2NumsTo2Chars (tempHM_V) + convert2NumsTo2Chars(tempMM_V); + + printf ("dayStr=%s timeStr=%s file=%s\n", dayStr.c_str(), timeStr.c_str(), imageFileName.c_str() ); + } else { printf ("\tError: No matched Satellite file format.\n" ); @@ -1720,7 +1751,20 @@ string findOneDataFile ( string dataDir, string nameType, string dayMidStr, stri //get format if ( firstDataFile == 0 ) { - for (i=tmp_str.size()-1; i>=0; i--) + //May need to change depending on the file naming + int datePos = tmp_str.size()-1; + + if ( tmp_str.find( ".ncf" ) !=string::npos ) + { + datePos = datePos - 4 ; + } + else if ( tmp_str.find( ".nc" ) !=string::npos ) + { + datePos = datePos - 3; + } + + + for (i=datePos; i>=0; i--) { if ( isdigit(tmp_str.at(i)) == 0 ) break; //is not number j++; @@ -2319,3 +2363,71 @@ int dayofweek( string dateStr ) return day; } + + +/************************************************************************/ +/* 52. Fill day time arrays for WRF Netcdf file output */ +/************************************************************************/ +map fillTimeArrays_Day (float *times, char *timesStr, string startTime, string endTime, int timeStep, int startMins) +{ + int numTimeSteps = 0; //every timeStep minutes + int mRange0; //starting minute input; + string tmp_str; + string time_str; + long minutesPass; + string datetimeStr; //date and time string at the step + map timeIndexHash; + + //get starting minutes + string tmp_hhmm = startTime.substr(8, 4); + string tmp_date = startTime.substr(0, 8); + + //set the first time step + datetimeStr = string (startTime ); + + while (! isStepTimePassEndTime (datetimeStr, endTime) ) + { + times[numTimeSteps] = numTimeSteps * 24 * 60 * 60; + + //printf ( "\t Fill time times=%f step=%d day time string=%s\n",times[numTimeSteps],numTimeSteps,datetimeStr.c_str() ); + tmp_str = datetimeStr.substr(0, 4); //get year + time_str = tmp_str; + time_str.append ( "-" ); + + tmp_str = datetimeStr.substr(4, 2); //get month + time_str.append ( tmp_str ); + time_str.append ( "-" ); + + tmp_str = datetimeStr.substr(6, 2); //get day + time_str.append ( tmp_str ); + time_str.append ( "_" ); + + tmp_str = datetimeStr.substr(8, 2); //get hour + time_str.append ( tmp_str ); + time_str.append ( ":" ); + + tmp_str = datetimeStr.substr(10, 2); //get minute + time_str.append ( tmp_str ); + time_str.append ( ":00" ); + + strcat ( timesStr, time_str.c_str() ); //append the time string to time string array + + + //file timeIndexHash table + timeIndexHash[datetimeStr] = numTimeSteps; + + numTimeSteps++; + + + string tmp_date = datetimeStr.substr(0, 8); + + datetimeStr = getNextDayStr ( tmp_date ); + datetimeStr.append (tmp_hhmm ); + } + + printf ( "total time step = %d\n", numTimeSteps); + + + return timeIndexHash; +} +