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"