From 2bcfbcc484bab314f6ff29faa847d26785584eb1 Mon Sep 17 00:00:00 2001 From: Melissa Sulprizio Date: Wed, 24 Apr 2024 09:19:59 -0400 Subject: [PATCH] Add option for 0.125x0.15625 resolution following updates from Xiaolin Wang The grid resolution option is expanded to include 0.125x0.15625. These modifications were obtained from Xiaolin Wang. These simulations require using the feature/0.125x0.15625_resolution in the GCClassic/src/GEOS-Chem repository (see also https://github.com/geoschem/geos-chem/pull/1980). Signed-off-by: Melissa Sulprizio --- config.yml | 3 +- .../config.harvard-cannon.12km.yml | 217 ++++++++++++++++++ .../config.harvard-cannon.global_inv.yml | 3 +- envs/Harvard-Cannon/config.harvard-cannon.yml | 3 +- src/components/prior_component/prior.sh | 8 +- src/components/setup_component/setup.sh | 7 +- .../statevector_component/aggregation.py | 9 +- .../statevector_component/statevector.sh | 12 +- src/components/template_component/template.sh | 13 +- src/inversion_scripts/imi_preview.py | 10 +- src/inversion_scripts/invert.py | 7 +- src/utilities/make_state_vector_file.py | 24 +- src/utilities/utils.py | 7 +- 13 files changed, 294 insertions(+), 29 deletions(-) create mode 100644 envs/Harvard-Cannon/config.harvard-cannon.12km.yml diff --git a/config.yml b/config.yml index aa9496c7..21426913 100644 --- a/config.yml +++ b/config.yml @@ -82,7 +82,8 @@ Gamma: 1.0 PrecomputedJacobian: false ## Grid -## Options are 0.25x0.3125 (GEOSFP only), 0.5x0.625, 2.0x2.5, or 4.0x5.0 +## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only), +## 0.5x0.625, 2.0x2.5, or 4.0x5.0 Res: "0.25x0.3125" ## Meteorology diff --git a/envs/Harvard-Cannon/config.harvard-cannon.12km.yml b/envs/Harvard-Cannon/config.harvard-cannon.12km.yml new file mode 100644 index 00000000..c87fa087 --- /dev/null +++ b/envs/Harvard-Cannon/config.harvard-cannon.12km.yml @@ -0,0 +1,217 @@ +## IMI configuration file +## Documentation @ https://imi.readthedocs.io/en/latest/getting-started/imi-config-file.html + +## General +RunName: "Test_IMI_Philadelphia" +isAWS: false +UseSlurm: true +SafeMode: true +S3Upload: false + +## Period of interest +StartDate: 20220101 +EndDate: 20230101 +SpinupMonths: 1 + +## Use blended TROPOMI+GOSAT data (true)? Or use operational TROPOMI data (false)? +BlendedTROPOMI: false + +## Is this a regional inversion? Set to false for global inversion +isRegional: true + +## Select two character region ID (for using pre-cropped meteorological fields) +## Current options are listed below with ([lat],[lon]) bounds: +## "AF" : Africa ([-37,40], [-20,53]) +## "AS" : Asia ([-11,55],[60,150]) +## "EU" : Europe ([33,61],[-30,70]) +## "ME" : Middle East ([12,50], [-20,70]) +## "NA" : North America ([10,70],[-140,-40]) +## "OC" : Oceania ([-50,5], [110,180]) +## "RU" : Russia ([41,83], [19,180]) +## "SA" : South America ([-59,16], [-88,-31]) +## "" : Use for global global simulation or custom regions +## For example, if the region of interest is in Europe ([33,61],[-30,70]), select "EU". +RegionID: "NA" + +## Region of interest +## These lat/lon bounds are only used if CreateAutomaticRectilinearStateVectorFile: true +## Otherwise lat/lon bounds are determined from StateVectorFile +LonMin: -105 +LonMax: -103 +LatMin: 31 +LatMax: 33 + +## Kalman filter options +KalmanMode: false +UpdateFreqDays: 7 +NudgeFactor: 0.1 + +## State vector +CreateAutomaticRectilinearStateVectorFile: true +nBufferClusters: 8 +BufferDeg: 5 +LandThreshold: 0.25 +OffshoreEmisThreshold: 0 +OptimizeBCs: false +OptimizeOH: false + +## Point source datasets +## Used for visualizations and state vector clustering +PointSourceDatasets: ["SRON"] + +## Clustering Options +ReducedDimensionStateVector: false +DynamicKFClustering: false +ClusteringMethod: "kmeans" +NumberOfElements: 45 +ForcedNativeResolutionElements: + - [31.5, -104] + +## Custom state vector +StateVectorFile: "/path/to/StateVector.nc" +ShapeFile: "resources/shapefiles/tmp.shp" + +## Inversion +## Note PriorError and PriorErrorOH are relative fractions (e.g. 0.5 = 50%) +## and PriorErrorBCs is in ppb +PriorError: 0.5 +PriorErrorBCs: 10.0 +PriorErrorOH: 0.5 +ObsError: 15 +Gamma: 1.0 +PrecomputedJacobian: false + +## Grid +## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only), +## 0.5x0.625, 2.0x2.5, or 4.0x5.0 +Res: "0.125x0.15625" + +## Meteorology +## Options are GEOSFP or MERRA2 +Met: "GEOSFP" + +## Setup modules +## Turn on/off different steps in setting up the inversion +RunSetup: true +SetupTemplateRundir: true +SetupSpinupRun: false +SetupJacobianRuns: false +SetupInversion: false +SetupPosteriorRun: false + +## Run modules +## Turn on/off different steps in performing the inversion +DoPriorEmis: true +DoSpinup: false +DoJacobian: false +DoInversion: false +DoPosterior: false + +## IMI preview +## NOTE: RunSetup must be true to run preview +DoPreview: true +DOFSThreshold: 0 + +## Resource allocation settings for slurm jobs +PriorMemory: 50000 +SimulationCPUs: 16 +SimulationMemory: 200000 +JacobianCPUs: 16 +JacobianMemory: 100000 +RequestedTime: "0-24:00" +SchedulerPartition: "sapphire,huce_cascade,huce_intel,seas_compute,shared" + +## Max number of simultaneous Jacobian runs from the job array (-1: no limit) +MaxSimultaneousRuns: -1 + +## Specify number of Jacobians runs to perform +## The default is -1 which will create and submit a jacobian run +## for each state vector element. Specifying a value >= 1 will +## combine state vector elements into fewer runs. +NumJacobianRuns: -1 + +##==================================================================== +## +## Advanced Settings (optional) +## +##==================================================================== + +## These settings are intended for advanced users who wish to: +## a. modify additional GEOS-Chem options, or +## b. run the IMI on a local cluster. +## They can be ignored for any standard cloud application of the IMI. + +##-------------------------------------------------------------------- +## Additional settings for GEOS-Chem simulations +##-------------------------------------------------------------------- + +## Jacobian settings +## Note PerturbValue and PerturbValueOH are relative scale factors and +## PerturbValueBCs is in ppb +PerturbValue: 1.5 +PerturbValueOH: 1.5 +PerturbValueBCs: 10.0 + +## Apply scale factors from a previous inversion? +UseEmisSF: false +UseOHSF: false + +## Save out hourly diagnostics from GEOS-Chem? +## For use in satellite operators via post-processing -- required for TROPOMI +## inversions +HourlyCH4: true + +## Turn on planeflight diagnostic in GEOS-Chem? +## For use in comparing GEOS-Chem against planeflight data. The path +## to those data must be specified in input.geos. +PLANEFLIGHT: false + +## Turn on old observation operators in GEOS-Chem? +## These will save out text files comparing GEOS-Chem to observations, but have +## to be manually incorporated into the IMI +GOSAT: false +TCCON: false +AIRS: false + +##------------------------------------------------------------------ +## Settings for running on local cluster +##------------------------------------------------------------------ + +## Path for IMI runs and output +OutputPath: "/n/holyscratch01/jacob_lab/$USER" + +## Path to GEOS-Chem input data +DataPath: "/n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/gcdata/ExtData" + +## Path to TROPOMI Data +DataPathTROPOMI: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi" + +## Conda environment file +## See envs/README to create the Conda environment specified below +CondaEnv: "imi_env" + +## GEOS-Chem environment file (with fortran compiler, netcdf libraries, etc.) +## NOTE: Copy your own file in the envs/ directory within the IMI +GEOSChemEnv: "envs/Harvard-Cannon/gcclassic.rocky+gnu12.minimal.env" + +## Download initial restart file from AWS S3? +## NOTE: Must have AWS CLI enabled +RestartDownload: false + +## Path to initial GEOS-Chem restart file + prefix +## ("YYYYMMDD_0000z.nc4" will be appended) +RestartFilePrefix: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi-boundary-conditions/v2023-10/GEOSChem.BoundaryConditions." +RestartFilePreviewPrefix: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi-boundary-conditions/v2023-10/GEOSChem.BoundaryConditions." + +## Path to GEOS-Chem boundary condition files (for regional simulations) +## BCversion will be appended to the end of this path. ${BCpath}/${BCversion} +BCpath: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi-boundary-conditions" +BCversion: "v2023-10" + +## Options to download missing GEOS-Chem input data from AWS S3 +## NOTE: Must have AWS CLI enabled +PreviewDryRun: false +SpinupDryrun: false +ProductionDryRun: false +PosteriorDryRun: false +BCdryrun: false diff --git a/envs/Harvard-Cannon/config.harvard-cannon.global_inv.yml b/envs/Harvard-Cannon/config.harvard-cannon.global_inv.yml index 9c9d4f49..7a4dafb2 100644 --- a/envs/Harvard-Cannon/config.harvard-cannon.global_inv.yml +++ b/envs/Harvard-Cannon/config.harvard-cannon.global_inv.yml @@ -82,7 +82,8 @@ Gamma: 1.0 PrecomputedJacobian: false ## Grid -## Options are 0.25x0.3125 (GEOSFP only), 0.5x0.625, 2.0x2.5, or 4.0x5.0 +## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only), +## 0.5x0.625, 2.0x2.5, or 4.0x5.0 Res: "2.0x2.5" ## Meteorology diff --git a/envs/Harvard-Cannon/config.harvard-cannon.yml b/envs/Harvard-Cannon/config.harvard-cannon.yml index 35cfd266..62a28db8 100644 --- a/envs/Harvard-Cannon/config.harvard-cannon.yml +++ b/envs/Harvard-Cannon/config.harvard-cannon.yml @@ -82,7 +82,8 @@ Gamma: 1.0 PrecomputedJacobian: false ## Grid -## Options are 0.25x0.3125 (GEOSFP only), 0.5x0.625, 2.0x2.5, or 4.0x5.0 +## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only), +## 0.5x0.625, 2.0x2.5, or 4.0x5.0 Res: "0.25x0.3125" ## Meteorology diff --git a/src/components/prior_component/prior.sh b/src/components/prior_component/prior.sh index 87095d45..1b12b480 100644 --- a/src/components/prior_component/prior.sh +++ b/src/components/prior_component/prior.sh @@ -40,10 +40,12 @@ run_prior() { resnum="3" elif [ "$Res" == "0.25x0.3125" ]; then resnum="4" + elif [ "$Res" == "0.125x0.15625" ]; then + resnum="4" else - printf "\nERROR: Grid resolution ${Res} is not supported by the IMI. " - printf "\n Options are 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n" - exit 1 + printf "\nERROR: Grid resolution ${Res} is not supported by the IMI." + printf "\n Options are 0.125x0.15625, 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n" + exit 1 fi HEMCOconfig=${RunTemplate}/HEMCO_Config.rc diff --git a/src/components/setup_component/setup.sh b/src/components/setup_component/setup.sh index e26554b4..355410bb 100644 --- a/src/components/setup_component/setup.sh +++ b/src/components/setup_component/setup.sh @@ -97,7 +97,10 @@ setup_imi() { exit 1 fi - if [ "$Res" == "0.25x0.3125" ]; then + if [ "$Res" == "0.125x0.15625" ]; then + gridDir="${Res}" + gridFile="0125x015625" + elif [ "$Res" == "0.25x0.3125" ]; then gridDir="${Res}" gridFile="025x03125" elif [ "$Res" == "0.5x0.625" ]; then @@ -111,7 +114,7 @@ setup_imi() { gridFile="4x5" else printf "\nERROR: Grid resolution ${Res} is not supported by the IMI. " - printf "\n Options are 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n" + printf "\n Options are 0.125x0.15625, 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n" exit 1 fi # Use cropped met for regional simulations instead of using global met diff --git a/src/components/statevector_component/aggregation.py b/src/components/statevector_component/aggregation.py index c3195287..a1bc4e95 100755 --- a/src/components/statevector_component/aggregation.py +++ b/src/components/statevector_component/aggregation.py @@ -224,7 +224,9 @@ def get_max_aggregation_level(config, sensitivities, desired_element_num): desired_element_num int : desired number of state vector elements Returns: int : max gridcells per cluster """ - if config["Res"] == "0.25x0.3125": + if config["Res"] == "0.125x0.15625": + max_aggregation_level = 512 + elif config["Res"] == "0.25x0.3125": max_aggregation_level = 256 elif config["Res"] == "0.5x0.625": max_aggregation_level = 64 @@ -328,7 +330,10 @@ def force_native_res_pixels(config, clusters, sensitivities): ) return sensitivities - if config["Res"] == "0.25x0.3125": + if config["Res"] == "0.125x0.15625": + lat_step = 0.125 + lon_step = 0.15625 + elif config["Res"] == "0.25x0.3125": lat_step = 0.25 lon_step = 0.3125 elif config["Res"] == "0.5x0.625": diff --git a/src/components/statevector_component/statevector.sh b/src/components/statevector_component/statevector.sh index 02a7d294..62366a5d 100644 --- a/src/components/statevector_component/statevector.sh +++ b/src/components/statevector_component/statevector.sh @@ -12,12 +12,18 @@ create_statevector() { # Use GEOS-FP or MERRA-2 CN file to determine ocean/land grid boxes if "$isRegional"; then - LandCoverFile="${DataPath}/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${RegionID}.${LandCoverFileExtension}" + if [ "$Res" = "0.125x0.15625" ]; then + LandCoverFile="/n/home00/xlwang/xlwang/methane_inversion/InputData/land_cover/data/IMERG_land_sea_mask_0125x015625.nc" + HemcoDiagFile="/n/home00/xlwang/xlwang/Test_12km/gc_0125x015625_NA_47L_geosfp_CH4_emission/OutputDir/HEMCO_diagnostics.202205010000.nc" + else + LandCoverFile="${DataPath}/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${RegionID}.${LandCoverFileExtension}" + HemcoDiagFile="${DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridFile}.20190101.nc" + fi else LandCoverFile="${DataPath}/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${LandCoverFileExtension}" + HemcoDiagFile="${DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridFile}.20190101.nc" fi - HemcoDiagFile="${DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridFile}.20190101.nc" - + if "$isAWS"; then # Download land cover and HEMCO diagnostics files s3_lc_path="s3://gcgrid/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${RegionID}.${LandCoverFileExtension}" diff --git a/src/components/template_component/template.sh b/src/components/template_component/template.sh index 8aa7a00a..60b0c9ee 100644 --- a/src/components/template_component/template.sh +++ b/src/components/template_component/template.sh @@ -53,9 +53,11 @@ setup_template() { else cmd="3\n${metNum}\n4\n1\n2\n${RunDirs}\n${runDir}\nn\n" fi + elif [ "$Res" == "0.125x0.15625" ]; then + cmd="3\n${metNum}\n5\n4\n2\n${RunDirs}\n${runDir}\nn\n" #regional run else printf "\nERROR: Grid resolution ${Res} is not supported by the IMI. " - printf "\n Options are 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n" + printf "\n Options are 0.125x0.15625, 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n" exit 1 fi @@ -115,6 +117,15 @@ setup_template() { sed -i "s/$OLD/$NEW/g" geoschem_config.yml fi + # Modify the METDIR for 0.125x0.15625 simulation + if [ "$Res" = "0.125x0.15625" ]; then + sed -i -e "s:GEOS_0.25x0.3125\/GEOS_FP:GEOS_0.25x0.3125_NA\/GEOS_FP:g" HEMCO_Config.rc.gmao_metfields_0125 + OLD="/n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/gcdata/ExtData/GEOS_0.125x0.15625/GEOS_FP" + NEW="/n/holylfs05/LABS/jacob_lab/Users/xlwang/methane_inversion/InputData/GEOS_0.125x0.15625_NA/GEOS_FP_DerivedWinds" + sed -i "s|$OLD|$NEW|g" HEMCO_Config.rc.gmao_metfields_0125 + sed -i '/METDIR/d' HEMCO_Config.rc + fi + # Modify HEMCO_Config.rc based on settings in config.yml # Use cropped met fields (add the region to both METDIR and the met files) if "$isRegional"; then diff --git a/src/inversion_scripts/imi_preview.py b/src/inversion_scripts/imi_preview.py index 2d5e5457..72047e08 100755 --- a/src/inversion_scripts/imi_preview.py +++ b/src/inversion_scripts/imi_preview.py @@ -165,7 +165,9 @@ def imi_preview( ] inversion_area_km = calculate_area_in_km(coords) - if config["Res"] == "0.25x0.3125": + if config["Res"] == "0.125x0.15625": + res_factor = 2 + elif config["Res"] == "0.25x0.3125": res_factor = 1 elif config["Res"] == "0.5x0.625": res_factor = 0.5 @@ -492,7 +494,11 @@ def estimate_averaging_kernel( # Set resolution specific variables # L_native = Rough length scale of native state vector element [m] - if config["Res"] == "0.25x0.3125": + if config["Res"] == "0.125x0.15625": + L_native = 12 * 1000 + lat_step = 0.125 + lon_step = 0.15625 + elif config["Res"] == "0.25x0.3125": L_native = 25 * 1000 lat_step = 0.25 lon_step = 0.3125 diff --git a/src/inversion_scripts/invert.py b/src/inversion_scripts/invert.py index 6e49e923..1bcd48a2 100644 --- a/src/inversion_scripts/invert.py +++ b/src/inversion_scripts/invert.py @@ -56,8 +56,11 @@ def do_inversion( # Need to ignore data in the GEOS-Chem 3 3 3 3 buffer zone # Shave off one or two degrees of latitude/longitude from each side of the domain # ~1 degree if 0.25x0.3125 resolution, ~2 degrees if 0.5x0.6125 resolution - # This assumes 0.25x0.3125 and 0.5x0.625 simulations are always regional - if "0.25x0.3125" in res: + # This assumes 0.125x0.15625, 0.25x0.3125, and 0.5x0.625 simulations are always regional + if "0.125x0.15625" in res: + degx = 4 * 0.15625 + degy = 4 * 0.125 + elif "0.25x0.3125" in res: degx = 4 * 0.3125 degy = 4 * 0.25 elif "0.5x0.625" in res: diff --git a/src/utilities/make_state_vector_file.py b/src/utilities/make_state_vector_file.py index 525e31df..833daf50 100644 --- a/src/utilities/make_state_vector_file.py +++ b/src/utilities/make_state_vector_file.py @@ -131,16 +131,20 @@ def make_state_vector_file( lc = xr.load_dataset(land_cover_pth) hd = xr.load_dataset(hemco_diag_pth) - # Require hemco diags on same global grid as land cover map - # TODO remove this offset once the HEMCO standalone files - # are regenerated with recent bugfix that corrects the offset - if config["Res"] == "0.25x0.3125": - hd["lon"] = hd["lon"] - 0.03125 - elif config["Res"] == "0.5x0.625": - hd["lon"] = hd["lon"] - 0.0625 - - # Select / group fields together - lc = (lc["FRLAKE"] + lc["FRLAND"] + lc["FRLANDIC"]).drop("time").squeeze() + # Group together + if config['Res'] == '0.125x0.15625': + lc = lc["landseamask"] #100% = all water and 0% = all land + lc = np.round( -(lc/100.-1), decimals=5) + else: + # Require hemco diags on same global grid as land cover map + # TODO remove this offset once the HEMCO standalone files + # are regenerated with recent bugfix that corrects the offset + if config["Res"] == "0.25x0.3125": + hd["lon"] = hd["lon"] - 0.03125 + elif config["Res"] == "0.5x0.625": + hd["lon"] = hd["lon"] - 0.0625 + lc = (lc["FRLAKE"] + lc["FRLAND"] + lc["FRLANDIC"]).drop("time").squeeze() + hd = (hd["EmisCH4_Oil"] + hd["EmisCH4_Gas"]).drop("time").squeeze() # Check compatibility of region of interest diff --git a/src/utilities/utils.py b/src/utilities/utils.py index ea6e4ba3..7a8bf232 100644 --- a/src/utilities/utils.py +++ b/src/utilities/utils.py @@ -28,7 +28,10 @@ def download_landcover_files(config): elif config["Res"] == "0.25x0.3125": gridDir= "0.25x0.3125" gridFile= "025x03125" - + elif config["Res"] == "0.125x0.15625": + gridDir= "0.125x0.15625" + gridFile= "0125x015625" + if len(config["RegionID"]) == 2: LandCoverFile = f"{DataPath}/GEOS_{gridDir}_{config['RegionID']}/{metDir}/{constYr}/01/{config['Met']}.{constYr}0101.CN.{gridFile}.{config['RegionID']}.{LandCoverFileExtension}" s3_lc_path = f"s3://gcgrid/GEOS_{gridDir}_{config['RegionID']}/{metDir}/{constYr}/01/{config['Met']}.{constYr}0101.CN.{gridFile}.{config['RegionID']}.{LandCoverFileExtension}" @@ -62,6 +65,8 @@ def download_hemcodiags_files(config): gridFile = "05x0625" elif config["Res"] == "0.25x0.3125": gridFile = "025x03125" + elif config["Res"] == "0.125x0.15625": + gridFile = "0125x015625" HemcoDiagFile = f"{DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.{gridFile}.20190101.nc" s3_hd_path = f"s3://gcgrid/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.{gridFile}.20190101.nc"