From f281ba83f4a3ccb243650386a1db36ba7ad6be8c Mon Sep 17 00:00:00 2001 From: djvaron Date: Mon, 29 May 2023 17:37:38 -0400 Subject: [PATCH] Feature/offshore emissions (#109) * Implement inclusion of offshore emissions in automatic state vector. Including new config variables allowing threshold for emissions to be included in sv. --- .gitignore | 3 + config.yml | 1 + .../getting-started/imi-config-file.rst | 2 + .../statevector_component/statevector.sh | 10 +- .../statevector_from_shapefile.ipynb | 45 +++++---- src/utilities/make_state_vector_file.py | 94 +++++++++---------- src/utilities/sanitize_input_yaml.py | 1 + src/utilities/utils.py | 34 ++++++- 8 files changed, 117 insertions(+), 73 deletions(-) diff --git a/.gitignore b/.gitignore index 7a12c949..e924251e 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,9 @@ GEOS-Chem GCClassic docs/build/ .DS_STORE +src/inversion_scripts/__pycache__ +src/utilities/__pycache__ +src/notebooks/.ipynb_checkpoints __pycache__/ .ipynb_checkpoints/ slurm-*.out diff --git a/config.yml b/config.yml index 28a99e95..cfd2c1e2 100644 --- a/config.yml +++ b/config.yml @@ -43,6 +43,7 @@ CreateAutomaticRectilinearStateVectorFile: true nBufferClusters: 8 BufferDeg: 5 LandThreshold: 0.25 +OffshoreEmisThreshold: 0 ## Clustering Options ReducedDimensionStateVector: false diff --git a/docs/source/getting-started/imi-config-file.rst b/docs/source/getting-started/imi-config-file.rst index 75819104..c1ae5098 100644 --- a/docs/source/getting-started/imi-config-file.rst +++ b/docs/source/getting-started/imi-config-file.rst @@ -66,6 +66,8 @@ State vector - Width of the buffer elements, in degrees; will not be used if ``CreateAutomaticRectilinearStateVectorFile`` is ``false``. Default is ``5`` (~500 km). * - ``LandThreshold`` - Land-cover fraction below which to exclude GEOS-Chem grid cells from the state vector when creating the state vector file. Default value is ``0.25``. + * - ``OffshoreEmisThreshold`` + - Offshore GEOS-Chem grid cells with oil/gas emissions above this threshold will be included in the state vector. Default value is ``0``. Clustering Options ^^^^^^^^^^^^^^^^^^ diff --git a/src/components/statevector_component/statevector.sh b/src/components/statevector_component/statevector.sh index 2bd73b9c..bd1fca35 100644 --- a/src/components/statevector_component/statevector.sh +++ b/src/components/statevector_component/statevector.sh @@ -12,11 +12,14 @@ create_statevector() { # Use GEOS-FP or MERRA-2 CN file to determine ocean/land grid boxes LandCoverFile="${DataPath}/GEOS_${gridDir}/${metDir}/${constYr}/01/${metUC}.${constYr}0101.CN.${gridRes}.${NestedRegion}.${LandCoverFileExtension}" + HemcoDiagFile="${DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridRes}.20190101.nc" if "$isAWS"; then - # Download land cover file + # Download land cover and Hemco diagnostics files s3_lc_path="s3://gcgrid/GEOS_${gridDir}/${metDir}/${constYr}/01/${metUC}.${constYr}0101.CN.${gridRes}.${NestedRegion}.${LandCoverFileExtension}" aws s3 cp --request-payer=requester ${s3_lc_path} ${LandCoverFile} + s3_hd_path="s3://gcgrid/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridRes}.20190101.nc" + aws s3 cp --request-payer=requester ${s3_hd_path} ${HemcoDiagFile} fi # Output path and filename for state vector file @@ -29,8 +32,11 @@ create_statevector() { cp ${InversionPath}/src/utilities/make_state_vector_file.py . chmod 755 make_state_vector_file.py + # Get config path + config_path=${InversionPath}/${ConfigFile} + printf "\nCalling make_state_vector_file.py\n" - python make_state_vector_file.py $LandCoverFile $StateVectorFName $LatMin $LatMax $LonMin $LonMax $BufferDeg $LandThreshold $nBufferClusters + python make_state_vector_file.py $config_path $LandCoverFile $HemcoDiagFile $StateVectorFName printf "\n=== DONE CREATING RECTANGULAR STATE VECTOR FILE ===\n" } diff --git a/src/notebooks/statevector_from_shapefile.ipynb b/src/notebooks/statevector_from_shapefile.ipynb index e11409a0..8b9baa84 100644 --- a/src/notebooks/statevector_from_shapefile.ipynb +++ b/src/notebooks/statevector_from_shapefile.ipynb @@ -31,6 +31,7 @@ "import sys\n", "import xarray as xr \n", "import numpy as np\n", + "import fiona\n", "import regionmask\n", "import xesmf as xe\n", "import geopandas as gpd\n", @@ -42,7 +43,7 @@ "import yaml\n", "import warnings; warnings.filterwarnings(action='ignore')\n", "sys.path.append('/home/ubuntu/integrated_methane_inversion/src/')\n", - "from utilities.utils import download_landcover_files\n", + "from utilities.utils import download_landcover_files, download_hemcodiags_files\n", "%matplotlib inline" ] }, @@ -78,7 +79,8 @@ ], "source": [ "# Download landcover files from s3 if not already present\n", - "download_landcover_files(config)" + "download_landcover_files(config)\n", + "download_hemcodiags_files(config)" ] }, { @@ -178,7 +180,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAADnCAYAAACKYEokAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa3ElEQVR4nO3daZBV1b338e/fBg2DBJBGCUhAAyISBSEgaAShClHuoz6OINdEnC4qopjSOCUqmmiRq1GJCIpDNE7YV40PNJh6wImEBkFRaJAIRKSBVIvNkCg00P2/L86hPUAPp7v3mfb+fap24Z7XedE/19prrb3N3RERCaNDMl0AEZFUUcCJSGgp4EQktBRwIhJaCjgRCa0mte0cMWKEb9myJV1lEZEaLF269G13H9GYa5x5Rgv/uqwiuft9Wt7o+2WDWgNuy5YtLFmyJF1lEZEamFm7xl7j67IKFr/dOalj8zp83uj7ZYNaA05EwsOBSiozXYy0UsCJRITj7PHkmqhhoYATiRDV4EQklBynImJTMxVwIhFSiQJORELIgQoFnIiElWpwIhJKDuzRMzgRCSPH1UQVkZByqIhWvingRKIiNpMhWhRwIpFhVGCZLkRaKeBEIiLWyaCAE5EQio2DU8CJSEhVqgYnImGkGpyIhJZjVETsKwXR+rUiEVfpltRSFzN7xsxKzWxFDfvNzB4zszVm9qmZnXzA/jwz+9jMZgX006qlgBOJCMfY7XlJLUl4Dqjtmw1nAd3iyzXAEwfsvxFY1YCfUS8KOJGIiA30PSSppc5rub8PlNVyyLnA8x5TBLQ2sw4AZtYJGAnMaPyvqp2ewYlESD06GdqZWeIXp5509yfrcauOwIaE9ZL4ts3AI8CtwOH1uF6DKOBEIsLdqPCkG21b3L1fI25XXZK6mf0HUOruS81sSCOunxQ1UUUipBJLaglACXB0wnonYBNwKnCOmX0BvAIMNbM/BXHD6ijgRCIi1snQJKklAG8BP4v3pp4CbHf3ze5+u7t3cvcuwChgvrv/ZxA3rI6aqCIRsa+TIQhm9jIwhNizuhLgbqApgLtPAwqBs4E1wLfA2EBuXE8KOJEIqQhoqpa7j65jvwPX13HMu8C7gRSoBgo4kYiI4kwGBZxIhFQm34saCgo4kYiITbZXwIlICDnGnuSmYYWGAk4kItypz0DfUFDAiURGYIN4c4YCTiQiHNXgRCTE1MkgIqHkJPcyyzBRwIlEROyzgdH6k4/WrxWJNH34WURCytFMBhEJMdXgRCSU3E01OBEJp1gng6ZqiUgo1eubDKGggBOJiFgng57BiUhIaSaDiISSZjKISKgF9dGZXKGAE4kId9hTqYATkRCKNVEVcCISUprJICKhpGEiIhJiaqKKSIjpmwwiEkqxXlTNRRWRENJAXxEJNTVRRSSU1IsqIqGmXlQRCSV3Y68CTkTCSk1UEQmlKD6Di1Z9VSTiKt2SWupiZs+YWamZrahhv5nZY2a2xsw+NbOT49uPNrN3zGyVmRWb2Y0B/8T9KOBEImLfOLggAg54DhhRy/6zgG7x5Rrgifj2vcAv3P144BTgejPr2eAfVQcFnEiEVGJJLXVx9/eBsloOORd43mOKgNZm1sHdN7v7R/Fr/AtYBXQM4KdVS8/gRCLCHfam74WXHYENCesl8W2b920wsy5AH2BRqgqhgBOJkHp0MrQzsyUJ60+6+5P1uFV1N/KqnWYtgf8BbnL3HfW4br0o4EQiop5zUbe4e79G3K4EODphvROwCcDMmhILtxfd/fVG3KNOegYnEiHultQSgLeAn8V7U08Btrv7ZjMz4Glglbs/HMSNaqManEiEBDXZ3sxeBoYQa8qWAHcDTQHcfRpQCJwNrAG+BcbGTz0VuAxYbmbL4tvucPfCQAp2AAWcSES4BzfQ191H17Hfgeur2b6A6p/PpYQCTiQyjAp9NlBEwiqg52s5QwEnEhFRnIuqgBOJCo89h4sSBVwWO+axhzJdhP2sm/CLTBdBGkmvLJeMy7Zg26c+5VIYZh9XJ4NkWraGW33t+x0KuuyiJqpkTFjCTbJX1HpRo1VfzWJhDbew/q5c5J7WqVpZQTW4LKAQkHTRMBFJK4WbpJOewYlIKDlGpXpRJV1Ue5N0i1gFTgGXKQo3STuPXi+qAi4DohRuGgeXZSJWhVPApVmUwk2yj2pwkhIKNsk0ByorFXASMIWbZAUHVIOToCjYJNtoHJw0moJNspYCThpD4fYd9aBmm3DNM02GAi5ACjfJeqrBiUgoObh6UaUhVHuT3KCAE5GwUhNVpPHUwZClFHBSX2qeSk7I0YG+ZtYduAX4IQmZ5e5D6zpXASeBU+0te+XoQN/XgGnAU0BFfU5UwIlESW72ou519ycacqICTgKjmlv2s9yswf0/M7sOeAMo37fR3cvqOlEBJ4FQuOUAJ1c7GX4e//eWhG0OHFPXiQo4aRQFWy6xnOxkcPeuDT1XAScNpnDLQTlYgzOzpsC1wOnxTe8C0919T13nRusTOymiP3TJGZVJLtnlCaAvMDW+9I1vq5NqcNJg+8b/KeBzRI6OgwN+4u4nJazPN7NPkjlRNThptGMee0iDnXOEeXJLndcxe8bMSs1sRQ37zcweM7M1ZvapmZ2csG+Ema2O77stiWJXmNmxCecfQ5Lj4RRwEhgFXQ7wJJe6PQeMqGX/WUC3+HIN8SalmeUBj8f39wRGm1nPOu51C/COmb1rZu8B84Gkmg0KuIComfYdhVz4ufv7QG3j0M4FnveYIqC1mXUA+gNr3H2du+8GXokfW9u95hELygnx5Th3fyeZcirgJCUUctmpHk3Udma2JGG5pp636ghsSFgviW+rafvBZTUbGv/3fGAk8CPgWGBkfFud1MkgEhVOfaZqbXH3fo24W3U38lq2V2cwsebo/6nhnNfrKoQCLkDrJvxCNZcExzz2kJru2SZ94+BKgKMT1jsBm4BDa9h+EHe/O/7v2IYWQgEnEiFpnIv6FjDezF4BBgDb3X2zmX0FdDOzrsBGYBRwaW0XMrObq9m8HVjq7stqO1fP4AKmGotktYB6Uc3sZWAhcJyZlZjZlWY2zszGxQ8pBNYBa4i95ug6AHffC4wH3gZWATPdvbiO2/UDxvHdM7xrgCHAU2Z2a20nqgaXAmqqfkfN1CwTUA3O3UfXsd+B62vYV0gsAJN1BHCyu/8bwMzuBgqITd1aCkyu6UQFnEhEJDuINwt1BnYnrO8BfujuO82svIZzADVRU0a1FslKlZbckl1eAorM7O547e2vwMtm1gJYWduJCjhJOTXXs0dQU7XSyd3vA64GthHrXBjn7pPc/Rt3H1PbuQq4FFItTrJOcFO10q0ZsMPdHwHWx3th66SASzGFnGSNJGtv2VaDizdLfwncHt/UFPhTMucq4NJAIadmatbIzRrc/wXOAb4BcPdNwOHJnKiASxOFnGQDq0xuyTK748NOHCDeuZAUBVwarZvwCwWdSD2YmQGzzGw6sTeSXA38f2KDh+ukgMuAqIacmqlZIMeaqPGa23nEBvb+D3Ac8Gt3n5LM+RromyH7Qk5/9JI2WdiBkKSFwDZ3v6XOIw+ggMswBZ2kVW4G3BnAf5nZeuIdDQDufmJdJyrgskRUgk5zUzMsNwPurIaeqIDLMol//GEPO0kvIyt7SOvk7usbeq4CLosFXdPJlsBULS5DcvcZXIMp4CJEtUPJ0SZqgyngIioqz/zkABELOI2Di7hMDT5WsGZGLs5FbQzV4ARQjS4yQhReyVANTvajh/8h5jk7F7XBFHByEIVciOXYVK3GUsBJtRRy4RS1Z3AKOKlRqjsg9LwvAyJWg1Mng9RJ4+dCImThlQzV4KRe1HRNv1/+8peBXMeIXhNVNTipNw0pSR93Z/LkGr9rXG9hCq9kqAYnDaY3FKeemfGXv/wluAtG7BmcAk4aTSGXWrt27QruYhELODVRJRAN6YhQMCbnpJNOCuZCIXu+lgwFnASuuuDS87qG69y5c3AXU8CJBE+1tewQpmlYyVDAiWS5rVu3BnYtNVFFJCvMmjWLmTNnsm7dumAuGLIOhGQo4ESy1Jw5c3jhhReCvWjEAk7DRESyVN++fWnVqlVg14viTAYFnEgWeuedd/jtb3/LxIkT2bNnT2DXtUpPaqnzOmYjzGy1ma0xs9uq2d/GzN4ws0/NbLGZ9UrYN9HMis1shZm9bGbfC+wHHkABJ5Jlpk6dymWXXcZDDz3EPffcQ5MmAT1JSnaQbx35ZmZ5wOPEvlfaExhtZj0POOwOYFn848w/Ax6Nn9sRmAD0c/deQB4wqvE/rnp6BieSJebNm8fkyZNZv3497733Hscee2zg9wio+dkfWOPu6wDM7BXgXGBlwjE9gQcA3P0zM+tiZkfG9zUBmpnZHqA5sCmQUlVDNTiRDPvyyy8ZPXo0V155JaNHj2bJkiUpCTcgqKlaHYENCesl8W2JPgHOBzCz/sAPgU7uvhH4b+BLYDOw3d0DnGy7PwWcSCN98MEHzJgxg/Lycq699lrMjK+//rraY8vLy5k9ezY33XQTgwcP5uKLL6Z3794cd9xxFBcXc/nll9OyZcuUlbUenQztzGxJwnJN4mWqufSBsfgg0MbMlgE3AB8De82sDbHaXlfgB0ALM/vPgH9mFTVRRRpp4sSJLF26lA8//JBDDz0UgHHjxlFQUADAxo0badq0KQUFBVx33XX7nfvwww/z+OOPk5+fn57CJt9E3eLu/WrYVwIcnbDeiQOame6+AxgLYGYG/CO+nAn8w92/iu97HRgE/CnpktWDAk6kkaZMmcKgQYN46qmncI8lSEFBAd26dePzzz+nY8cDW28wevRopk+fzuGHH56+gnpgU7U+BLqZWVdgI7FOgksTDzCz1sC37r4buAp43913mNmXwClm1hzYCQwDlgRSqmqoiSrSSCeddBJt2rSpCrfbb7+d3bt389prrx107IwZMygrK+Oll15Kb7gR3Dg4d98LjAfeBlYBM9292MzGmdm4+GHHA8Vm9hmx3tYb4+cuAgqAj4DlxDLoyeB/bYxqcCINUF5eznnnncfvfvc7evXqxcKFC+nRowcA7du354033uCSSy4BoEuXLrz55pvBvfaoMTyYblR3LwQKD9g2LeG/FwLdajj3buDuQApSBwWcSAM88sgjzJ07l7lz57Jy5Uq6du3KqFGj+Nvf/sZdd93FN998Q7Nmzfjkk0/o1q3av/OMCNMshWQo4EQaYNiwYQC0adOGsWPHYmb8/e9/p6ysDIj1rJ522mmZLOLBIjjZXs/gROpp5syZtG3blg4dOjB9+nSOOOIIioqKqsLtwQcfzL5wi7PK5JawUA1OpA7uzs6dO/nss8+49957eeuttzj//PPZvHkz99xzDytXxgbwv/rqq1x44YUcckj21hvCFF7JUMCJ1GH8+PFMnTp1v22vv/46ACtXrmT8+PHce++9tG3bNhPFS54TWCdDrlDAidSirKyMqVOn0qpVK3bs2AHEeklLS0sB+Otf/8qgQYMyWcR6iVonQ/bWpUWyQKtWrbjllluqwm3IkCGUlpbSo0cPioqKcircAH02UES+06RJEx588EG++uornnvuOd59910g1kvarl27zBaunvYN9I0S1eBEEkybNo0jjzySjh07YmYMHz6cyZMnM2vWrKqhIStWrMi5cAPAk3vZZTIvvMwVCjiRBG3atKG0tJRNm2Jzx1u2bMntt9/OxIkTmTdvHvfddx8nnHBChkvZCBFroirgRBJccskluDvXXnstQ4cOZfz48axatYo777wToOrfXBW1bzLoGZxINR566CGmTZtW1SwF2LJlC7E3/+QoB0LU/EyGanAi1WjWrBmdOnWiefPmALzwwgscccQRGS5VACLWRFUNTqQa7s6kSZP49ttvOf744xkzZkymixSIMDU/k6EanETaokWLePTRR6ve5bbP4sWL+eKLL4BYz2pON00TRK0XVTU4ibTNmzdz0003sWLFCp566inmzZvHl19+SWFhIf/+979p27Ytp59+eqaLGYyQNT+ToYCTSDvnnHOA2Jt2ly5dyscff7zf/ilTpmSiWCkRG+gbrYRTE1Ui7eWXX67678RwO/LI2Cc8zz333LSXKaUqk1xCQjU4ibRRo0bRokULNm3aREFBAUVFRezcuZNt27bRtWtXWrRokekiBipqNbhaA279+vU899xz9OnThyVLlnDmmWfSqVOndJVNJOXy8vI477zzALjuuuuoqKhg0KBBLF68mIsuuiizhQuansHtr1mzZhQWFjJ27Fgg9maF4uJihZyE1syZM1m8eDFA6GpvEK4e0mTU+gyuffv2zJw5k40bNzJs2DB27NhR9aI/kbBZvXo1d9xxB8cccwwAe/fuzXCJUsA9uSUkkupk+MEPfsDbb78NwI033pjSAolkytq1a9m2bRsLFiwAYm8NCRWP3jcZku5FzcvLY9asWbRu3fqgQZEiYdCiRQu2bdvGwIEDAQ4aMhIKqsHVbOTIkbRs2TJ8/2cTIfbMGaBnz54AVTMZQiVic1HrPQ7u1ltv5ZJLLmHt2rXs3r2bFStWYGYMGjSILVu2pKKMImnRv3//qilaV111FQA7d+7McKmCZZWVSS1hUe+Au+GGG7jqqqsYMGAAhx12GD/+8Y8BWLhwIfn5+VRUVAReSJF06dOnD2VlZQwfPhyAuXPnZrhEAXIiN9C3QTMZbr75ZkpKSiguLsbdWbVqFVOmTKF79+5MnDiRPXv2BF1OkbRo0qQJP/3pT9m1axcAl19+eWYLFCDDMU9uCYsGT9X63ve+V/WsokePHowfP54PPviA5cuXM2HCBHVESM7Kz8/ntttuq/qaVqiaqepkaLj27dvz5ptvMmfOnKrBkiK55oknnmDTpk3ceuutADzzzDMZLlGAFHCN8/3vf5/TTjuN++67j/Xr1wd9eZGUu+GGGzjqqKOq1n/9619nsDQB0jO4YLz44ovMnj2bLl26qNNBcs7vf/97vv76a/Lz8xkzZgxlZWWZLlJg1IvaSGvXrt1v/Yorrgj6FiIplZeXxwsvvADEvrIFUFJSkskiBSTJ5qmaqDXr27cv8N2LAlu3bh30LURS7uKLL6Znz55s374diLVKcp6jgGuMzz//nO3bt9O/f3+uv/56du/ezcMPPxzkLUTSwszo378/K1euBODpp5/OcIkComdwDbNt2zaGDx/O9OnTWbRoEWZG06ZNycvLC+oWIml19dVX88ADDwDfveE31wU1Ds7MRpjZajNbY2a3VbO/jZm9YWafmtliM+uVsK+1mRWY2WdmtsrMBgb8M6sEFnBTpkzhjDPO4JprrgnqkiIZNWjQILp27QpQ9YaRnBdAE9XM8oDHgbOAnsBoM+t5wGF3AMvc/UTgZ8CjCfseBea6ew/gJGBVQL/uIIEF3NatW+nWrVtQlxPJCs8++yx//OMfadmyZe4PXneHisrkltr1B9a4+zp33w28Ahz48YqewLzYbf0zoIuZHWlmrYDTgafj+3a7+7YAf+V+Agu43r178/777wd1OZGsMHjwYI466ij69u0bjm+jJl+Da2dmSxKWxKZZR2BDwnpJfFuiT4DzAcysP/BDoBNwDPAV8KyZfWxmM8wsZa9ODizgVq9eTevWrfnqq69qPOaOO+6gsLAwqFuKpMWGDRvo2PHAv98clXzAbXH3fgnLkwlXqS7pD6zePgi0MbNlwA3Ax8BeYp9JOBl4wt37AN8ABz3DC0pgAXfqqafyyiuv0LlzZ5599tmD9r/++us88MADbNu2LahbiqTF3r17adq0aaaL0XgOVHpyS+1KgKMT1jsBm/a7lfsOdx/r7r2JPYPLB/4RP7fE3RfFDy0gFngpEVjA7XsL6q5du7jiiito27Ytq1ev5l//+hc///nPueCCCwC49NJLg7qlSFrMnDmTYcOGZboYAXDwyuSW2n0IdDOzrmZ2KDAKeCvxgHhP6aHx1auA9+Oh909gg5kdF983DFgZ3G/cX2DfRW3evDmnnnoqS5Ysoby8nK1bt9KjR4/9jpk2bRrffvstzZs3D+q2Iim1du1aVqxYwYUXXpjpojSek0wHQt2Xcd9rZuOBt4E84Bl3LzazcfH904DjgefNrIJYgF2ZcIkbgBfjAbgOGNvoQtUgsBrcYYcdxoIFC9i1axfuzpVXXnnQ/nHjxpGfn8+cOXOCuq1ISv3hD3/g6quvrnqdec4LaCaDuxe6e3d3P9bdfxPfNi0ebrj7Qnfv5u493P18d9+acO6y+HO9E939vMR9QUvJZHuAGTNm8MEHH1Stl5eXM3ToUI499limT5+eqtuKBKa0tJTnn38+XGM7NVUrOKeddhrLly+vWp8/fz7Lly/nz3/+M2ZGUVFRKm8v0ij3338/Y8aMoXPnzpkuSkA02T5wvXr1wt1xd9555x3OOuusqn0DBw5k1qxZeqWSZJ2FCxfy2muv8atf/SrTRQmOA5WVyS0hkfKASzRkyBAKCwtxdyorKykoKGDSpEl07tyZiy66iLFjx4bza+KSU/aNBHjsscfIz8/PdHGCFbEaXGC9qPVlZlxwwQVccMEFLFq0iFNOOQWAk08+mcGDB3PiiSdmqmgScZMmTaJXr15cdNFFmS5KwDyQXtRckrGAS7Tv04MAEyZMoGXLlvTu3Ztp06ZxwgknZLBkEjV79+5l+vTpIf2qPXjdY9xCJa1N1Jo0b96cysrKqoe53bt3Z8GCBfTr109TuyStioqK6Ny5c4g6Fg4QzEyGnJEVNTiINVm/+OILunfvzkcffQTEnoWMHDky99/iIDmjsLCQs88+O9PFSJ2I/S1lRQ1uHzNj0aLYFLXLLruM9u3bA7ExdCLpEOqAc1cvaqa1bduWyZMns3r1akpLS4HYLAiRVNu4cSMbNmyo6vAKpYj1omZdwAHcfPPNVZP3AYqLizNYGomKOXPmMHz48BC/Zt/xioqklrDIyoDLy8vjkUce4Sc/+UnVukiqzZ49m5EjR2a6GKkT3OuSckZWBtw+jz76KD/60Y8oKipi69aUzccVoby8nPnz53PmmWdmuiipFczrknJGVgfcgAED6NOnD2PHjqVt27bcddddmtYlKbFgwQKOP/748M1cSOCAV3pSS1hkdcAdcsghvPrqq9x5550A/OY3v2HYsGEUFRXRoUMHevXqxZ49ezJcSgmDUPee7uOBvfAyZ2R1wEFs6Mj9999f1aP63nvvMXDgQP75z39SXFzMSy+9lOESShgsW7aMAQMGZLoYKRe1TgarbRCtmS1JY1lEpGZb3H1EYy5gZnOBdum6XzaoNeBERHJZ1jdRRUQaSgEnIqGlgBOR0FLAiUhoKeBEJLT+F+jsaoqD9zL0AAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAADnCAYAAACKYEokAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa3ElEQVR4nO3daZBV1b338e/fBg2DBJBGCUhAAyISBSEgaAShClHuoz6OINdEnC4qopjSOCUqmmiRq1GJCIpDNE7YV40PNJh6wImEBkFRaJAIRKSBVIvNkCg00P2/L86hPUAPp7v3mfb+fap24Z7XedE/19prrb3N3RERCaNDMl0AEZFUUcCJSGgp4EQktBRwIhJaCjgRCa0mte0cMWKEb9myJV1lEZEaLF269G13H9GYa5x5Rgv/uqwiuft9Wt7o+2WDWgNuy5YtLFmyJF1lEZEamFm7xl7j67IKFr/dOalj8zp83uj7ZYNaA05EwsOBSiozXYy0UsCJRITj7PHkmqhhoYATiRDV4EQklBynImJTMxVwIhFSiQJORELIgQoFnIiElWpwIhJKDuzRMzgRCSPH1UQVkZByqIhWvingRKIiNpMhWhRwIpFhVGCZLkRaKeBEIiLWyaCAE5EQio2DU8CJSEhVqgYnImGkGpyIhJZjVETsKwXR+rUiEVfpltRSFzN7xsxKzWxFDfvNzB4zszVm9qmZnXzA/jwz+9jMZgX006qlgBOJCMfY7XlJLUl4Dqjtmw1nAd3iyzXAEwfsvxFY1YCfUS8KOJGIiA30PSSppc5rub8PlNVyyLnA8x5TBLQ2sw4AZtYJGAnMaPyvqp2ewYlESD06GdqZWeIXp5509yfrcauOwIaE9ZL4ts3AI8CtwOH1uF6DKOBEIsLdqPCkG21b3L1fI25XXZK6mf0HUOruS81sSCOunxQ1UUUipBJLaglACXB0wnonYBNwKnCOmX0BvAIMNbM/BXHD6ijgRCIi1snQJKklAG8BP4v3pp4CbHf3ze5+u7t3cvcuwChgvrv/ZxA3rI6aqCIRsa+TIQhm9jIwhNizuhLgbqApgLtPAwqBs4E1wLfA2EBuXE8KOJEIqQhoqpa7j65jvwPX13HMu8C7gRSoBgo4kYiI4kwGBZxIhFQm34saCgo4kYiITbZXwIlICDnGnuSmYYWGAk4kItypz0DfUFDAiURGYIN4c4YCTiQiHNXgRCTE1MkgIqHkJPcyyzBRwIlEROyzgdH6k4/WrxWJNH34WURCytFMBhEJMdXgRCSU3E01OBEJp1gng6ZqiUgo1eubDKGggBOJiFgng57BiUhIaSaDiISSZjKISKgF9dGZXKGAE4kId9hTqYATkRCKNVEVcCISUprJICKhpGEiIhJiaqKKSIjpmwwiEkqxXlTNRRWRENJAXxEJNTVRRSSU1IsqIqGmXlQRCSV3Y68CTkTCSk1UEQmlKD6Di1Z9VSTiKt2SWupiZs+YWamZrahhv5nZY2a2xsw+NbOT49uPNrN3zGyVmRWb2Y0B/8T9KOBEImLfOLggAg54DhhRy/6zgG7x5Rrgifj2vcAv3P144BTgejPr2eAfVQcFnEiEVGJJLXVx9/eBsloOORd43mOKgNZm1sHdN7v7R/Fr/AtYBXQM4KdVS8/gRCLCHfam74WXHYENCesl8W2b920wsy5AH2BRqgqhgBOJkHp0MrQzsyUJ60+6+5P1uFV1N/KqnWYtgf8BbnL3HfW4br0o4EQiop5zUbe4e79G3K4EODphvROwCcDMmhILtxfd/fVG3KNOegYnEiHultQSgLeAn8V7U08Btrv7ZjMz4Glglbs/HMSNaqManEiEBDXZ3sxeBoYQa8qWAHcDTQHcfRpQCJwNrAG+BcbGTz0VuAxYbmbL4tvucPfCQAp2AAWcSES4BzfQ191H17Hfgeur2b6A6p/PpYQCTiQyjAp9NlBEwiqg52s5QwEnEhFRnIuqgBOJCo89h4sSBVwWO+axhzJdhP2sm/CLTBdBGkmvLJeMy7Zg26c+5VIYZh9XJ4NkWraGW33t+x0KuuyiJqpkTFjCTbJX1HpRo1VfzWJhDbew/q5c5J7WqVpZQTW4LKAQkHTRMBFJK4WbpJOewYlIKDlGpXpRJV1Ue5N0i1gFTgGXKQo3STuPXi+qAi4DohRuGgeXZSJWhVPApVmUwk2yj2pwkhIKNsk0ByorFXASMIWbZAUHVIOToCjYJNtoHJw0moJNspYCThpD4fYd9aBmm3DNM02GAi5ACjfJeqrBiUgoObh6UaUhVHuT3KCAE5GwUhNVpPHUwZClFHBSX2qeSk7I0YG+ZtYduAX4IQmZ5e5D6zpXASeBU+0te+XoQN/XgGnAU0BFfU5UwIlESW72ou519ycacqICTgKjmlv2s9yswf0/M7sOeAMo37fR3cvqOlEBJ4FQuOUAJ1c7GX4e//eWhG0OHFPXiQo4aRQFWy6xnOxkcPeuDT1XAScNpnDLQTlYgzOzpsC1wOnxTe8C0919T13nRusTOymiP3TJGZVJLtnlCaAvMDW+9I1vq5NqcNJg+8b/KeBzRI6OgwN+4u4nJazPN7NPkjlRNThptGMee0iDnXOEeXJLndcxe8bMSs1sRQ37zcweM7M1ZvapmZ2csG+Ema2O77stiWJXmNmxCecfQ5Lj4RRwEhgFXQ7wJJe6PQeMqGX/WUC3+HIN8SalmeUBj8f39wRGm1nPOu51C/COmb1rZu8B84Gkmg0KuIComfYdhVz4ufv7QG3j0M4FnveYIqC1mXUA+gNr3H2du+8GXokfW9u95hELygnx5Th3fyeZcirgJCUUctmpHk3Udma2JGG5pp636ghsSFgviW+rafvBZTUbGv/3fGAk8CPgWGBkfFud1MkgEhVOfaZqbXH3fo24W3U38lq2V2cwsebo/6nhnNfrKoQCLkDrJvxCNZcExzz2kJru2SZ94+BKgKMT1jsBm4BDa9h+EHe/O/7v2IYWQgEnEiFpnIv6FjDezF4BBgDb3X2zmX0FdDOzrsBGYBRwaW0XMrObq9m8HVjq7stqO1fP4AKmGotktYB6Uc3sZWAhcJyZlZjZlWY2zszGxQ8pBNYBa4i95ug6AHffC4wH3gZWATPdvbiO2/UDxvHdM7xrgCHAU2Z2a20nqgaXAmqqfkfN1CwTUA3O3UfXsd+B62vYV0gsAJN1BHCyu/8bwMzuBgqITd1aCkyu6UQFnEhEJDuINwt1BnYnrO8BfujuO82svIZzADVRU0a1FslKlZbckl1eAorM7O547e2vwMtm1gJYWduJCjhJOTXXs0dQU7XSyd3vA64GthHrXBjn7pPc/Rt3H1PbuQq4FFItTrJOcFO10q0ZsMPdHwHWx3th66SASzGFnGSNJGtv2VaDizdLfwncHt/UFPhTMucq4NJAIadmatbIzRrc/wXOAb4BcPdNwOHJnKiASxOFnGQDq0xuyTK748NOHCDeuZAUBVwarZvwCwWdSD2YmQGzzGw6sTeSXA38f2KDh+ukgMuAqIacmqlZIMeaqPGa23nEBvb+D3Ac8Gt3n5LM+RromyH7Qk5/9JI2WdiBkKSFwDZ3v6XOIw+ggMswBZ2kVW4G3BnAf5nZeuIdDQDufmJdJyrgskRUgk5zUzMsNwPurIaeqIDLMol//GEPO0kvIyt7SOvk7usbeq4CLosFXdPJlsBULS5DcvcZXIMp4CJEtUPJ0SZqgyngIioqz/zkABELOI2Di7hMDT5WsGZGLs5FbQzV4ARQjS4yQhReyVANTvajh/8h5jk7F7XBFHByEIVciOXYVK3GUsBJtRRy4RS1Z3AKOKlRqjsg9LwvAyJWg1Mng9RJ4+dCImThlQzV4KRe1HRNv1/+8peBXMeIXhNVNTipNw0pSR93Z/LkGr9rXG9hCq9kqAYnDaY3FKeemfGXv/wluAtG7BmcAk4aTSGXWrt27QruYhELODVRJRAN6YhQMCbnpJNOCuZCIXu+lgwFnASuuuDS87qG69y5c3AXU8CJBE+1tewQpmlYyVDAiWS5rVu3BnYtNVFFJCvMmjWLmTNnsm7dumAuGLIOhGQo4ESy1Jw5c3jhhReCvWjEAk7DRESyVN++fWnVqlVg14viTAYFnEgWeuedd/jtb3/LxIkT2bNnT2DXtUpPaqnzOmYjzGy1ma0xs9uq2d/GzN4ws0/NbLGZ9UrYN9HMis1shZm9bGbfC+wHHkABJ5Jlpk6dymWXXcZDDz3EPffcQ5MmAT1JSnaQbx35ZmZ5wOPEvlfaExhtZj0POOwOYFn848w/Ax6Nn9sRmAD0c/deQB4wqvE/rnp6BieSJebNm8fkyZNZv3497733Hscee2zg9wio+dkfWOPu6wDM7BXgXGBlwjE9gQcA3P0zM+tiZkfG9zUBmpnZHqA5sCmQUlVDNTiRDPvyyy8ZPXo0V155JaNHj2bJkiUpCTcgqKlaHYENCesl8W2JPgHOBzCz/sAPgU7uvhH4b+BLYDOw3d0DnGy7PwWcSCN98MEHzJgxg/Lycq699lrMjK+//rraY8vLy5k9ezY33XQTgwcP5uKLL6Z3794cd9xxFBcXc/nll9OyZcuUlbUenQztzGxJwnJN4mWqufSBsfgg0MbMlgE3AB8De82sDbHaXlfgB0ALM/vPgH9mFTVRRRpp4sSJLF26lA8//JBDDz0UgHHjxlFQUADAxo0badq0KQUFBVx33XX7nfvwww/z+OOPk5+fn57CJt9E3eLu/WrYVwIcnbDeiQOame6+AxgLYGYG/CO+nAn8w92/iu97HRgE/CnpktWDAk6kkaZMmcKgQYN46qmncI8lSEFBAd26dePzzz+nY8cDW28wevRopk+fzuGHH56+gnpgU7U+BLqZWVdgI7FOgksTDzCz1sC37r4buAp43913mNmXwClm1hzYCQwDlgRSqmqoiSrSSCeddBJt2rSpCrfbb7+d3bt389prrx107IwZMygrK+Oll15Kb7gR3Dg4d98LjAfeBlYBM9292MzGmdm4+GHHA8Vm9hmx3tYb4+cuAgqAj4DlxDLoyeB/bYxqcCINUF5eznnnncfvfvc7evXqxcKFC+nRowcA7du354033uCSSy4BoEuXLrz55pvBvfaoMTyYblR3LwQKD9g2LeG/FwLdajj3buDuQApSBwWcSAM88sgjzJ07l7lz57Jy5Uq6du3KqFGj+Nvf/sZdd93FN998Q7Nmzfjkk0/o1q3av/OMCNMshWQo4EQaYNiwYQC0adOGsWPHYmb8/e9/p6ysDIj1rJ522mmZLOLBIjjZXs/gROpp5syZtG3blg4dOjB9+nSOOOIIioqKqsLtwQcfzL5wi7PK5JawUA1OpA7uzs6dO/nss8+49957eeuttzj//PPZvHkz99xzDytXxgbwv/rqq1x44YUcckj21hvCFF7JUMCJ1GH8+PFMnTp1v22vv/46ACtXrmT8+PHce++9tG3bNhPFS54TWCdDrlDAidSirKyMqVOn0qpVK3bs2AHEeklLS0sB+Otf/8qgQYMyWcR6iVonQ/bWpUWyQKtWrbjllluqwm3IkCGUlpbSo0cPioqKcircAH02UES+06RJEx588EG++uornnvuOd59910g1kvarl27zBaunvYN9I0S1eBEEkybNo0jjzySjh07YmYMHz6cyZMnM2vWrKqhIStWrMi5cAPAk3vZZTIvvMwVCjiRBG3atKG0tJRNm2Jzx1u2bMntt9/OxIkTmTdvHvfddx8nnHBChkvZCBFroirgRBJccskluDvXXnstQ4cOZfz48axatYo777wToOrfXBW1bzLoGZxINR566CGmTZtW1SwF2LJlC7E3/+QoB0LU/EyGanAi1WjWrBmdOnWiefPmALzwwgscccQRGS5VACLWRFUNTqQa7s6kSZP49ttvOf744xkzZkymixSIMDU/k6EanETaokWLePTRR6ve5bbP4sWL+eKLL4BYz2pON00TRK0XVTU4ibTNmzdz0003sWLFCp566inmzZvHl19+SWFhIf/+979p27Ytp59+eqaLGYyQNT+ToYCTSDvnnHOA2Jt2ly5dyscff7zf/ilTpmSiWCkRG+gbrYRTE1Ui7eWXX67678RwO/LI2Cc8zz333LSXKaUqk1xCQjU4ibRRo0bRokULNm3aREFBAUVFRezcuZNt27bRtWtXWrRokekiBipqNbhaA279+vU899xz9OnThyVLlnDmmWfSqVOndJVNJOXy8vI477zzALjuuuuoqKhg0KBBLF68mIsuuiizhQuansHtr1mzZhQWFjJ27Fgg9maF4uJihZyE1syZM1m8eDFA6GpvEK4e0mTU+gyuffv2zJw5k40bNzJs2DB27NhR9aI/kbBZvXo1d9xxB8cccwwAe/fuzXCJUsA9uSUkkupk+MEPfsDbb78NwI033pjSAolkytq1a9m2bRsLFiwAYm8NCRWP3jcZku5FzcvLY9asWbRu3fqgQZEiYdCiRQu2bdvGwIEDAQ4aMhIKqsHVbOTIkbRs2TJ8/2cTIfbMGaBnz54AVTMZQiVic1HrPQ7u1ltv5ZJLLmHt2rXs3r2bFStWYGYMGjSILVu2pKKMImnRv3//qilaV111FQA7d+7McKmCZZWVSS1hUe+Au+GGG7jqqqsYMGAAhx12GD/+8Y8BWLhwIfn5+VRUVAReSJF06dOnD2VlZQwfPhyAuXPnZrhEAXIiN9C3QTMZbr75ZkpKSiguLsbdWbVqFVOmTKF79+5MnDiRPXv2BF1OkbRo0qQJP/3pT9m1axcAl19+eWYLFCDDMU9uCYsGT9X63ve+V/WsokePHowfP54PPviA5cuXM2HCBHVESM7Kz8/ntttuq/qaVqiaqepkaLj27dvz5ptvMmfOnKrBkiK55oknnmDTpk3ceuutADzzzDMZLlGAFHCN8/3vf5/TTjuN++67j/Xr1wd9eZGUu+GGGzjqqKOq1n/9619nsDQB0jO4YLz44ovMnj2bLl26qNNBcs7vf/97vv76a/Lz8xkzZgxlZWWZLlJg1IvaSGvXrt1v/Yorrgj6FiIplZeXxwsvvADEvrIFUFJSkskiBSTJ5qmaqDXr27cv8N2LAlu3bh30LURS7uKLL6Znz55s374diLVKcp6jgGuMzz//nO3bt9O/f3+uv/56du/ezcMPPxzkLUTSwszo378/K1euBODpp5/OcIkComdwDbNt2zaGDx/O9OnTWbRoEWZG06ZNycvLC+oWIml19dVX88ADDwDfveE31wU1Ds7MRpjZajNbY2a3VbO/jZm9YWafmtliM+uVsK+1mRWY2WdmtsrMBgb8M6sEFnBTpkzhjDPO4JprrgnqkiIZNWjQILp27QpQ9YaRnBdAE9XM8oDHgbOAnsBoM+t5wGF3AMvc/UTgZ8CjCfseBea6ew/gJGBVQL/uIIEF3NatW+nWrVtQlxPJCs8++yx//OMfadmyZe4PXneHisrkltr1B9a4+zp33w28Ahz48YqewLzYbf0zoIuZHWlmrYDTgafj+3a7+7YAf+V+Agu43r178/777wd1OZGsMHjwYI466ij69u0bjm+jJl+Da2dmSxKWxKZZR2BDwnpJfFuiT4DzAcysP/BDoBNwDPAV8KyZfWxmM8wsZa9ODizgVq9eTevWrfnqq69qPOaOO+6gsLAwqFuKpMWGDRvo2PHAv98clXzAbXH3fgnLkwlXqS7pD6zePgi0MbNlwA3Ax8BeYp9JOBl4wt37AN8ABz3DC0pgAXfqqafyyiuv0LlzZ5599tmD9r/++us88MADbNu2LahbiqTF3r17adq0aaaL0XgOVHpyS+1KgKMT1jsBm/a7lfsOdx/r7r2JPYPLB/4RP7fE3RfFDy0gFngpEVjA7XsL6q5du7jiiito27Ytq1ev5l//+hc///nPueCCCwC49NJLg7qlSFrMnDmTYcOGZboYAXDwyuSW2n0IdDOzrmZ2KDAKeCvxgHhP6aHx1auA9+Oh909gg5kdF983DFgZ3G/cX2DfRW3evDmnnnoqS5Ysoby8nK1bt9KjR4/9jpk2bRrffvstzZs3D+q2Iim1du1aVqxYwYUXXpjpojSek0wHQt2Xcd9rZuOBt4E84Bl3LzazcfH904DjgefNrIJYgF2ZcIkbgBfjAbgOGNvoQtUgsBrcYYcdxoIFC9i1axfuzpVXXnnQ/nHjxpGfn8+cOXOCuq1ISv3hD3/g6quvrnqdec4LaCaDuxe6e3d3P9bdfxPfNi0ebrj7Qnfv5u493P18d9+acO6y+HO9E939vMR9QUvJZHuAGTNm8MEHH1Stl5eXM3ToUI499limT5+eqtuKBKa0tJTnn38+XGM7NVUrOKeddhrLly+vWp8/fz7Lly/nz3/+M2ZGUVFRKm8v0ij3338/Y8aMoXPnzpkuSkA02T5wvXr1wt1xd9555x3OOuusqn0DBw5k1qxZeqWSZJ2FCxfy2muv8atf/SrTRQmOA5WVyS0hkfKASzRkyBAKCwtxdyorKykoKGDSpEl07tyZiy66iLFjx4bza+KSU/aNBHjsscfIz8/PdHGCFbEaXGC9qPVlZlxwwQVccMEFLFq0iFNOOQWAk08+mcGDB3PiiSdmqmgScZMmTaJXr15cdNFFmS5KwDyQXtRckrGAS7Tv04MAEyZMoGXLlvTu3Ztp06ZxwgknZLBkEjV79+5l+vTpIf2qPXjdY9xCJa1N1Jo0b96cysrKqoe53bt3Z8GCBfTr109TuyStioqK6Ny5c4g6Fg4QzEyGnJEVNTiINVm/+OILunfvzkcffQTEnoWMHDky99/iIDmjsLCQs88+O9PFSJ2I/S1lRQ1uHzNj0aLYFLXLLruM9u3bA7ExdCLpEOqAc1cvaqa1bduWyZMns3r1akpLS4HYLAiRVNu4cSMbNmyo6vAKpYj1omZdwAHcfPPNVZP3AYqLizNYGomKOXPmMHz48BC/Zt/xioqklrDIyoDLy8vjkUce4Sc/+UnVukiqzZ49m5EjR2a6GKkT3OuSckZWBtw+jz76KD/60Y8oKipi69aUzccVoby8nPnz53PmmWdmuiipFczrknJGVgfcgAED6NOnD2PHjqVt27bcddddmtYlKbFgwQKOP/748M1cSOCAV3pSS1hkdcAdcsghvPrqq9x5550A/OY3v2HYsGEUFRXRoUMHevXqxZ49ezJcSgmDUPee7uOBvfAyZ2R1wEFs6Mj9999f1aP63nvvMXDgQP75z39SXFzMSy+9lOESShgsW7aMAQMGZLoYKRe1TgarbRCtmS1JY1lEpGZb3H1EYy5gZnOBdum6XzaoNeBERHJZ1jdRRUQaSgEnIqGlgBOR0FLAiUhoKeBEJLT+F+jsaoqD9zL0AAAAAElFTkSuQmCC", "text/plain": [ "
" ] @@ -263,7 +265,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAADnCAYAAAApbXvLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPp0lEQVR4nO3dbVAV1R8H8O/yJIU4BlohJIoIKEmlgECRozSCZIo2E9WEk5pMlpYSmcw0OaX51IT6lxxK4kVpzKiMoQlpT0oaOsJEGSJQMCFh5pUHUbhwgf2/YLyFwuVyuXt37/H7eQV7F/Y3Dl/POXvOnpVkWQYRicFB7QKIyHoYaCKBMNBEAmGgiQTCQBMJxMnUh3FxcbJOp7NVLURkppKSkqOyLMfdetxkoHU6HYqLi5WriogsIknSqL6Os8tNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSiJPaBdDg+P3vQ6v9rurX3rDa7yJtYAtNJBAGmkgg7HJrmDW712r8/luxi688ttBEAmGgiQTCQBMJhGNoDbH1mJbEwxaaSCAMNJFAGGgigXAMrRKOl0kJbKGJBMJAEwmEXW4bYRebbIEtNJFAGGgigTDQRALhGFoBHC+TWhhoUgyff7Y9drmJBMJAEwmEgSYSCANNJBAGmkggDDSRQDhtRVbBKSptYAtNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJpK7KKvh4Z5VSW7bGFJhIIA00kEHa5STG3dsPZBVceA62Avv5wuS0R2QK73EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQDgPbSO3zk1zXpqUwBaaSCAMNJFA2OVWCZeHkhLYQhMJhIEmEgi73GQz3NVEeQy0hpjzx23uOJvTZHcmdrmJBMJAEwmEgSYSCMfQdsbSm0jWHJ+TdrGFJhIIA00kEAaaSCAcQ5MR15fbP7bQRAJhoIkEwi43mcTpLu35448/+v2MLTSRxjU1NUGSJGzevBmtra1ISEjo91wGmkjjuru7AQBpaWkYN24cHBz6jy0DTaRxHh4e+PDDDzFy5Ehs2rQJhw8f7vdcjqFpyDjOVt7Vq1cxYsQIhISEwNPTs9/z2EIT2YGUlBTU1tYiPDwcw4YN6/c8BprIDnh6eiIjIwMA4Ozs3O95DDSRnQgMDBzwHI6hySa4d9jQZGdnY+nSpXBzc8ONGzf6PY8tNJGGRUZGQpIkpKamAgBu3LiBjRs39ns+A02kYX5+fnB1dUVjYyOGDRuGixcvIi0trd/zGWgijZJlGa6urvD09ERWVhb0ej18fHxM/gzH0EQadPnyZbzwwgu4evUqfvzxR4wfP96sn2MLTaQhsizj559/xowZMxAWFobTp0+bHWaALTSRZlRUVOCVV15BWVkZNmzYgJdeemnQv4MtNJGVVVRU4MKFC2afX1tbi9WrV+Oxxx7DU089hbq6OovCDDDQRFaxdetWODk5oaCgANu3b8ekSZOwf/9+HDx4EHv27EFHRwcAoL29Hbt27cK0adMgSRIkSUJQUBA6OztRVlaGVatWwcnJ8o4zu9xEViBJErq6urBo0SLce++9AIBnnnnG+Hl2djaio6Px3nvv9fq52NhYFBQUQJIk69Qhy3K/H4aGhsrFxcVWuRCRyLq7u+Ho6Njr2IIFCzB//ny8+OKLt50/d+5cZGZmwtvb26LrSZJUIsty6K3H2eUmGgKdTgeDwQAHBwd8/vnnxuNVVVXIzc2Fm5tbr/NLS0shyzIOHz5scZhNYZebyEI6nQ6jR48GABgMBsTFxRk/CwoKQldXF4CehyoOHTqEgIAAxWtiC01koVOnThm/fuedd1BRUYE5c+bAzc3NGOazZ8/iwoULNgkzwDE0kUVu5mbGjBnw8fFBTk4OPDw80NDQYDzHYDAM6Y61KRxDE1lJbW0tPDw88Oabb6Kurg6PPPIIgoKCjGHetWsXmpubFQuzKWyhicxw6dIlLF68GI6OjsjPzzced3V1hV6vBwCsWrUK77//Pu6++27F6+mvheZNMSIzvPXWWzh69Ohtx/V6PYKCgvDtt98qctd6sBhoIjNUVFT0eXzatGk4c+bMbXPQauEYmsgMZ86cMU5LzZo1CwCwfv16FBcXaybMAFtoottcu3YNBQUFcHd3R3NzM1xcXBAaGgofHx/Mnz8feXl58PPzw9tvv612qbcxGeibr+AgupMUFRXh2WefNX7v4+ODuro6+Pn5wd3dHQBQXl6uVnkmmexym3qHDpGoYmNjUVhYiLCwMDzxxBMoLS3FypUrMW3aNPzyyy8oLCyEi4uL2mX2iYkl6kN0dDSOHj0KNzc3+Pr6IiAgAPv378fatWsRHR2tdnn94jw0kQnt7e3w8PBAa2srAKCrq0sTPVeuFCOywN69e43vkvr11181EWZTtF0dkQ3Isoxz5871eXzz5s1obGxEQkICpkyZokJ1g8NA0x1Pp9MhJCQE8+bNA9Azu1NaWopjx47hn3/+AQBs2LBBzRLNZnLa6uLFi7aqg0g1o0ePRnh4OA4fPowtW7Zg7dq1AABvb280NzcjKCgIwcHBKldpHpM3xSRJkk+dOoWoqCgblkRkWy0tLRgxYkS/n+fk5PSal9YCi26KOTk5YcWKFWhpaVGuMiKVDR8+HIsXLzYuGrlVTEyMjSuynMlAP/TQQwgLC8P06dOxZcsW/PDDD2hra0NjY2OvB7mJ7JkkScjOzsa1a9fQ2dmJ6upqbNy40bhG++Y2Q/ZgwHnos2fP4tixY/j6669RVFSEc+fOobW1Fc7Ozvj000+RlJRkw3KJlNfR0WGcqgL+3Z1ESyyeh5YkCbGxsdi2bRtOnz6N33//HUlJSTAYDFi0aJEy1RKpRJZlfPTRR8bvtT7vfKtBV+vl5YXPPvsMycnJAIDz589bvSgitRgMBqSkpCA7OxsTJ060uweULP7vJzMzE0DP6hkiUTg7OwMAlixZgsDAQJWrGTyLAy1JEjIyMlBQUGDNeohU1dnZafz65pZC/z2mdUMaICQkJKCgoKBXqPft24fJkyezK052ydnZGZWVlXB3dzeuHPvtt99Ursp8Q9qxxNvbG3l5eViwYAFGjhyJlpYW1NfXAwCCg4PR2NiIkSNHWqNOIpuZOHEioqKijH/LmZmZxiGm1g35Fl5kZCT+/PNPbN26Fbm5uTh79izmzJkDFxcXxMbGoqmpyQplEtlWfHw8SkpKMGHCBHz88cdql2M2q9yTHzZsGObNm4eIiAiEhoYiPz8fer0e4eHhiIqK4kozsjuRkZHIzMzEJ598AgAoKytTuSLzKDbJJkkSdu7ciUcffRRpaWlKXYZIEevWrQMATJ8+HQDw8ssvq1mO2RSfNV+8eDG+//571NbWKn0pIqtJSkqCJEmor69HdHQ0Tp48qXZJZlE80DU1NSgvL4evry+qqqqUvhyRVSQmJmL27NkICAgwvrDdHp5fUDzQq1evNn4dEBCA9vZ2pS9JNGQODg7Ys2cPgH8Xm+Tk5KhZklkUDXR5eTmuXLkCAJgwYQIA9Fr0TqRlo0aNQmJiIoqKigAA6enpKlc0MEUDvX79egBAfn4+Kisr8ddffyl5OSKrmz17Nurq6jB16lRUV1erXc6AFHsVzo4dO1BaWoqGhgbcc889AIAxY8YodTkiRcTExGDZsmXo7u7Gfffdh/b2dk33MhVpoVtaWrBu3TocOXLEGGYie+Tr64vXXnsNAHD58mXN39hVJNDXrl2Di4sLxo0bp8SvJ7KpZcuWIT09HQ888AAkSVK7HJMUCfSYMWPQ3d2t2Rd6EQ3G5MmTsWTJEjQ0NGj+kUpFAt3Y2IirV6/CwcGh3+1bjh8/zhVkZDcuXryI+++/H05O2n4DsyKBdnNzAwBMmjQJjz/+OAwGQ6/PDQYDZs6caVeL3unO1tnZqfkwAwoF+r93AU+ePAkXFxdkZWUB6HkU7eYeyBUVFUpcnsjq9u/fj1mzZqldxoAUm4fOzc2Fk5MTXF1dAfTcWJAkCcuXL4der8fcuXOh1+uVujyR1XR0dODAgQN28YCGYoFeuHAhDAYD2tra+nwR2FdffYWxY8di+fLlSpVAZBUHDx6Er68vQkJC1C5lQDbZo/TBBx9EW1sbZs6c2ev41q1bkZmZia6uLluUQTRo3d3d2LRpE1599VW1SzGLzTYddnV1xTfffIOnn37aeGzNmjUAel65k5ycrMkNzenO9sUXXxg38LAHNt1F3NHREQcOHIAsy9DpdL2exNq9ezcWLlyIS5cu2bIkon5duXIFqamp2LFjh+YXlNyk2msBPD09kZ6eDlmWIcsyampqMGrUKAQHByMiIgL+/v6oqalRqzwirFixAklJSYiIiFC7FLMN+G6r4uJiG5YDXL9+3fgWwJiYGCxduhTPPfecTWsgysvLw5o1a1BaWoq77rpL7XJuY/G7rWxt+PDhmDp1KgDgu+++Q3JyMry8vLBt2zaVK6M7SVZWFtatW6fJMJuiuUADQElJCbZv3w6gp8X++++/kZKSgpUrV/LGGSlOr9fjxIkTiIuLU7uUQdNkoAHg9ddfx969e3sdy8jIQHNzs0oV0Z3ixIkTCAkJgYeHh9qlDJpmAw0Azz//PPz8/ODv74/U1FQAPYvkiZSUn5+PJ598Uu0yLKLpQAPAl19+iYaGBuzbtw9Az+tsiZQiyzKOHDmC+Ph4tUuxiOYDPWXKFOTk5Bj39c7NzVW5IhJZVVUV2tra7GKZZ180H2igZ6O23bt3A4Dd/kOTfcjPz0d8fLzdLCS5lV0EGuh5k4G7uzveffddHDp0SO1ySFBHjhyx2/EzMMDCkoCAALmpqQnV1dUYPny4Dcvq239vVkiShKqqKuN+30RDdf36dXh5eaG+vt64uEmrLFpY0tzcjCtXriA8PFwTzy7Hx8ejsrISQM/NC39/f6Snp0OSJEiSxDvgNCTHjx9HeHi45sNsislA37yjXF5ejjNnztikoIFMnDgRsiwjMjISAPDGG28AAAIDA/HBBx+oWRrZufPnz+Phhx9Wu4whMRloJycnJCYmYuzYsdDpdJpapfXTTz+hsLDQ+H1FRQV27typYkVk72pqajB+/Hi1yxgSk2NoSZJs+2QGEZlLJ8vybWtTTQaaiOyL3UxbEdHAGGgigTDQRAJhoIkEwkATCeT/lMOEy78WhG0AAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAADnCAYAAAApbXvLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPp0lEQVR4nO3dbVAV1R8H8O/yJIU4BlohJIoIKEmlgECRozSCZIo2E9WEk5pMlpYSmcw0OaX51IT6lxxK4kVpzKiMoQlpT0oaOsJEGSJQMCFh5pUHUbhwgf2/YLyFwuVyuXt37/H7eQV7F/Y3Dl/POXvOnpVkWQYRicFB7QKIyHoYaCKBMNBEAmGgiQTCQBMJxMnUh3FxcbJOp7NVLURkppKSkqOyLMfdetxkoHU6HYqLi5WriogsIknSqL6Os8tNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJAEwmEgSYSiJPaBdDg+P3vQ6v9rurX3rDa7yJtYAtNJBAGmkgg7HJrmDW712r8/luxi688ttBEAmGgiQTCQBMJhGNoDbH1mJbEwxaaSCAMNJFAGGgigXAMrRKOl0kJbKGJBMJAEwmEXW4bYRebbIEtNJFAGGgigTDQRALhGFoBHC+TWhhoUgyff7Y9drmJBMJAEwmEgSYSCANNJBAGmkggDDSRQDhtRVbBKSptYAtNJBAGmkggDDSRQBhoIoEw0EQCYaCJBMJpK7KKvh4Z5VSW7bGFJhIIA00kEHa5STG3dsPZBVceA62Avv5wuS0R2QK73EQCYaCJBMJAEwmEgSYSCANNJBAGmkggDDSRQDgPbSO3zk1zXpqUwBaaSCAMNJFA2OVWCZeHkhLYQhMJhIEmEgi73GQz3NVEeQy0hpjzx23uOJvTZHcmdrmJBMJAEwmEgSYSCMfQdsbSm0jWHJ+TdrGFJhIIA00kEAaaSCAcQ5MR15fbP7bQRAJhoIkEwi43mcTpLu35448/+v2MLTSRxjU1NUGSJGzevBmtra1ISEjo91wGmkjjuru7AQBpaWkYN24cHBz6jy0DTaRxHh4e+PDDDzFy5Ehs2rQJhw8f7vdcjqFpyDjOVt7Vq1cxYsQIhISEwNPTs9/z2EIT2YGUlBTU1tYiPDwcw4YN6/c8BprIDnh6eiIjIwMA4Ozs3O95DDSRnQgMDBzwHI6hySa4d9jQZGdnY+nSpXBzc8ONGzf6PY8tNJGGRUZGQpIkpKamAgBu3LiBjRs39ns+A02kYX5+fnB1dUVjYyOGDRuGixcvIi0trd/zGWgijZJlGa6urvD09ERWVhb0ej18fHxM/gzH0EQadPnyZbzwwgu4evUqfvzxR4wfP96sn2MLTaQhsizj559/xowZMxAWFobTp0+bHWaALTSRZlRUVOCVV15BWVkZNmzYgJdeemnQv4MtNJGVVVRU4MKFC2afX1tbi9WrV+Oxxx7DU089hbq6OovCDDDQRFaxdetWODk5oaCgANu3b8ekSZOwf/9+HDx4EHv27EFHRwcAoL29Hbt27cK0adMgSRIkSUJQUBA6OztRVlaGVatWwcnJ8o4zu9xEViBJErq6urBo0SLce++9AIBnnnnG+Hl2djaio6Px3nvv9fq52NhYFBQUQJIk69Qhy3K/H4aGhsrFxcVWuRCRyLq7u+Ho6Njr2IIFCzB//ny8+OKLt50/d+5cZGZmwtvb26LrSZJUIsty6K3H2eUmGgKdTgeDwQAHBwd8/vnnxuNVVVXIzc2Fm5tbr/NLS0shyzIOHz5scZhNYZebyEI6nQ6jR48GABgMBsTFxRk/CwoKQldXF4CehyoOHTqEgIAAxWtiC01koVOnThm/fuedd1BRUYE5c+bAzc3NGOazZ8/iwoULNgkzwDE0kUVu5mbGjBnw8fFBTk4OPDw80NDQYDzHYDAM6Y61KRxDE1lJbW0tPDw88Oabb6Kurg6PPPIIgoKCjGHetWsXmpubFQuzKWyhicxw6dIlLF68GI6OjsjPzzced3V1hV6vBwCsWrUK77//Pu6++27F6+mvheZNMSIzvPXWWzh69Ohtx/V6PYKCgvDtt98qctd6sBhoIjNUVFT0eXzatGk4c+bMbXPQauEYmsgMZ86cMU5LzZo1CwCwfv16FBcXaybMAFtoottcu3YNBQUFcHd3R3NzM1xcXBAaGgofHx/Mnz8feXl58PPzw9tvv612qbcxGeibr+AgupMUFRXh2WefNX7v4+ODuro6+Pn5wd3dHQBQXl6uVnkmmexym3qHDpGoYmNjUVhYiLCwMDzxxBMoLS3FypUrMW3aNPzyyy8oLCyEi4uL2mX2iYkl6kN0dDSOHj0KNzc3+Pr6IiAgAPv378fatWsRHR2tdnn94jw0kQnt7e3w8PBAa2srAKCrq0sTPVeuFCOywN69e43vkvr11181EWZTtF0dkQ3Isoxz5871eXzz5s1obGxEQkICpkyZokJ1g8NA0x1Pp9MhJCQE8+bNA9Azu1NaWopjx47hn3/+AQBs2LBBzRLNZnLa6uLFi7aqg0g1o0ePRnh4OA4fPowtW7Zg7dq1AABvb280NzcjKCgIwcHBKldpHpM3xSRJkk+dOoWoqCgblkRkWy0tLRgxYkS/n+fk5PSal9YCi26KOTk5YcWKFWhpaVGuMiKVDR8+HIsXLzYuGrlVTEyMjSuynMlAP/TQQwgLC8P06dOxZcsW/PDDD2hra0NjY2OvB7mJ7JkkScjOzsa1a9fQ2dmJ6upqbNy40bhG++Y2Q/ZgwHnos2fP4tixY/j6669RVFSEc+fOobW1Fc7Ozvj000+RlJRkw3KJlNfR0WGcqgL+3Z1ESyyeh5YkCbGxsdi2bRtOnz6N33//HUlJSTAYDFi0aJEy1RKpRJZlfPTRR8bvtT7vfKtBV+vl5YXPPvsMycnJAIDz589bvSgitRgMBqSkpCA7OxsTJ060uweULP7vJzMzE0DP6hkiUTg7OwMAlixZgsDAQJWrGTyLAy1JEjIyMlBQUGDNeohU1dnZafz65pZC/z2mdUMaICQkJKCgoKBXqPft24fJkyezK052ydnZGZWVlXB3dzeuHPvtt99Ursp8Q9qxxNvbG3l5eViwYAFGjhyJlpYW1NfXAwCCg4PR2NiIkSNHWqNOIpuZOHEioqKijH/LmZmZxiGm1g35Fl5kZCT+/PNPbN26Fbm5uTh79izmzJkDFxcXxMbGoqmpyQplEtlWfHw8SkpKMGHCBHz88cdql2M2q9yTHzZsGObNm4eIiAiEhoYiPz8fer0e4eHhiIqK4kozsjuRkZHIzMzEJ598AgAoKytTuSLzKDbJJkkSdu7ciUcffRRpaWlKXYZIEevWrQMATJ8+HQDw8ssvq1mO2RSfNV+8eDG+//571NbWKn0pIqtJSkqCJEmor69HdHQ0Tp48qXZJZlE80DU1NSgvL4evry+qqqqUvhyRVSQmJmL27NkICAgwvrDdHp5fUDzQq1evNn4dEBCA9vZ2pS9JNGQODg7Ys2cPgH8Xm+Tk5KhZklkUDXR5eTmuXLkCAJgwYQIA9Fr0TqRlo0aNQmJiIoqKigAA6enpKlc0MEUDvX79egBAfn4+Kisr8ddffyl5OSKrmz17Nurq6jB16lRUV1erXc6AFHsVzo4dO1BaWoqGhgbcc889AIAxY8YodTkiRcTExGDZsmXo7u7Gfffdh/b2dk33MhVpoVtaWrBu3TocOXLEGGYie+Tr64vXXnsNAHD58mXN39hVJNDXrl2Di4sLxo0bp8SvJ7KpZcuWIT09HQ888AAkSVK7HJMUCfSYMWPQ3d2t2Rd6EQ3G5MmTsWTJEjQ0NGj+kUpFAt3Y2IirV6/CwcGh3+1bjh8/zhVkZDcuXryI+++/H05O2n4DsyKBdnNzAwBMmjQJjz/+OAwGQ6/PDQYDZs6caVeL3unO1tnZqfkwAwoF+r93AU+ePAkXFxdkZWUB6HkU7eYeyBUVFUpcnsjq9u/fj1mzZqldxoAUm4fOzc2Fk5MTXF1dAfTcWJAkCcuXL4der8fcuXOh1+uVujyR1XR0dODAgQN28YCGYoFeuHAhDAYD2tra+nwR2FdffYWxY8di+fLlSpVAZBUHDx6Er68vQkJC1C5lQDbZo/TBBx9EW1sbZs6c2ev41q1bkZmZia6uLluUQTRo3d3d2LRpE1599VW1SzGLzTYddnV1xTfffIOnn37aeGzNmjUAel65k5ycrMkNzenO9sUXXxg38LAHNt1F3NHREQcOHIAsy9DpdL2exNq9ezcWLlyIS5cu2bIkon5duXIFqamp2LFjh+YXlNyk2msBPD09kZ6eDlmWIcsyampqMGrUKAQHByMiIgL+/v6oqalRqzwirFixAklJSYiIiFC7FLMN+G6r4uJiG5YDXL9+3fgWwJiYGCxduhTPPfecTWsgysvLw5o1a1BaWoq77rpL7XJuY/G7rWxt+PDhmDp1KgDgu+++Q3JyMry8vLBt2zaVK6M7SVZWFtatW6fJMJuiuUADQElJCbZv3w6gp8X++++/kZKSgpUrV/LGGSlOr9fjxIkTiIuLU7uUQdNkoAHg9ddfx969e3sdy8jIQHNzs0oV0Z3ixIkTCAkJgYeHh9qlDJpmAw0Azz//PPz8/ODv74/U1FQAPYvkiZSUn5+PJ598Uu0yLKLpQAPAl19+iYaGBuzbtw9Az+tsiZQiyzKOHDmC+Ph4tUuxiOYDPWXKFOTk5Bj39c7NzVW5IhJZVVUV2tra7GKZZ180H2igZ6O23bt3A4Dd/kOTfcjPz0d8fLzdLCS5lV0EGuh5k4G7uzveffddHDp0SO1ySFBHjhyx2/EzMMDCkoCAALmpqQnV1dUYPny4Dcvq239vVkiShKqqKuN+30RDdf36dXh5eaG+vt64uEmrLFpY0tzcjCtXriA8PFwTzy7Hx8ejsrISQM/NC39/f6Snp0OSJEiSxDvgNCTHjx9HeHi45sNsislA37yjXF5ejjNnztikoIFMnDgRsiwjMjISAPDGG28AAAIDA/HBBx+oWRrZufPnz+Phhx9Wu4whMRloJycnJCYmYuzYsdDpdJpapfXTTz+hsLDQ+H1FRQV27typYkVk72pqajB+/Hi1yxgSk2NoSZJs+2QGEZlLJ8vybWtTTQaaiOyL3UxbEdHAGGgigTDQRAJhoIkEwkATCeT/lMOEy78WhG0AAAAASUVORK5CYII=", "text/plain": [ "
" ] @@ -328,7 +330,9 @@ "\n", "gridDir = f'{config[\"Res\"]}_{config[\"NestedRegion\"]}'\n", "land_cover_file = f'{met_token1}.{constYr}0101.CN.{gridRes}.{config[\"NestedRegion\"]}.nc'\n", - "land_cover_pth = f'/home/ubuntu/ExtData/GEOS_{gridDir}/{met_token2}/{constYr}/01/{land_cover_file}'" + "land_cover_pth = f'/home/ubuntu/ExtData/GEOS_{gridDir}/{met_token2}/{constYr}/01/{land_cover_file}'\n", + "hemco_diag_file = f'HEMCO_sa_diagnostics.{gridRes}.20190101.nc'\n", + "hemco_diag_pth = f'/home/ubuntu/ExtData/HEMCO/CH4/v2023-04/HEMCO_SA_Output/{hemco_diag_file}'" ] }, { @@ -339,19 +343,24 @@ "source": [ "# Load land cover data\n", "lc = xr.load_dataset(land_cover_pth)\n", + "hd = xr.load_dataset(hemco_diag_pth)\n", + "hd[\"lon\"] = hd[\"lon\"] - 0.03125 # initially offset by 0.03125 degrees\n", "\n", "# Group together\n", "lc = (lc['FRLAKE'] + lc['FRLAND'] + lc['FRLANDIC']).drop('time').squeeze()\n", - " \n", + "hd = (hd[\"EmisCH4_Oil\"] + hd[\"EmisCH4_Gas\"]).drop(\"time\").squeeze()\n", + "\n", "# Subset the area of interest\n", "lc = lc.isel(lon=lc.lon>=lon_min, lat=lc.lat>=lat_min)\n", "lc = lc.isel(lon=lc.lon<=lon_max, lat=lc.lat<=lat_max)\n", + "hd = hd.isel(lon=hd.lon >= lon_min, lat=hd.lat >= lat_min)\n", + "hd = hd.isel(lon=hd.lon <= lon_max, lat=hd.lat <= lat_max)\n", "\n", - "# Set pixels over water to 0\n", + "# Set pixels over water to 0, unless there are offshore emissions\n", "if config['LandThreshold']:\n", - " # Where there is no land, replace with 0\n", - " land = lc.where(lc > config['LandThreshold'])\n", - " state_vector.values[land.isnull()] = 0" + " # Where there is neither land nor emissions, replace with 0\n", + " land = lc.where((lc > config['LandThreshold']) | (hd > config['OffshoreEmisThreshold']))\n", + " state_vector.values[land.isnull().values] = 0" ] }, { @@ -368,7 +377,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAADnCAYAAAApbXvLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQyElEQVR4nO3deZhVxZ3G8e9lF2hZFJCgCIwBGgRRkWCkI9oOShOFx2VEFMI4YnADXCIGCIMoRkEU2cKmDsGMCMQJsgyYgGJjFAkZoiC4EIw6YFhlCyJ0n/xRp6VPV/elt3PP6er38zw891Z1dfdPH17qnq0q4XkeIuKGKlEXICLlR4EWcYgCLeIQBVrEIQq0iEOqJftiRiLh7U9VJSJSbJthped51xTsTxro/cCi0EoSkdJKhzML69dHbhGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHVIu6ACmZtAdHW33zd8wPtHfvy7TGrDh6efB73upbvoVJLGiGFnGIAi3iEH3kjrH09b+z+kaMetrqG7jizUA7UWO9Neb87asD7SUX/Is15tObBlp9D2QOCrRr1dtrjWk+c6rVN/zWiwLt268fZI2R8qcZWsQhCrSIQxRoEYfoGDpGrtqRG2jPu2eBNea9lbdZff9z57RA+y831LfGTDyUFmg/dctfrTFL9o6y+haPDLb/NrKRNWZU5gdWX7dj91p9Ej7N0CIOUaBFHKJAizhEx9ARmZf1hNX35a4Vgfa1dc+1xkxfssvqO5zRLtDuMuxRa8zXf3o20D6r6mJrzI4tv7YLza4XaHYeOscacukFja2+Ork59s+S0GmGFnGIAi3iEH3kTpGDI7YF2q1WLLXGrGvcLdC+9fT51pjqOUOsvvH1/iPQnjz6Y2vMkKuDt4zW6jfeGjNiaK7VN5WDgfbw+e2sMbUb2Zet2mwaEGgftkZIGDRDizhEgRZxiAIt4hAdQ4cgjS+svs8zgsfDzy/9mTXmoZteDLSrjDnbGjP228esvs4fBG/rbDTxEWvMsQcbBNr9hv/AGvNuz/Osvq1/DN7qOaL7PmvMmdumWH2H37jK6pPwaYYWcYgCLeIQBVrEITqGDsFP77jT6ptQN3hb5RUstMbMvHZdoN3/X9+2xjw6r4HV139K+0B79n961pj9n+wOtFf3q22NGdr2HKvvwTeDz0/2bW7X/WG9NKtPoqEZWsQhCrSIQ/SROwT3zWlo9aWffX2gXaOKfYnok0WJQDt3wFBrzLITW6y+PZ0zAu0WIxLWmO+//KNA+/+/Z9dY47RbrL6x2cFF+1s2b2+NafO4/fskGpqhRRyiQIs4RIEWcYiOoUPwzTtjrL5DXY8G2mPu7WSN2TBrVaD92yZv2T/7fft4deAzvwi0d0+7yBoz+Q/Bxx5bTL/bGnNrQ/uRyhl/3RBoz/3eP6wxL1g9EhXN0CIOUaBFHKJAizhEx9AhWD/BXnWzT6srAu0F9aZbY5buez/Q7vmqfXtm7SP277so/Xig3WqTvQrn9tbBZYrO6zHPGvNUTXtFzy4Xzwy0f1jgXADA+Mn2skgPM9kuVEKnGVrEIQq0iEP0kTsEt736ktX3ZYN3Au1VT9xvjZnxafVA+8YBGdaYm7fbK4FecknPQDu3nX1JalL97ED7goR9m+fkYfZKK4caBX/fY73tlUhaVXnI6nv4+HOB9viquj00FTRDizhEgRZxiAIt4hAdQ6fIbQeDu1ls+f04a8zcVw8E2mM6HLfG3LHJvvVy2ofBlU2OfnXIGtP1gfWB9o4h1a0x50ywj8/rPx7c1O7ngwdbY375gr0x/YKqfaw+CZ9maBGHKNAiDlGgRRyS8Dx7hcg85ycS3qIUFlOZTN45xupr0r1+oP3nuvYSRF2HVbX6HhkYXC30nRtHWmMWJ4K3bE6obx9De0PsHTCOpQeXHLr/rt9ZY27p1Mfqa5Ib/HuVc7euQ5endNjgeV7ngv2aoUUcokCLOESXrSIypOkYq2/W2BmB9msdJlhjLvkox+q7LCf4cXbtXfZh1LE2awLtQU9cYY05bZS9mXuXxNhAe9oZn1ljEjUfsPraffFkoG1vCS9h0Awt4hAFWsQhCrSIQ3QMHSN3jg7eVmlveQfpvR62+p4dEWx/3N2+RNRp5cpAe1iDQdaYjOH271t1dnCJlCG0ssYsGDjM6ptV4Hrn4bX2z67bze6TstEMLeIQBVrEIQq0iEN066d8572f2UsXHe5+eqDdvrd9/fqlPr+3+hotejrQHndkkjXmozrDSlagfEe3fopUAgq0iEMUaBGH6BhakppG70B7ym/+Zo25ZfpGq2/+29uCHTkHrTFbql5YptoqMx1Di1QCCrSIQ3TrpyR1D8FVP7feao9Z8KTdx/27Ak1vrL3x3tayFFaJfZ7ka5qhRWLuIJAOzAaOAvcmGatAi8Rcrv/6DHAVyUOrQIvEXH1gOHA6cD9g7yx+ko6hpcw2P1JY76WB1sRu9uXRXuGU46SvgTpAG0zAi6IZWqQC+AmwE/g3oEaScQq0SAXQABjlv++QZJwCLVJBtCzGGB1DS0r0WqudM8rit5gZ+jTMpauiaIYWibG+mGvQ4/32UcyZ7qIo0CIxdg5QE3NzSQ3gDQpfPDKPAi0SUx4mxGc2a8acOXM45nl09zzaJnlCUsfQIjG0B3gYc/05Ozubli2Lc0pMM7RIrHjAh8AAzOWp+VDsMINmaJHY2A6MBT4FhgA3leJnKNAi5Ww7Zqa19xgp3A5gLrAE+CmQffw41aqVLpoKtEg5mANMwjw4sRp4BXgWqIq51HQN5gTXt8AizHXlD/3vrQXcgAn0GVDqMIMCLVIuEkAO8AjQ0O/Lf734VeBi7CelugGz/O8vDzopJlIO/t1/3Q/kLY94FfBL//06gmHuDryJWbSgPO+hU6BFymA/cBwTpKfy9a8AXs/N5byFCwPjN27ciOd5vOF5XO5fUy74pyz0kVuklPYDP/TffwBk5PtaLyCnipkv27Rpw2uvvUbr1q1Dr0kztEgp/Tnf+ymYs9sZQG3M8TTA+vXr2bp1a0rCDJqhRUrFA67EnOg6C3Ni6xXgQL4xx8tw+am0NEOLlNAOoCswAfg70BbzrHJemEcD6ynb5afS0gwtUgy7gJGYGfAtv+9FzJNQE/32AGAY5pnlqCjQIsUwEVhbSP8xzB1hLwBNUlpR4RRokWL4rIj+9phj5/ZlvNxUXnQMLVIMr3DyslRX/3UI5jbOqpFUVDjN0CIFHMYcJ9fx31fHzMRNMGe2V2NWErkrqgKTSBro3GRfFHHURuDBfO2zgK8wIc7bcm9pimsqrqSB1udxqYy6AfMwC/PVxZwQmw7sBlb6X+sYk2PmgpRZkUJ0xjw4URvIBFpgwjzI/1pcKdAiRaiH2fHRAx73+4ZFVk3xKNAiSSzh5F5Si4l/YOJen0joPODjIvpnY9bEzgRS83hF2SjQUuntB3oDd/vtXGAL8Dawz+8bGkFdpZH0LPfOVFUhEqGGQEfMrhSzMcfNYK47H8Lc2vn9aEorsaQz9D7g/1JUiEhUjgDv+++fydf/d/913Msvl+uqImFKGuiqwGOY/2ARV9UGrsfcGVaYzMzMFFZTNgkvyb825ycSXntgA+YYoyNwAeYJEw+on4ICRVIpB3OouRx4DnM8nSwjUUkkEhs8z7MuiZ8y0AsxJweyMbfEfYJZZ7gaMA64LoxqRSL0LWbiylORAp30I3etiy8m3fO4w/OY63m8glnN8DrgBDA8nFpFIuMB/52vXdEuA5W43saY5Urz9t35tFzLEYnWcczf73HAuVS8B5RK/Q/Qo/7rR+VUiEgcVPdfR2LWCatoShTo/Kft0z2PqVOnkh1WZSIROJHvfd6SQidOnChsaCyV6RChT58+ZHNy0TSA/wV+jD6KS8VUHfN3uA5whd+3adOm6AoqoTIFulmzZkzFfDzJAi4HHsDs7XMt5h5YkYqmBXAh5vlngBkzZkRXTAmV+STehcAq4CHMdbuFmLWXqmOeHVWopSL6EbAJaA7MnDkz4mqKr8yBbut5dPQ87vY8+noeN3oes4C/AB2AfuhOM6l4OmEWBsw7+bt58+boiimB0C6zJYBRmBn8mVOMFYmbKf5rB/918ODBUZVSIqFfN78eszfujrB/kUg56o2ZlHYDGRkZrF1b2DL78RN6oL/EnCTLpOjFykXipidwmf96dba5OPtuIsFW/09chR7oJ/O974m5T1Yk7qpgVv2Ek4sGLIuolpIINdDbOLniQ3P/tUYRY0XipgFmEtrot+dGV0qxhRroX/mvMzEX69eE+ctEQnAZZqGDdsAXEddSHKFthfNrYCvwLmY5VDAPdohUJF0x+z3nAmdgDhnj/CkzlBn6CDAVmMHJMItURM2A2/z3e4n/id1QZui8Db6ahfHDRVLsJqAp5lNnfM9vG6HM0I0xH1G2hfHDRVLsPOAG4ADxf6QylEAfAL72f3hRi7e8h+4gk4pjJ3Am8d9/OZRA52252Qvoj1kFIr/jwE+ABWH8cpEQ5BCvjd2LEkqg858F3IBZLXSh354PXOK/Xx7GLxcJwUrgB1EXUQyhXYd+DvPxpKbfHg2kY55eOQZ0919F4u5b4HWgb9SFFEMogW7reQzxPI57Ht94HosLGfMmcCUwJowCRMrRHzBnudtEXUgxpGSV0taY2+cKfmR5CPPMaU4qihAphVzMflf9oi6kmFK27HBN4HmgR76+p/3X8zEfyeO3nLlUdksx54SujLqQYkrpOuJVMcfWW4A/AgPzfW0hcB+wK5UFiSSxD5gA/Jz431CSJyWX1Yrare9S4EX//apEghmYXTnOxVzHfh44OwX1iRTmMczfx74x3AqnKLHZ6aMZ5n/gKszWnp8Dv6BiPIMq7lmFebjovqgLKaHYBDpPHcyjamCe1BqNWUX0v6IqSCqlRcA9QK2oCymh2AQ6/64cmz2PSZMmAfAPYA9mv6HH0YkzCV+Lo0fZkJbG7Xv3xnpz98LEJtAFDR06lAkF+n4DHIqiGKlU1qxZQ8eOHWnYsGHUpZRYbAMNZkudczDLF93u930VXTlSSSxfvpxevXpFXUapxDrQYBZKOIBZwgigUYS1iPs8YNmyZWRlZUVdSqnEOtBtPY/rPI8FK1ey0+97PdKKxHWfAUePHqVjx45Rl1IqsQ50nh49ejB79mzA3EYqEpa3gKysLBIxXns7mQoRaID+/fuTlpbGdGB11MWIs9ZAhT1+hlME+uDBgzRu3JjDhw+nqp4i1axZk/nz57MWc32wHebmE5HycgSzyWJmZmbUpZRa0kAfOHCA3bt306VLF7755ptU1VSkrKys706OecDVmBtO0v0/O4v4PpHieA+zGEdaWlrUpZRa0kA3bdoUgC1btrBu3bqUFHQqLTAPd3Ty20/5ry2BFyKoR9yxDWgbdRFllDTQ1apV4+abb6Z58+bs2bMHL0Z3zbwMzMvX3g68FFEt4oYvqfgPAyWShTSRSPwphbWISPHt8TzvmoKdSQMtIhVLhblsJSKnpkCLOESBFnGIAi3iEAVaxCH/BIiMeD1780HgAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAADnCAYAAAApbXvLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQyElEQVR4nO3deZhVxZ3G8e9lF2hZFJCgCIwBGgRRkWCkI9oOShOFx2VEFMI4YnADXCIGCIMoRkEU2cKmDsGMCMQJsgyYgGJjFAkZoiC4EIw6YFhlCyJ0n/xRp6VPV/elt3PP6er38zw891Z1dfdPH17qnq0q4XkeIuKGKlEXICLlR4EWcYgCLeIQBVrEIQq0iEOqJftiRiLh7U9VJSJSbJthped51xTsTxro/cCi0EoSkdJKhzML69dHbhGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHKNAiDlGgRRyiQIs4RIEWcYgCLeIQBVrEIQq0iEMUaBGHVIu6ACmZtAdHW33zd8wPtHfvy7TGrDh6efB73upbvoVJLGiGFnGIAi3iEH3kjrH09b+z+kaMetrqG7jizUA7UWO9Neb87asD7SUX/Is15tObBlp9D2QOCrRr1dtrjWk+c6rVN/zWiwLt268fZI2R8qcZWsQhCrSIQxRoEYfoGDpGrtqRG2jPu2eBNea9lbdZff9z57RA+y831LfGTDyUFmg/dctfrTFL9o6y+haPDLb/NrKRNWZU5gdWX7dj91p9Ej7N0CIOUaBFHKJAizhEx9ARmZf1hNX35a4Vgfa1dc+1xkxfssvqO5zRLtDuMuxRa8zXf3o20D6r6mJrzI4tv7YLza4XaHYeOscacukFja2+Ork59s+S0GmGFnGIAi3iEH3kTpGDI7YF2q1WLLXGrGvcLdC+9fT51pjqOUOsvvH1/iPQnjz6Y2vMkKuDt4zW6jfeGjNiaK7VN5WDgfbw+e2sMbUb2Zet2mwaEGgftkZIGDRDizhEgRZxiAIt4hAdQ4cgjS+svs8zgsfDzy/9mTXmoZteDLSrjDnbGjP228esvs4fBG/rbDTxEWvMsQcbBNr9hv/AGvNuz/Osvq1/DN7qOaL7PmvMmdumWH2H37jK6pPwaYYWcYgCLeIQBVrEITqGDsFP77jT6ptQN3hb5RUstMbMvHZdoN3/X9+2xjw6r4HV139K+0B79n961pj9n+wOtFf3q22NGdr2HKvvwTeDz0/2bW7X/WG9NKtPoqEZWsQhCrSIQ/SROwT3zWlo9aWffX2gXaOKfYnok0WJQDt3wFBrzLITW6y+PZ0zAu0WIxLWmO+//KNA+/+/Z9dY47RbrL6x2cFF+1s2b2+NafO4/fskGpqhRRyiQIs4RIEWcYiOoUPwzTtjrL5DXY8G2mPu7WSN2TBrVaD92yZv2T/7fft4deAzvwi0d0+7yBoz+Q/Bxx5bTL/bGnNrQ/uRyhl/3RBoz/3eP6wxL1g9EhXN0CIOUaBFHKJAizhEx9AhWD/BXnWzT6srAu0F9aZbY5buez/Q7vmqfXtm7SP277so/Xig3WqTvQrn9tbBZYrO6zHPGvNUTXtFzy4Xzwy0f1jgXADA+Mn2skgPM9kuVEKnGVrEIQq0iEP0kTsEt736ktX3ZYN3Au1VT9xvjZnxafVA+8YBGdaYm7fbK4FecknPQDu3nX1JalL97ED7goR9m+fkYfZKK4caBX/fY73tlUhaVXnI6nv4+HOB9viquj00FTRDizhEgRZxiAIt4hAdQ6fIbQeDu1ls+f04a8zcVw8E2mM6HLfG3LHJvvVy2ofBlU2OfnXIGtP1gfWB9o4h1a0x50ywj8/rPx7c1O7ngwdbY375gr0x/YKqfaw+CZ9maBGHKNAiDlGgRRyS8Dx7hcg85ycS3qIUFlOZTN45xupr0r1+oP3nuvYSRF2HVbX6HhkYXC30nRtHWmMWJ4K3bE6obx9De0PsHTCOpQeXHLr/rt9ZY27p1Mfqa5Ib/HuVc7euQ5endNjgeV7ngv2aoUUcokCLOESXrSIypOkYq2/W2BmB9msdJlhjLvkox+q7LCf4cXbtXfZh1LE2awLtQU9cYY05bZS9mXuXxNhAe9oZn1ljEjUfsPraffFkoG1vCS9h0Awt4hAFWsQhCrSIQ3QMHSN3jg7eVmlveQfpvR62+p4dEWx/3N2+RNRp5cpAe1iDQdaYjOH271t1dnCJlCG0ssYsGDjM6ptV4Hrn4bX2z67bze6TstEMLeIQBVrEIQq0iEN066d8572f2UsXHe5+eqDdvrd9/fqlPr+3+hotejrQHndkkjXmozrDSlagfEe3fopUAgq0iEMUaBGH6BhakppG70B7ym/+Zo25ZfpGq2/+29uCHTkHrTFbql5YptoqMx1Di1QCCrSIQ3TrpyR1D8FVP7feao9Z8KTdx/27Ak1vrL3x3tayFFaJfZ7ka5qhRWLuIJAOzAaOAvcmGatAi8Rcrv/6DHAVyUOrQIvEXH1gOHA6cD9g7yx+ko6hpcw2P1JY76WB1sRu9uXRXuGU46SvgTpAG0zAi6IZWqQC+AmwE/g3oEaScQq0SAXQABjlv++QZJwCLVJBtCzGGB1DS0r0WqudM8rit5gZ+jTMpauiaIYWibG+mGvQ4/32UcyZ7qIo0CIxdg5QE3NzSQ3gDQpfPDKPAi0SUx4mxGc2a8acOXM45nl09zzaJnlCUsfQIjG0B3gYc/05Ozubli2Lc0pMM7RIrHjAh8AAzOWp+VDsMINmaJHY2A6MBT4FhgA3leJnKNAi5Ww7Zqa19xgp3A5gLrAE+CmQffw41aqVLpoKtEg5mANMwjw4sRp4BXgWqIq51HQN5gTXt8AizHXlD/3vrQXcgAn0GVDqMIMCLVIuEkAO8AjQ0O/Lf734VeBi7CelugGz/O8vDzopJlIO/t1/3Q/kLY94FfBL//06gmHuDryJWbSgPO+hU6BFymA/cBwTpKfy9a8AXs/N5byFCwPjN27ciOd5vOF5XO5fUy74pyz0kVuklPYDP/TffwBk5PtaLyCnipkv27Rpw2uvvUbr1q1Dr0kztEgp/Tnf+ymYs9sZQG3M8TTA+vXr2bp1a0rCDJqhRUrFA67EnOg6C3Ni6xXgQL4xx8tw+am0NEOLlNAOoCswAfg70BbzrHJemEcD6ynb5afS0gwtUgy7gJGYGfAtv+9FzJNQE/32AGAY5pnlqCjQIsUwEVhbSP8xzB1hLwBNUlpR4RRokWL4rIj+9phj5/ZlvNxUXnQMLVIMr3DyslRX/3UI5jbOqpFUVDjN0CIFHMYcJ9fx31fHzMRNMGe2V2NWErkrqgKTSBro3GRfFHHURuDBfO2zgK8wIc7bcm9pimsqrqSB1udxqYy6AfMwC/PVxZwQmw7sBlb6X+sYk2PmgpRZkUJ0xjw4URvIBFpgwjzI/1pcKdAiRaiH2fHRAx73+4ZFVk3xKNAiSSzh5F5Si4l/YOJen0joPODjIvpnY9bEzgRS83hF2SjQUuntB3oDd/vtXGAL8Dawz+8bGkFdpZH0LPfOVFUhEqGGQEfMrhSzMcfNYK47H8Lc2vn9aEorsaQz9D7g/1JUiEhUjgDv+++fydf/d/913Msvl+uqImFKGuiqwGOY/2ARV9UGrsfcGVaYzMzMFFZTNgkvyb825ycSXntgA+YYoyNwAeYJEw+on4ICRVIpB3OouRx4DnM8nSwjUUkkEhs8z7MuiZ8y0AsxJweyMbfEfYJZZ7gaMA64LoxqRSL0LWbiylORAp30I3etiy8m3fO4w/OY63m8glnN8DrgBDA8nFpFIuMB/52vXdEuA5W43saY5Urz9t35tFzLEYnWcczf73HAuVS8B5RK/Q/Qo/7rR+VUiEgcVPdfR2LWCatoShTo/Kft0z2PqVOnkh1WZSIROJHvfd6SQidOnChsaCyV6RChT58+ZHNy0TSA/wV+jD6KS8VUHfN3uA5whd+3adOm6AoqoTIFulmzZkzFfDzJAi4HHsDs7XMt5h5YkYqmBXAh5vlngBkzZkRXTAmV+STehcAq4CHMdbuFmLWXqmOeHVWopSL6EbAJaA7MnDkz4mqKr8yBbut5dPQ87vY8+noeN3oes4C/AB2AfuhOM6l4OmEWBsw7+bt58+boiimB0C6zJYBRmBn8mVOMFYmbKf5rB/918ODBUZVSIqFfN78eszfujrB/kUg56o2ZlHYDGRkZrF1b2DL78RN6oL/EnCTLpOjFykXipidwmf96dba5OPtuIsFW/09chR7oJ/O974m5T1Yk7qpgVv2Ek4sGLIuolpIINdDbOLniQ3P/tUYRY0XipgFmEtrot+dGV0qxhRroX/mvMzEX69eE+ctEQnAZZqGDdsAXEddSHKFthfNrYCvwLmY5VDAPdohUJF0x+z3nAmdgDhnj/CkzlBn6CDAVmMHJMItURM2A2/z3e4n/id1QZui8Db6ahfHDRVLsJqAp5lNnfM9vG6HM0I0xH1G2hfHDRVLsPOAG4ADxf6QylEAfAL72f3hRi7e8h+4gk4pjJ3Am8d9/OZRA52252Qvoj1kFIr/jwE+ABWH8cpEQ5BCvjd2LEkqg858F3IBZLXSh354PXOK/Xx7GLxcJwUrgB1EXUQyhXYd+DvPxpKbfHg2kY55eOQZ0919F4u5b4HWgb9SFFEMogW7reQzxPI57Ht94HosLGfMmcCUwJowCRMrRHzBnudtEXUgxpGSV0taY2+cKfmR5CPPMaU4qihAphVzMflf9oi6kmFK27HBN4HmgR76+p/3X8zEfyeO3nLlUdksx54SujLqQYkrpOuJVMcfWW4A/AgPzfW0hcB+wK5UFiSSxD5gA/Jz431CSJyWX1Yrare9S4EX//apEghmYXTnOxVzHfh44OwX1iRTmMczfx74x3AqnKLHZ6aMZ5n/gKszWnp8Dv6BiPIMq7lmFebjovqgLKaHYBDpPHcyjamCe1BqNWUX0v6IqSCqlRcA9QK2oCymh2AQ6/64cmz2PSZMmAfAPYA9mv6HH0YkzCV+Lo0fZkJbG7Xv3xnpz98LEJtAFDR06lAkF+n4DHIqiGKlU1qxZQ8eOHWnYsGHUpZRYbAMNZkudczDLF93u930VXTlSSSxfvpxevXpFXUapxDrQYBZKOIBZwgigUYS1iPs8YNmyZWRlZUVdSqnEOtBtPY/rPI8FK1ey0+97PdKKxHWfAUePHqVjx45Rl1IqsQ50nh49ejB79mzA3EYqEpa3gKysLBIxXns7mQoRaID+/fuTlpbGdGB11MWIs9ZAhT1+hlME+uDBgzRu3JjDhw+nqp4i1axZk/nz57MWc32wHebmE5HycgSzyWJmZmbUpZRa0kAfOHCA3bt306VLF7755ptU1VSkrKys706OecDVmBtO0v0/O4v4PpHieA+zGEdaWlrUpZRa0kA3bdoUgC1btrBu3bqUFHQqLTAPd3Ty20/5ry2BFyKoR9yxDWgbdRFllDTQ1apV4+abb6Z58+bs2bMHL0Z3zbwMzMvX3g68FFEt4oYvqfgPAyWShTSRSPwphbWISPHt8TzvmoKdSQMtIhVLhblsJSKnpkCLOESBFnGIAi3iEAVaxCH/BIiMeD1780HgAAAAAElFTkSuQmCC", "text/plain": [ "
" ] @@ -390,10 +399,10 @@ " state_vector[r,c] = count\n", " count += 1\n", "\n", - "# Now set pixels over water to NaN\n", + "# Now set pixels over water to NaN, unless there are offshore emissions\n", "if config['LandThreshold']:\n", - " # Where there is no land, replace with NaN\n", - " state_vector = state_vector.where(lc > config['LandThreshold'])\n", + " # Where there is no land or offshore emissions, replace with NaN\n", + " state_vector = state_vector.where((lc > config['LandThreshold']) | (hd > config['OffshoreEmisThreshold']))\n", " \n", "# Plot\n", "f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))\n", @@ -420,7 +429,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAADnCAYAAAApbXvLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVyElEQVR4nO3deVxVdfoH8M/lsiiKKKCBqIColCiQOTq5pOYalVnua5mmNTVKZWlOP/ctc1fcy3FNHckldwtFwR3TEHfcERXR2GXxnt8fGDOHLyDbvefer5/3X32f+8B9Xq+Zj+cczqZTFAVEJAcrrQcgorLDQBNJhIEmkggDTSQRBppIItaFfVjbpZaSnvXYVLNIwzbbVesRzI6Dq07rEaQSffWPPYqidMpbLzTQ6VmP8f6r3Yw3laS8HozUegSz03pUOa1HkMqL3V1d8qtzl5tIIgw0kUQYaCKJMNBEEmGgiSTCQBNJpNDTVvRsPEVF5oRbaCKJMNBEEmGgiSTCQBNJhIEmkggDTSQRBppIIgw0kUQYaCKJMNBEEmGgiSTCQBNJhIEmkggDTSQR3j5JJnFgmvg4aD4JtOxxC00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATScRa6wHo+VVv/XitR5AOt9BEEmGgiSTCQBNJpNBj6ESDA3antzLVLJapwlGh9Enq3432dX03egi1LXZBqnWc/3ahZ8D9t9U/0yRI6DGmIZ7fmfT7nlfcQhNJhIEmkghPWxnBojLaDR/Y7z9CbeTsr4XahgSDam1wniP0fH25p2r9eFhNoWeP7QWhFvlSOdV6QlaU0DP62xlCbWDoV+qCp9BCRsAtNJFEGGgiiTDQRBLRKYpS4If2Dt5KvUZTTDjO8yO/Y+pWQdGq9daTT4SekWGPxF+2x1G1bDx8udDyqr+tat3FYZzQ8/57DYTalB4jVOuT4fWEnhof3RJqZ5Irq9Zrsq4KPVRyViGLIhVFaSzUtRiGiIyDgSaSCANNJBGeh9bI0aSBQu1G9nDVesHoJKGn6nfNhFpcuWWq9ZY2t4WecwdbqNahOvHy0NstxJ+7NdJTtR5QuYfQc/77XUKtSVCeWyOz2go9VPa4hSaSCANNJBHucpvImT43VGvfkWFCz649sar1hcNVhZ7Rl6cKNZfA+ar1x+VThZ53r9VVrTtvF3f5bSMyhNqqA+rd8D2/fCX0zLv2g1BbeVunWrdy5C63KXALTSQRBppIIgw0kUR4DG0Enr7xQm334ZfU6wk/Cz2GA/9SrVdWXSf0fOg3V6h5z/NRraN3eQo99Y84q9bHu78m9KwZV0OoHQhsqlp3mSVe5ln1j1FCLd0xUaiR8XELTSQRBppIIgw0kUR4DG0Eo8aKjw6afdZJtfYv31vomXBIfa7WULun0PNldBWh9v0h9WOCwj4Vv3/+QPVlpT0jFgg9He+Kv7uH22rVutO+vULP9D5rhBq2jhVrZHTcQhNJhIEmkgh3uY1g1K/i0zPD1qufsvntZ/WFnuRz6l3sldUHCz2Ha1wUv9BR/dRPW6s6Qkv02r6qtWHAcKFnR/Z5obbgQLpqPVu8OhQtuHttNriFJpIIA00kEQaaSCI8hjaCwFsbhNrsV5ao1tbKNLHn3DDV+oP654Serb1/FGpdqrup1pE99gk9EaM/Va0vfiKefnL7YIRQi+xup1rHd/UVerD//8QaaYJbaCKJMNBEEmGgiSTCY2gjSI5qI9ReHtFOtZ7U+Huh50Woj6FvGsSXuw8NOSjUNjouVK2Db4rfPybP99nPEN+Q2ehJE6E2y6WTat01c5DQM6hhqFBbFtVBqJHxcQtNJBEGmkgi3OU2gl9urxdq2R9/rlpP/XGj0LPxg5uqdY/j5YSeB+3ThVr3PHvPq8qJL4WfM6CXar25pfj0TkPqMKE2cJn6ktHhnx0Xevx14p1jXVbPUa23BPQXeqjscQtNJBEGmkgiDDSRRHgMbSJJjRap1jcbPhB6Xl54XbX+9Z8Gocd332Sh1uuNKNV6aqL4QrvBv6iPhYOdxZfCp99NFmrNjqvfinFnmI3Qs3mQ+DeDz6+MUxe6xQg9VPa4hSaSCANNJBEGmkgiPIY2kYhzj1XrEcGThJ6glurbJX++XUnoeaF1ZaFW7UP1Oe1/LRa//8hp9QvfwwdmCj3643FC7efe6ieRTh8j/vJL9r2Emgd0eSqfiENRmeMWmkgiDDSRRHSKohT4ob2Dt1Kv0RQTjvN8G+6g3nVt+I3Ysy5U3OWd3TBJtZ54Ubysc+co9S5w+AHxf/cwH/El9DWnqO/cKh8rnkp7UzdRqAU7qw8XmiXk87RSKjGrkEWRiqI0FupaDENExsFAE0mEgSaSCE9bmZG5yepLKPsHjxF6Frq3FmqGDrtV61kbfYSeHqPVa6ttQ4WegI5dhVpQlY9U65YjhRaMXNRUqC2crH7pfevDQUJP5kyeyipr3EITSYSBJpIIA00kER5Dm5H+7uIxc16TYzsLtVniSyoFT67VU61n4oDQM6KdeEy78lf1/0Xco4QWXGy2X6jZd8tWrb8S7+gExEN2KiVuoYkkwkATSYSBJpIIj6E1UpTjZVOb0a+nUEtfd0m13qKLFXrmr70h1HrHnVatl/cUH/U7OVsoUSlxC00kEQaaSCLc5TYRc9zFLoqUPurTXZ69zws9MxeJt0ZufDvPE0vqiG/lwE+rSjXb8yomJbHAz7iFJjJzf2ZmwCpkEaZdPIW07Cy8e2R3gb0MNJGZMyDnYRSjzx6D1+41sNLlfV7bfzHQRGbOybYcZjRshso2tpji+3dsa/ZGgb08hjYCSz1eLrGa4u2a08+oj73fHPWq+HNdeftkUSVkPkYlG1v4OTrD2VZ8K+lfuIUmsgBf1PXHzbQUNN0fAjsrfYF9DDSRBXC2K4f5AS0BALablxTYx0ATWQifipWf2cNj6FJ67o6XS2hH73piDb+p1r2qi2/goBzro37HiPBfYG9jg7SsrAL7GGgiM9Z53Q84FRcLR7ucP4SlZWVhVIvXMS08NN9+7nITmbFajlVgp7dGYsZj2On1OD4kCJ81bVFgf6Fb6Dpu5bH124ZlPqQl27TiXa1HoOeEoiiws7aGi+sLGD9+PAYNGvTfD2dOyPdnuMtNZIbiU1MwbOdmPHqcjkNHD8PLy6tIP8ddbiIzoigKzt6LQ7cNK+HvWh3b+gwqcpgBbqGJzEbMwwcY/dtOXHoQj6+at0Efv0bF/h0MNFEZi3n4AIoC1HF2KVJ/bFIilkUexc/nozCsaQsciLkMa+uSRZOBJioDC49H4LvwUKx4txf2xVzC6jORWPxWN+itrJCWlYm3fXxhq9cjIzsb68/+jvVRpxF1Pw4AUM7aGr0avIzQDz6Bi32FEocZYKCJyoROp8MTRUHQrq1wsa8AAPh4+6bczzecPY0m7rUw5+hB1c+18vTGmvf6QFfILZHFwT+KEZWBoY1z7iZ7mJ6GSwk5L+rrVOdFzO70DgDg8K3rqjC3q10XJ4YEYW3XvmUWZoBb6GLrNnCzas3z0s+3h2lpcLCzg41ej7lvdMHwXVsAAIc+/AzNl89DSEgIsHtrbv/p06fh7+9vtHm4hSYqoYdpafBbNANecyYj22BAG886uZ+1XhEMKysrdO/eHT4+Prh48SIURTFqmAEGmqjETty5lfvfMyL2I+bRA7Tx9Ia9jQ2eKDmPDTpx4gQuXLiAevXEm1OMgbvcRCWgKAo6eNdDU/dacHWohAXHI7D2j1N49Dg9tycrK6tUf7EuCQaaqJhikxLRYdUS9GoQgLjkJLTzrouz9+IQ8ygBADClbSDefamhycMMMNBERXIvJRlf7NkGvU6H0GtXAABLIo/CTm+NyQdz7use3KgpRrZ4HeVtbDSbk4EmKoIph35D2PUYoZ7xJBt1nFzwU7d+cHOopMFkagw0URHEPHyQb73hC27Y3mcQas2aaOKJ8se/chMVwfa+g9Ha0xsA0LymJwDgq+atsavfR9BbmU+MuIUmyiM5IwP7r11BBVtbJGdkwEavh7+rG9wcKqGDtw/2xlyEh2MVDP/7a1qPKig00AbFYKo5iMxGZNwt/GNHSO7arWIlxKUkwcOxCirY2gIA9g/8h1bjFarQQFvpzGdXguRjrk/5bO1ZByE938eEsH2oZGuHhW91xawjB3E/NRk7Lp1HSI/34TV7ktZj5ouJJcpH0xoeWNu1L+xtbNF06VzUruKEHZfO49MmzdG0pofW4xWIgSYqQOVy5bHwra4wAPi/0JxXuI5s8bq2Qz0DA01UiM3no2Cnz3mX1L4BQwt9las5YKDpuacoCs7H38u3Hnw8AokZj9Gxjg9eqvqCBtMVT6GBVp6+aJpIZg/T09B+1RIM3LweAGBQFETfv4uwGzFISEsFAHzdvI2WIxZZoX/ljrt/11RzEGnG2b4CAlyrY9/VSwg+HoGph3KuzXat6ICkzAzUcXKBj0s1jacsmkIDnZD4EJHRv+MV35dNNQ+RyaVkZuD03TsAkBtmALibkgwAmBg8H+69zPMUW16F7nJb6/UYEzwZKU93O4hkVMHGFj19A1Dx6UUjebVt29bEE5WcTlEKPk72q+erNKjrixNnI9G1/Tvw92mIRi/5IyMzA4qioHKlyqab1EzxmWIlZ44XljwxGBCbnIhtF6IxPXw/DFBQWEa0otPpIhVFaZy3/oxruXWYMnwsDkZGIOxEOKb/OAcXr11GekY6bKyt8d0XE/Feu87GmpnI5PRWVnCtWKnA17Wau0K30I0bN1ZOnjyZu76x7xzuJ8Rj6g+zsPnXbQCA63ujjT6kpeFWu2jMcQutKAqWnzqG8Qf2AgCsdDo8MZjfPQ0FbaGLfR66mnNVzP56KnoHdgcAXL5xpQzGIzIPWQYDxh/Yi5kdO8OrshMMZri7XZgSX1gyZfhYAMD5q5fKbBgirdk8vbf5yz3b4O3krPE0xVes+6E92tdXrRcsWIDQLXvRuU1gmQ5FpJXs/9m9dq2Y80ih7OxsTR74VxKluvSzS5cuCDsRjv3HD+XWtoftRrvBb3NXnCySjV6Pgx9+ioq2tmjvnfMs7bNnz2o8VdGV6p8dd3d3LB0/H0PHD0OlCpWQmp6Kewn3AQDtP3oHZ34+AseK2j84jag4aldxxivVa+Le0wtLFi9ejMWLF2s8VdGU+uaMV+oHIGL1r/jmoy+waMwcbFuwAa3/1hK2NjYY8M0QJKYklcWcRCb1ulcdRN2Pg0flKliyZInW4xRZqQ8M/jqurvdmQG7Nr14DKIqCcQunoGtQX2yZtx4Vn75ik8gSvFK9Bsbu34MN3fqj56bViI6Ohq+vr9ZjPZPRbp/U6XQY/+m/0Nj3ZUz/cbaxvobIKGYeDgMABLi5AwA+/vhjLccpMqPfD929w7s4fPoYYu/fMfZXEZWZri/5QQfgXmoyWrZsifDwcK1HKhKjB/rW3VhcuXkVzfu1x7XYG8b+OqIy0flFX7zm6Y3XfgzGO+UrAwDOfvo1YkeMReyIsdoOVwijB3ri4u9y/7vNwEBkZGYa+yuJSs1Kp8P8N3Iu4bV++giirRfM//SVUQN95WYMEhIfAgA8qtcEANgVcIsakblxsrdHZx9fnLpzGwCwLPKoxhM9m1Evf5m3Nufc3YpJi9GqcXPEP8r//UBE5uo1j9rYc+UiGlZzRZQFPMHHaIH+cfNqnIu5gDMhh+Ho4AgAeMHZMh7jQvSXFh5e+HrfdhgUBS72FZCRnQ07M74M1Ci73ClpqZi9KhgrJi7KDTORJapRqTI+fLkJAOBBWiqu/flQ44kKZ5R/alLSUmBrbYMaru7G+PVEJtXHrxHcKzliWeRRmPdTuY20hX7BuRoMigFXboovyCayNPWcq6JngwD8+Tgd3k4uWo9TKKNsoROTE/Eo6U/odFZQFAW6fN42cOTMcRw8GYGRgz43xghkZszx6STFcSc5CVUrVIS1Gb0LOj9Gma58OXsAQLvBb6PHlwOQlZ2l+jwrOwu9vxqIdTs3GuPricrcE4MB1hbwNlajTPi/55pPnD2FuoEBWL9rEwBgzfYNaNClKQAg9Icdxvh6ojK3/dI5NK/lpfUYz2S0f3IWj5kDa7017GztAACjZo+FZwdffDtvAjIyM9C2aStkZGYY6+uJykxmdjZ2XDqHfv6vaD3KMxXrqZ8ltXfpVnQc2iXfz/q+1ROTh40p9XeYu+f9SaCWfAy97UI01p/9Heu69cutuc8Yr+FEZfjUz5Lw8aqLC9tP4VX/Jqr6qMFfYu32DXjy5IkpxiAqNoOiYMHxcLwf8DetRykSkx3ll7O1w5ppy9GpRfvc2rTlMwEA3m/44ZvZY83yDQX0fNt8Pgq2ej06PH2+mLkz6Z/t9Ho9Fo+Zg+t7o/H7pggMeu/93M9+2rUJQ8cPx/2EeFOORFSghLRUTArbh/FtOuV76tUcmeSi1LyP/wUADwDLuzbDcvwbABC+eh+Cf1qC9kM6w8vdE4+SHmHN1OWo6VbDFCMSCb79bRfeq++Hzut+0HqUIjObE2s1Xd0x7fMJOLzmV5y+8Adu3LmFUXPGYut+ntoi09tz5SKi4+9hRLPWWo9SLGYT6L9UKF8BDerkbNEjfj+K0XPG4W+9WmF5yEqNJ6PnyU9Rp/D5q6+hvI2N1qMUi9ncB/a/u+VRl6Mxd+5cBAUFITU9DanpaZi0ZDpu3b2Ncf8YbTHHM2SZnCd9g2NLZ2PDiaNwcnLSepxiMbst9F+GDx+OuaO+U9VWbl2HpNRkjSai50VYWBj8/PwsLsyAGQcaAN55/S3UcqsJz+q1MKTbQABAXLz5PzWCLNvOnTvx5ptvaj1GiZh1oAFg6bh5+DM5ETsO7gYAVHOqqvFEJDNFUbBjxw4EBlrmCxjNOtAe7euj45Au2LBpI2LvxwEAdh3ap/FUJLNrjx4iPT0dfn5+Wo9SImYd6L906NABy5YtAwC8WNsyrtghyxR67TICAwMt9g+vFhFoAOjfvz8cHBwwd81C7DsSqvU4JKnfrl622ONn4BmBTkpKQrVq1ZCSkmKqeQpkZ2eH9evX4+DJCHw09p/w6tgAN+7c1HoskkhqZiZOxcWibdu2Wo9SYoUGOjExEfHx8WjSpAkeP35sqpkKFBgYiP0rdgLI+eNFqw/ewPJN/4ZnB194dvDFnafH2UQlceTWdQS4VYeDg4PWo5RYoYF2c3MDAJw/fx7Hjh0zyUDP4uXuget7o9GofgAAYNLS7wEAtWt4Yel/Vmg4GVm6Swnx8K3qqvUYpVJooK2trdGzZ0/UqlULDx48MKvbG3+esxYbZ67KXV+9fQ3/3rpWw4nI0t1K+hM1HStrPUapFPrEEp1OV/rHlRCRMTxQFKVT3mKhgSYiy2Ixp62I6NkYaCKJMNBEEmGgiSTCQBNJ5P8Bed3ZmLfeKSQAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAADnCAYAAAApbXvLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVyElEQVR4nO3deVxVdfoH8M/lsiiKKKCBqIColCiQOTq5pOYalVnua5mmNTVKZWlOP/ctc1fcy3FNHckldwtFwR3TEHfcERXR2GXxnt8fGDOHLyDbvefer5/3X32f+8B9Xq+Zj+cczqZTFAVEJAcrrQcgorLDQBNJhIEmkggDTSQRBppIItaFfVjbpZaSnvXYVLNIwzbbVesRzI6Dq07rEaQSffWPPYqidMpbLzTQ6VmP8f6r3Yw3laS8HozUegSz03pUOa1HkMqL3V1d8qtzl5tIIgw0kUQYaCKJMNBEEmGgiSTCQBNJpNDTVvRsPEVF5oRbaCKJMNBEEmGgiSTCQBNJhIEmkggDTSQRBppIIgw0kUQYaCKJMNBEEmGgiSTCQBNJhIEmkggDTSQR3j5JJnFgmvg4aD4JtOxxC00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATSYSBJpIIA00kEQaaSCIMNJFEGGgiiTDQRBJhoIkkwkATScRa6wHo+VVv/XitR5AOt9BEEmGgiSTCQBNJpNBj6ESDA3antzLVLJapwlGh9Enq3432dX03egi1LXZBqnWc/3ahZ8D9t9U/0yRI6DGmIZ7fmfT7nlfcQhNJhIEmkghPWxnBojLaDR/Y7z9CbeTsr4XahgSDam1wniP0fH25p2r9eFhNoWeP7QWhFvlSOdV6QlaU0DP62xlCbWDoV+qCp9BCRsAtNJFEGGgiiTDQRBLRKYpS4If2Dt5KvUZTTDjO8yO/Y+pWQdGq9daTT4SekWGPxF+2x1G1bDx8udDyqr+tat3FYZzQ8/57DYTalB4jVOuT4fWEnhof3RJqZ5Irq9Zrsq4KPVRyViGLIhVFaSzUtRiGiIyDgSaSCANNJBGeh9bI0aSBQu1G9nDVesHoJKGn6nfNhFpcuWWq9ZY2t4WecwdbqNahOvHy0NstxJ+7NdJTtR5QuYfQc/77XUKtSVCeWyOz2go9VPa4hSaSCANNJBHucpvImT43VGvfkWFCz649sar1hcNVhZ7Rl6cKNZfA+ar1x+VThZ53r9VVrTtvF3f5bSMyhNqqA+rd8D2/fCX0zLv2g1BbeVunWrdy5C63KXALTSQRBppIIgw0kUR4DG0Enr7xQm334ZfU6wk/Cz2GA/9SrVdWXSf0fOg3V6h5z/NRraN3eQo99Y84q9bHu78m9KwZV0OoHQhsqlp3mSVe5ln1j1FCLd0xUaiR8XELTSQRBppIIgw0kUR4DG0Eo8aKjw6afdZJtfYv31vomXBIfa7WULun0PNldBWh9v0h9WOCwj4Vv3/+QPVlpT0jFgg9He+Kv7uH22rVutO+vULP9D5rhBq2jhVrZHTcQhNJhIEmkgh3uY1g1K/i0zPD1qufsvntZ/WFnuRz6l3sldUHCz2Ha1wUv9BR/dRPW6s6Qkv02r6qtWHAcKFnR/Z5obbgQLpqPVu8OhQtuHttNriFJpIIA00kEQaaSCI8hjaCwFsbhNrsV5ao1tbKNLHn3DDV+oP654Serb1/FGpdqrup1pE99gk9EaM/Va0vfiKefnL7YIRQi+xup1rHd/UVerD//8QaaYJbaCKJMNBEEmGgiSTCY2gjSI5qI9ReHtFOtZ7U+Huh50Woj6FvGsSXuw8NOSjUNjouVK2Db4rfPybP99nPEN+Q2ehJE6E2y6WTat01c5DQM6hhqFBbFtVBqJHxcQtNJBEGmkgi3OU2gl9urxdq2R9/rlpP/XGj0LPxg5uqdY/j5YSeB+3ThVr3PHvPq8qJL4WfM6CXar25pfj0TkPqMKE2cJn6ktHhnx0Xevx14p1jXVbPUa23BPQXeqjscQtNJBEGmkgiDDSRRHgMbSJJjRap1jcbPhB6Xl54XbX+9Z8Gocd332Sh1uuNKNV6aqL4QrvBv6iPhYOdxZfCp99NFmrNjqvfinFnmI3Qs3mQ+DeDz6+MUxe6xQg9VPa4hSaSCANNJBEGmkgiPIY2kYhzj1XrEcGThJ6glurbJX++XUnoeaF1ZaFW7UP1Oe1/LRa//8hp9QvfwwdmCj3643FC7efe6ieRTh8j/vJL9r2Emgd0eSqfiENRmeMWmkgiDDSRRHSKohT4ob2Dt1Kv0RQTjvN8G+6g3nVt+I3Ysy5U3OWd3TBJtZ54Ubysc+co9S5w+AHxf/cwH/El9DWnqO/cKh8rnkp7UzdRqAU7qw8XmiXk87RSKjGrkEWRiqI0FupaDENExsFAE0mEgSaSCE9bmZG5yepLKPsHjxF6Frq3FmqGDrtV61kbfYSeHqPVa6ttQ4WegI5dhVpQlY9U65YjhRaMXNRUqC2crH7pfevDQUJP5kyeyipr3EITSYSBJpIIA00kER5Dm5H+7uIxc16TYzsLtVniSyoFT67VU61n4oDQM6KdeEy78lf1/0Xco4QWXGy2X6jZd8tWrb8S7+gExEN2KiVuoYkkwkATSYSBJpIIj6E1UpTjZVOb0a+nUEtfd0m13qKLFXrmr70h1HrHnVatl/cUH/U7OVsoUSlxC00kEQaaSCLc5TYRc9zFLoqUPurTXZ69zws9MxeJt0ZufDvPE0vqiG/lwE+rSjXb8yomJbHAz7iFJjJzf2ZmwCpkEaZdPIW07Cy8e2R3gb0MNJGZMyDnYRSjzx6D1+41sNLlfV7bfzHQRGbOybYcZjRshso2tpji+3dsa/ZGgb08hjYCSz1eLrGa4u2a08+oj73fHPWq+HNdeftkUSVkPkYlG1v4OTrD2VZ8K+lfuIUmsgBf1PXHzbQUNN0fAjsrfYF9DDSRBXC2K4f5AS0BALablxTYx0ATWQifipWf2cNj6FJ67o6XS2hH73piDb+p1r2qi2/goBzro37HiPBfYG9jg7SsrAL7GGgiM9Z53Q84FRcLR7ucP4SlZWVhVIvXMS08NN9+7nITmbFajlVgp7dGYsZj2On1OD4kCJ81bVFgf6Fb6Dpu5bH124ZlPqQl27TiXa1HoOeEoiiws7aGi+sLGD9+PAYNGvTfD2dOyPdnuMtNZIbiU1MwbOdmPHqcjkNHD8PLy6tIP8ddbiIzoigKzt6LQ7cNK+HvWh3b+gwqcpgBbqGJzEbMwwcY/dtOXHoQj6+at0Efv0bF/h0MNFEZi3n4AIoC1HF2KVJ/bFIilkUexc/nozCsaQsciLkMa+uSRZOBJioDC49H4LvwUKx4txf2xVzC6jORWPxWN+itrJCWlYm3fXxhq9cjIzsb68/+jvVRpxF1Pw4AUM7aGr0avIzQDz6Bi32FEocZYKCJyoROp8MTRUHQrq1wsa8AAPh4+6bczzecPY0m7rUw5+hB1c+18vTGmvf6QFfILZHFwT+KEZWBoY1z7iZ7mJ6GSwk5L+rrVOdFzO70DgDg8K3rqjC3q10XJ4YEYW3XvmUWZoBb6GLrNnCzas3z0s+3h2lpcLCzg41ej7lvdMHwXVsAAIc+/AzNl89DSEgIsHtrbv/p06fh7+9vtHm4hSYqoYdpafBbNANecyYj22BAG886uZ+1XhEMKysrdO/eHT4+Prh48SIURTFqmAEGmqjETty5lfvfMyL2I+bRA7Tx9Ia9jQ2eKDmPDTpx4gQuXLiAevXEm1OMgbvcRCWgKAo6eNdDU/dacHWohAXHI7D2j1N49Dg9tycrK6tUf7EuCQaaqJhikxLRYdUS9GoQgLjkJLTzrouz9+IQ8ygBADClbSDefamhycMMMNBERXIvJRlf7NkGvU6H0GtXAABLIo/CTm+NyQdz7use3KgpRrZ4HeVtbDSbk4EmKoIph35D2PUYoZ7xJBt1nFzwU7d+cHOopMFkagw0URHEPHyQb73hC27Y3mcQas2aaOKJ8se/chMVwfa+g9Ha0xsA0LymJwDgq+atsavfR9BbmU+MuIUmyiM5IwP7r11BBVtbJGdkwEavh7+rG9wcKqGDtw/2xlyEh2MVDP/7a1qPKig00AbFYKo5iMxGZNwt/GNHSO7arWIlxKUkwcOxCirY2gIA9g/8h1bjFarQQFvpzGdXguRjrk/5bO1ZByE938eEsH2oZGuHhW91xawjB3E/NRk7Lp1HSI/34TV7ktZj5ouJJcpH0xoeWNu1L+xtbNF06VzUruKEHZfO49MmzdG0pofW4xWIgSYqQOVy5bHwra4wAPi/0JxXuI5s8bq2Qz0DA01UiM3no2Cnz3mX1L4BQwt9las5YKDpuacoCs7H38u3Hnw8AokZj9Gxjg9eqvqCBtMVT6GBVp6+aJpIZg/T09B+1RIM3LweAGBQFETfv4uwGzFISEsFAHzdvI2WIxZZoX/ljrt/11RzEGnG2b4CAlyrY9/VSwg+HoGph3KuzXat6ICkzAzUcXKBj0s1jacsmkIDnZD4EJHRv+MV35dNNQ+RyaVkZuD03TsAkBtmALibkgwAmBg8H+69zPMUW16F7nJb6/UYEzwZKU93O4hkVMHGFj19A1Dx6UUjebVt29bEE5WcTlEKPk72q+erNKjrixNnI9G1/Tvw92mIRi/5IyMzA4qioHKlyqab1EzxmWIlZ44XljwxGBCbnIhtF6IxPXw/DFBQWEa0otPpIhVFaZy3/oxruXWYMnwsDkZGIOxEOKb/OAcXr11GekY6bKyt8d0XE/Feu87GmpnI5PRWVnCtWKnA17Wau0K30I0bN1ZOnjyZu76x7xzuJ8Rj6g+zsPnXbQCA63ujjT6kpeFWu2jMcQutKAqWnzqG8Qf2AgCsdDo8MZjfPQ0FbaGLfR66mnNVzP56KnoHdgcAXL5xpQzGIzIPWQYDxh/Yi5kdO8OrshMMZri7XZgSX1gyZfhYAMD5q5fKbBgirdk8vbf5yz3b4O3krPE0xVes+6E92tdXrRcsWIDQLXvRuU1gmQ5FpJXs/9m9dq2Y80ih7OxsTR74VxKluvSzS5cuCDsRjv3HD+XWtoftRrvBb3NXnCySjV6Pgx9+ioq2tmjvnfMs7bNnz2o8VdGV6p8dd3d3LB0/H0PHD0OlCpWQmp6Kewn3AQDtP3oHZ34+AseK2j84jag4aldxxivVa+Le0wtLFi9ejMWLF2s8VdGU+uaMV+oHIGL1r/jmoy+waMwcbFuwAa3/1hK2NjYY8M0QJKYklcWcRCb1ulcdRN2Pg0flKliyZInW4xRZqQ8M/jqurvdmQG7Nr14DKIqCcQunoGtQX2yZtx4Vn75ik8gSvFK9Bsbu34MN3fqj56bViI6Ohq+vr9ZjPZPRbp/U6XQY/+m/0Nj3ZUz/cbaxvobIKGYeDgMABLi5AwA+/vhjLccpMqPfD929w7s4fPoYYu/fMfZXEZWZri/5QQfgXmoyWrZsifDwcK1HKhKjB/rW3VhcuXkVzfu1x7XYG8b+OqIy0flFX7zm6Y3XfgzGO+UrAwDOfvo1YkeMReyIsdoOVwijB3ri4u9y/7vNwEBkZGYa+yuJSs1Kp8P8N3Iu4bV++giirRfM//SVUQN95WYMEhIfAgA8qtcEANgVcIsakblxsrdHZx9fnLpzGwCwLPKoxhM9m1Evf5m3Nufc3YpJi9GqcXPEP8r//UBE5uo1j9rYc+UiGlZzRZQFPMHHaIH+cfNqnIu5gDMhh+Ho4AgAeMHZMh7jQvSXFh5e+HrfdhgUBS72FZCRnQ07M74M1Ci73ClpqZi9KhgrJi7KDTORJapRqTI+fLkJAOBBWiqu/flQ44kKZ5R/alLSUmBrbYMaru7G+PVEJtXHrxHcKzliWeRRmPdTuY20hX7BuRoMigFXboovyCayNPWcq6JngwD8+Tgd3k4uWo9TKKNsoROTE/Eo6U/odFZQFAW6fN42cOTMcRw8GYGRgz43xghkZszx6STFcSc5CVUrVIS1Gb0LOj9Gma58OXsAQLvBb6PHlwOQlZ2l+jwrOwu9vxqIdTs3GuPricrcE4MB1hbwNlajTPi/55pPnD2FuoEBWL9rEwBgzfYNaNClKQAg9Icdxvh6ojK3/dI5NK/lpfUYz2S0f3IWj5kDa7017GztAACjZo+FZwdffDtvAjIyM9C2aStkZGYY6+uJykxmdjZ2XDqHfv6vaD3KMxXrqZ8ltXfpVnQc2iXfz/q+1ROTh40p9XeYu+f9SaCWfAy97UI01p/9Heu69cutuc8Yr+FEZfjUz5Lw8aqLC9tP4VX/Jqr6qMFfYu32DXjy5IkpxiAqNoOiYMHxcLwf8DetRykSkx3ll7O1w5ppy9GpRfvc2rTlMwEA3m/44ZvZY83yDQX0fNt8Pgq2ej06PH2+mLkz6Z/t9Ho9Fo+Zg+t7o/H7pggMeu/93M9+2rUJQ8cPx/2EeFOORFSghLRUTArbh/FtOuV76tUcmeSi1LyP/wUADwDLuzbDcvwbABC+eh+Cf1qC9kM6w8vdE4+SHmHN1OWo6VbDFCMSCb79bRfeq++Hzut+0HqUIjObE2s1Xd0x7fMJOLzmV5y+8Adu3LmFUXPGYut+ntoi09tz5SKi4+9hRLPWWo9SLGYT6L9UKF8BDerkbNEjfj+K0XPG4W+9WmF5yEqNJ6PnyU9Rp/D5q6+hvI2N1qMUi9ncB/a/u+VRl6Mxd+5cBAUFITU9DanpaZi0ZDpu3b2Ncf8YbTHHM2SZnCd9g2NLZ2PDiaNwcnLSepxiMbst9F+GDx+OuaO+U9VWbl2HpNRkjSai50VYWBj8/PwsLsyAGQcaAN55/S3UcqsJz+q1MKTbQABAXLz5PzWCLNvOnTvx5ptvaj1GiZh1oAFg6bh5+DM5ETsO7gYAVHOqqvFEJDNFUbBjxw4EBlrmCxjNOtAe7euj45Au2LBpI2LvxwEAdh3ap/FUJLNrjx4iPT0dfn5+Wo9SImYd6L906NABy5YtAwC8WNsyrtghyxR67TICAwMt9g+vFhFoAOjfvz8cHBwwd81C7DsSqvU4JKnfrl622ONn4BmBTkpKQrVq1ZCSkmKqeQpkZ2eH9evX4+DJCHw09p/w6tgAN+7c1HoskkhqZiZOxcWibdu2Wo9SYoUGOjExEfHx8WjSpAkeP35sqpkKFBgYiP0rdgLI+eNFqw/ewPJN/4ZnB194dvDFnafH2UQlceTWdQS4VYeDg4PWo5RYoYF2c3MDAJw/fx7Hjh0zyUDP4uXuget7o9GofgAAYNLS7wEAtWt4Yel/Vmg4GVm6Swnx8K3qqvUYpVJooK2trdGzZ0/UqlULDx48MKvbG3+esxYbZ67KXV+9fQ3/3rpWw4nI0t1K+hM1HStrPUapFPrEEp1OV/rHlRCRMTxQFKVT3mKhgSYiy2Ixp62I6NkYaCKJMNBEEmGgiSTCQBNJ5P8Bed3ZmLfeKSQAAAAASUVORK5CYII=", "text/plain": [ "
" ] @@ -477,7 +486,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Now set pixels over water to missing_value = -9999\n", + "# Now set pixels over water (without offshore emissions) to missing_value = -9999\n", "if config['LandThreshold']:\n", " state_vector.values[state_vector.isnull()] = -9999" ] @@ -957,7 +966,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -971,7 +980,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.6.10" } }, "nbformat": 4, diff --git a/src/utilities/make_state_vector_file.py b/src/utilities/make_state_vector_file.py index d507037d..b260a3ea 100644 --- a/src/utilities/make_state_vector_file.py +++ b/src/utilities/make_state_vector_file.py @@ -1,6 +1,7 @@ import numpy as np import xarray as xr from sklearn.cluster import KMeans +import yaml def get_nested_grid_bounds(land_cover_pth): @@ -45,42 +46,47 @@ def check_nested_grid_compatibility(lat_min, lat_max, lon_min, lon_max, land_cov def make_state_vector_file( - land_cover_pth, - save_pth, - lat_min, - lat_max, - lon_min, - lon_max, - buffer_deg=5, - land_threshold=0.25, - k_buffer_clust=8, + config_path, land_cover_pth, hemco_diag_pth, save_pth, ): """ Generates the state vector file for an analytical inversion. Arguments + config_path [str] : Path to configuration file land_cover_pth [str] : Path to land cover file + hemco_diag_pth [str] : Path to initial HEMCO diagnostics file save_pth [str] : Where to save the state vector file - lat_min [float] : Minimum latitude - lat_max [float] : Maximum latitude - lon_min [float] : Minimum longitude - lon_max [float] : Maximum longitude - buffer_deg [float] : Width of k-means buffer area in degrees - land_threshold [float] : Minimum land fraction to include pixel as a state vector element - k_buffer_clust [int] : Number of buffer clusters for k-means Returns ds_statevector [] : xarray dataset containing state vector field formatted for HEMCO Notes - - Land cover file looks like 'GEOSFP.20200101.CN.025x03125.NA.nc' + - Land cover file looks like 'GEOSFP.20200101.CN.025x03125.NA.nc' (or 0.5-deg equivalent) + - HEMCO diags file needs to be global, is used to include offshore emissions in state vector + - Land cover file and HEMCO diags file need to have the same grid resolution """ - # Load land cover data + # Get config + config = yaml.load(open(config_path), Loader=yaml.FullLoader) + lat_min = config["LatMin"] + lat_max = config["LatMax"] + lon_min = config["LonMin"] + lon_max = config["LonMax"] + buffer_deg = config["BufferDeg"] + land_threshold = config["LandThreshold"] + emis_threshold = config["OffshoreEmisThreshold"] + k_buffer_clust = config["nBufferClusters"] + + # Load land cover data and HEMCO diagnostics lc = xr.load_dataset(land_cover_pth) + hd = xr.load_dataset(hemco_diag_pth) - # Group fields together + # Require hemco diags on same global grid as land cover map + hd["lon"] = hd["lon"] - 0.03125 # initially offset by 0.03125 degrees + + # Select / group fields together 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 with nesting window compatible = check_nested_grid_compatibility( @@ -88,7 +94,7 @@ def make_state_vector_file( ) if not compatible: raise ValueError( - "Region of interest not contained within selected NestedRegion; see config.yml)." + "Region of interest not contained within selected NestedRegion; see config.yml." ) # Define bounds of inversion domain @@ -103,35 +109,30 @@ def make_state_vector_file( lat_min_inv_domain = np.max([lat_min - buffer_deg, minLat_allowed]) lat_max_inv_domain = np.min([lat_max + buffer_deg, maxLat_allowed]) - # Subset inversion domain using land cover file + # Subset inversion domain for land cover and hemco diagnostics fields lc = lc.isel(lon=lc.lon >= lon_min_inv_domain, lat=lc.lat >= lat_min_inv_domain) lc = lc.isel(lon=lc.lon <= lon_max_inv_domain, lat=lc.lat <= lat_max_inv_domain) + hd = hd.isel(lon=hd.lon >= lon_min_inv_domain, lat=hd.lat >= lat_min_inv_domain) + hd = hd.isel(lon=hd.lon <= lon_max_inv_domain, lat=hd.lat <= lat_max_inv_domain) - # Replace all values with NaN (to be filled later) + # Initialize state vector from land cover, replacing all values with NaN (to be filled later) statevector = lc.where(lc == -9999.0) # Set pixels in buffer areas to 0 statevector[:, (statevector.lon < lon_min) | (statevector.lon > lon_max)] = 0 statevector[(statevector.lat < lat_min) | (statevector.lat > lat_max), :] = 0 - # Also set pixels over water to 0 + # Also set pixels over water to 0, unless there are offshore emissions if land_threshold: - # Where there is no land, replace with 0 - land = lc.where(lc > land_threshold) - statevector.values[land.isnull().values] = 0 + # Where there is neither land nor emissions, replace with 0 + land = lc.where((lc > land_threshold) | (hd > emis_threshold)) + statevector.values[land.isnull().values] = -9999 # Fill in the remaining NaNs with state vector element values statevector.values[statevector.isnull().values] = np.arange( 1, statevector.isnull().sum() + 1 )[::-1] - # Now set pixels over water to missing_value = -9999 - if land_threshold: - # First, where there is no land, replace with NaN - statevector = statevector.where(lc > land_threshold) - # Fill with missing_value = -9999 - statevector.values[statevector.isnull().values] = -9999 - # Assign buffer pixels (the remaining 0's) to state vector # ------------------------------------------------------------------------- buffer_area = statevector.values == 0 @@ -175,7 +176,7 @@ def make_state_vector_file( ds_statevector.StateVector.attrs["_FillValue"] = -9999 # Save - if save_pth: + if save_pth is not None: print("Saving file {}".format(save_pth)) ds_statevector.to_netcdf(save_pth) @@ -185,24 +186,13 @@ def make_state_vector_file( if __name__ == "__main__": import sys - land_cover_pth = sys.argv[1] - save_pth = sys.argv[2] - lat_min = float(sys.argv[3]) - lat_max = float(sys.argv[4]) - lon_min = float(sys.argv[5]) - lon_max = float(sys.argv[6]) - buffer_deg = float(sys.argv[7]) - land_threshold = float(sys.argv[8]) - k_buffer_clust = int(sys.argv[9]) - + config_path = sys.argv[1] + land_cover_pth = sys.argv[2] + hemco_diag_pth = sys.argv[3] + save_pth = sys.argv[4] make_state_vector_file( - land_cover_pth, + config_path, + land_cover_pth, + hemco_diag_pth, save_pth, - lat_min, - lat_max, - lon_min, - lon_max, - buffer_deg, - land_threshold, - k_buffer_clust, ) diff --git a/src/utilities/sanitize_input_yaml.py b/src/utilities/sanitize_input_yaml.py index cb486238..4026d3bc 100644 --- a/src/utilities/sanitize_input_yaml.py +++ b/src/utilities/sanitize_input_yaml.py @@ -36,6 +36,7 @@ "nBufferClusters", "BufferDeg", "LandThreshold", + "OffshoreEmisThreshold", "ReducedDimensionStateVector", "StateVectorFile", "ShapeFile", diff --git a/src/utilities/utils.py b/src/utilities/utils.py index 7353451d..392caf20 100644 --- a/src/utilities/utils.py +++ b/src/utilities/utils.py @@ -9,7 +9,9 @@ def download_landcover_files(config): # conditionally set variables to create s3 and landcover file paths gridDir = ( - f"{config['Res']}_{config['NestedRegion']}" if len(config["NestedRegion"]) == 2 else "" + f"{config['Res']}_{config['NestedRegion']}" + if len(config["NestedRegion"]) == 2 + else "" ) if config["Met"] == "geosfp": @@ -45,3 +47,33 @@ def download_landcover_files(config): else results ) print(output) + + +def download_hemcodiags_files(config): + """ + Download global hemco diagnostics files from s3 given the config file + """ + DataPath = "/home/ubuntu/ExtData" + + if config["Res"] == "4x5": + gridRes = "${Res}" + elif config["Res"] == "2x2.5": + gridRes = "2x25" + elif config["Res"] == "0.5x0.625": + gridRes = "05x0625" + elif config["Res"] == "0.25x0.3125": + gridRes = "025x03125" + + HemcoDiagFile = f"{DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.{gridRes}.20190101.nc" + s3_hd_path = f"s3://gcgrid/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.{gridRes}.20190101.nc" + + # run the aws command to download the files + command = f"aws s3 cp --request-payer=requester {s3_hd_path} {HemcoDiagFile}" + results = subprocess.run(command.split(), capture_output=True, text=True) + + output = ( + "Successfully downloaded hemco diags files" + if results.returncode == 0 + else results + ) + print(output)