Plot_irrigated_fraction_and_irrigation_water_withdrawal_maps.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "351c8bae",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scipy.io as scio\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import cartopy.feature as cfeature\n",
    "import cartopy.crs as ccrs\n",
    "import netCDF4 as nc\n",
    "import numpy as np\n",
    "from Load_data import Data_from_nc\n",
    "import xarray as xr\n",
    "import math\n",
    "import gc\n",
    "\n",
    "def get_data_from_mat_for_calcu(file, variable):\n",
    "    var_dict = scio.loadmat(file)\n",
    "    var = var_dict[variable]\n",
    "    var = var[:, 29:]\n",
    "    return var.T\n",
    "\n",
    "def get_data_from_nc(file, variable):\n",
    "    file_obj = nc.Dataset(file)\n",
    "    data = file_obj.variables[variable]\n",
    "    var_data = np.array(data)\n",
    "    var_data = var_data[:, 29:, :]\n",
    "    var_data[var_data > 1000000] = np.nan\n",
    "    var_data = np.squeeze(var_data)\n",
    "    return var_data\n",
    "\n",
    "def plot_fre_single(ax, data_xarray, title, metric):\n",
    "    bwr = mpl.cm.get_cmap('Blues')\n",
    "    vmin=0\n",
    "    vmax=0.05\n",
    "    alpha=0.05\n",
    "    colors = [bwr(0),bwr(0.2),bwr(0.3), bwr(0.4), bwr(0.5), bwr(0.6), bwr(0.7), bwr(0.8),bwr(0.9), bwr(1.0)]\n",
    "    levels = [1, 10, 25, 50, 75, 100, 150, 200, 250]\n",
    "    cmap=mpl.colors.ListedColormap(colors)\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    \n",
    "    cb = plt.colorbar(im,fraction= 0.02, pad= 0.04, extend='both')\n",
    "    cb.ax.tick_params(labelsize=16)\n",
    "    ax.coastlines(color='white', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=18)\n",
    "    ax.set_title(metric, loc='left', fontsize=18)\n",
    "    \n",
    "\n",
    "def plot_fre_single_no_cb(ax, data_xarray, title, metric):\n",
    "    bwr = mpl.cm.get_cmap('Blues')\n",
    "    vmin=0\n",
    "    vmax=0.05\n",
    "    alpha=0.05\n",
    "    colors = [bwr(0),bwr(0.2),bwr(0.3), bwr(0.4), bwr(0.5), bwr(0.6), bwr(0.7), bwr(0.8),bwr(0.9), bwr(1.0)]\n",
    "    levels = [1, 10, 25, 50, 75, 100, 150, 200, 250]\n",
    "    cmap=mpl.colors.ListedColormap(colors)\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    \n",
    "    #cb = plt.colorbar(im,fraction= 0.02, pad= 0.04, extend='both')\n",
    "    ax.coastlines(color='white', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=18)\n",
    "    ax.set_title(metric, loc='left', fontsize=18)\n",
    "    \n",
    "def plot_fre_single_1(ax, data_xarray, title, metric):\n",
    "    bwr = mpl.cm.get_cmap('Reds')\n",
    "    vmin=0\n",
    "    vmax=0.05\n",
    "    alpha=0.05\n",
    "    colors = ['white','b','c', 'g', 'm', 'r', 'y', 'orange', 'pink', 'purple']\n",
    "    levels = [1, 4, 5, 7, 20, 33, 36, 38, 39]\n",
    "    cmap=mpl.colors.ListedColormap(colors)\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    \n",
    "    cb = plt.colorbar(im,fraction= 0.02, pad= 0.04, extend='both')\n",
    "    cb.ax.tick_params(labelsize=16)\n",
    "    ax.coastlines(color='white', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=18)\n",
    "    ax.set_title(metric, loc='left', fontsize=18)\n",
    "    \n",
    "def plot_fre_single_2(ax, data_xarray, title, metric):\n",
    "    bwr = mpl.cm.get_cmap('Reds')\n",
    "    vmin=0\n",
    "    vmax=0.05\n",
    "    alpha=0.05\n",
    "    colors = ['white',bwr(0.2),bwr(0.3), bwr(0.4), bwr(0.5), bwr(0.6), bwr(0.7), bwr(0.8),bwr(0.9), bwr(1.0)]\n",
    "    levels = [0.01, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40]\n",
    "    cmap=mpl.colors.ListedColormap(colors)\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    \n",
    "    cb = plt.colorbar(im,fraction= 0.02, pad= 0.04, extend='both')\n",
    "    cb.ax.tick_params(labelsize=16)\n",
    "    cb.set_ticklabels(['1%','5%','10%','15%','20%','25%','30%','35%','40%'])\n",
    "    ax.coastlines(color='white', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=18)\n",
    "    ax.set_title(metric, loc='left', fontsize=18)\n",
    "    \n",
    "def plot_land_single(ax, data_xarray, title, metric):\n",
    "    bwr = mpl.cm.get_cmap('Greens')\n",
    "    vmin=0\n",
    "    vmax=0.05\n",
    "    alpha=0.05\n",
    "    colors = ['white',bwr(0.2),bwr(0.3), bwr(0.4), bwr(0.5), bwr(0.6), bwr(0.7), bwr(0.8),bwr(0.9), bwr(1.0)]\n",
    "    levels = [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]\n",
    "    cmap=mpl.colors.ListedColormap(colors)\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    \n",
    "    cb = plt.colorbar(im,fraction= 0.02, pad= 0.04, extend='both')\n",
    "    \n",
    "    cb.ax.tick_params(labelsize=16)\n",
    "    cb.set_ticklabels(['1%','10%','20%','30%','40%','50%','60%','70%','80%'])\n",
    "    ax.coastlines(color='dimgray', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=18)\n",
    "    ax.set_title(metric, loc='left', fontsize=18)\n",
    "\n",
    "def get_data_from_mat(file, variable):\n",
    "    var_dict = scio.loadmat(file)\n",
    "    var = var_dict[variable]\n",
    "    var = var[:, 29:] \n",
    "    return var\n",
    "\n",
    "\n",
    "def get_data_from_mat_1(file, variable):\n",
    "    var_dict = scio.loadmat(file)\n",
    "    var = var_dict[variable]\n",
    "    return var\n",
    "\n",
    "ar6_region = get_data_from_mat_for_calcu('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/ar6_region.mat','ar6_region')\n",
    "\n",
    "irr_diff = get_data_from_mat_for_calcu('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/irr_diff_out.mat', 'irr_diff_out')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "907b7fba",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_surface = Data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/surfdata_irrigation_method.nc')   #load the data\n",
    "data_irrigation_method = data_surface.load_variable('irrigation_method')\n",
    "data_lon = data_surface.load_variable('LONGXY')\n",
    "data_lon = data_lon[0,:]\n",
    "data_lat = data_surface.load_variable('LATIXY')\n",
    "data_lat = data_lat[29:,0]\n",
    "data_lat = data_lat[:]\n",
    "\n",
    "IWW_IRR_1901_1930_CESM2 = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/QIRRIG/CESM2_IRR_YEARLYQIRRIG_timmean_1901_1930.nc', 'QIRRIG')\n",
    "IWW_IRR_1985_2014_CESM2 = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/QIRRIG/CESM2_IRR_YEARLYQIRRIG_timmean_1985_2014.nc', 'QIRRIG')\n",
    "\n",
    "IWW_IRR_1901_1930_CESM2_gw = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2_gw/QIRRIG/CESM2_gw_IRR_YEARLYQIRRIG_timmean_1901_1930.nc', 'QIRRIG')\n",
    "IWW_IRR_1985_2014_CESM2_gw = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2_gw/QIRRIG/CESM2_gw_IRR_YEARLYQIRRIG_timmean_1985_2014.nc', 'QIRRIG')\n",
    "\n",
    "IWW_IRR_1901_1930_NorESM = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/NorESM/QIRRIG/NorESM_IRR_YEARLYQIRRIG_timmean_1901_1930.nc', 'QIRRIG')\n",
    "IWW_IRR_1985_2014_NorESM = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/NorESM/QIRRIG/NorESM_IRR_YEARLYQIRRIG_timmean_1985_2014.nc', 'QIRRIG')\n",
    "\n",
    "IWW_IRR_1901_1930_IPSL = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/IPSL-CM6/IRR01/after_processing/IRR01_1901_2014_irr_Day_Month.nc_1901_1930_timmean_YEARLY_0.9x1.25', 'irr')\n",
    "IWW_IRR_1985_2014_IPSL = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/IPSL-CM6/IRR01/after_processing/IRR01_1901_2014_irr_Day_Month.nc_1985_2014_timmean_YEARLY_0.9x1.25', 'irr')\n",
    "\n",
    "IWW_IRR_1901_1930_E3SM = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/Water_Energy/IRR_QIRRIG_REAL_timmean_1901_1930.nc_0.9x1.25_YEARLY', 'QIRRIG_REAL')\n",
    "IWW_IRR_1985_2014_E3SM = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/Water_Energy/IRR_QIRRIG_REAL_timmean_1985_2014.nc_0.9x1.25_YEARLY', 'QIRRIG_REAL')\n",
    "\n",
    "IWW_IRR_1901_1930_CNRM = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/CNRM/CNRM-CM6-1_hist-irr_r1i1p1f2/after_processing/airrww_mon_CNRM-CM6-1_hist-irr_r1i1p1f2_gr_190101-201412.nc_yearmean_YEARLY_1901_1930_timmean_0.9x1.25', 'airrww')\n",
    "IWW_IRR_1985_2014_CNRM = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/CNRM/CNRM-CM6-1_hist-irr_r1i1p1f2/after_processing/airrww_mon_CNRM-CM6-1_hist-irr_r1i1p1f2_gr_190101-201412.nc_yearmean_YEARLY_1985_2014_timmean_0.9x1.25', 'airrww')\n",
    "\n",
    "# IWW_IRR_1901_1930_MIROC = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/MIROC/Water_energy/irrac_IRR01_1901_2014_yearmean.nc_1901_1930_timmean_0.9x1.25_YEARLY', 'irrac')\n",
    "# IWW_IRR_1985_2014_MIROC = get_data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/MIROC/Water_energy/irrac_IRR01_1901_2014_yearmean.nc_1985_2014_timmean_0.9x1.25_YEARLY', 'irrac')\n",
    "\n",
    "IWW_IRR_1901_1930 = (IWW_IRR_1901_1930_CESM2+IWW_IRR_1901_1930_CESM2_gw+IWW_IRR_1901_1930_NorESM+IWW_IRR_1901_1930_IPSL+IWW_IRR_1901_1930_E3SM+IWW_IRR_1901_1930_CNRM)/6\n",
    "IWW_IRR_1985_2014 = (IWW_IRR_1985_2014_CESM2+IWW_IRR_1985_2014_CESM2_gw+IWW_IRR_1985_2014_NorESM+IWW_IRR_1985_2014_IPSL+IWW_IRR_1985_2014_E3SM+IWW_IRR_1985_2014_CNRM)/6\n",
    "\n",
    "\n",
    "globe_irr_land_1901 = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/input_data/irrigated_GRID_1901.mat', 'irrigated_GRID_1901')\n",
    "globe_irr_land_1941 = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/input_data/irrigated_GRID_1941.mat', 'irrigated_GRID_1941')\n",
    "globe_irr_land_1981 = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/input_data/irrigated_GRID_1981.mat', 'irrigated_GRID_1981')\n",
    "globe_irr_land_2014 = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/input_data/irrigated_GRID_2014.mat', 'irrigated_GRID_2014')\n",
    "\n",
    "globe_irr_land_1901[np.isnan(IWW_IRR_1901_1930_CESM2.T)] = np.nan\n",
    "globe_irr_land_1941[np.isnan(IWW_IRR_1901_1930_CESM2.T)] = np.nan\n",
    "globe_irr_land_1981[np.isnan(IWW_IRR_1901_1930_CESM2.T)] = np.nan\n",
    "globe_irr_land_2014[np.isnan(IWW_IRR_1901_1930_CESM2.T)] = np.nan\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce30921c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_xr_1 = xr.DataArray(globe_irr_land_1901.T, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "fig = plt.figure(figsize=(20,8),dpi=300)\n",
    "fig.subplots_adjust(hspace=0.1, wspace=0.05, left = 0.05, right = 0.95, top = 0.95, bottom = 0.05)\n",
    "\n",
    "proj=ccrs.PlateCarree()\n",
    "\n",
    "\n",
    "irr_diff = (irr_diff + IWW_IRR_1901_1930 - IWW_IRR_1901_1930)/100\n",
    "data_xr_1 = xr.DataArray(irr_diff, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "ax = plt.subplot(221, projection=proj, frameon=False)\n",
    "ax.text(0.005, 0.95, 'a', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "plot_fre_single_2(ax, data_xr_1, 'irrigation expansion 1901-2014', 'Percentage of the grid area')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "\n",
    "ar6_region = ar6_region + IWW_IRR_1901_1930_CESM2 - IWW_IRR_1901_1930_CESM2 - 0.01\n",
    "data_xr_1 = xr.DataArray(ar6_region, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "ax = plt.subplot(222, projection=proj, frameon=False)\n",
    "ax.text(0.005, 0.95, 'b', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "plot_fre_single_1(ax, data_xr_1, '', 'AR6_region')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "\n",
    "data_xr_1 = xr.DataArray(IWW_IRR_1901_1930, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "ax = plt.subplot(223, projection=proj, frameon=False)\n",
    "ax.text(0.005, 0.95, 'c', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "plot_fre_single(ax, data_xr_1, 'tranirr 1901-1930', 'Irrigation water withdrawal ($\\mathregular{mm/yr}$)')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "\n",
    "\n",
    "data_xr_1 = xr.DataArray(IWW_IRR_1985_2014, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "ax = plt.subplot(224, projection=proj, frameon=False)\n",
    "ax.text(0.005, 0.95, 'd', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "plot_fre_single(ax, data_xr_1, 'tranirr 1985-2014', 'Irrigation water withdrawal ($\\mathregular{mm/yr}$)')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ddf4bab",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_xr_1 = xr.DataArray(globe_irr_land_1901.T, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "fig = plt.figure(figsize=(20,8),dpi=300)\n",
    "fig.subplots_adjust(hspace=0.1, wspace=0.05, left = 0.05, right = 0.95, top = 0.95, bottom = 0.05)\n",
    "\n",
    "proj=ccrs.PlateCarree()\n",
    "ax = plt.subplot(221, projection=proj, frameon=False)\n",
    "plot_land_single(ax, data_xr_1, '1901', 'irrigated fraction')\n",
    "ax.text(0.005, 0.95, 'a', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "\n",
    "data_xr_1 = xr.DataArray(globe_irr_land_1941.T, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "ax = plt.subplot(222, projection=proj, frameon=False)\n",
    "ax.text(0.005, 0.95, 'b', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "plot_land_single(ax, data_xr_1, '1941', 'irrigated fraction')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "\n",
    "data_xr_1 = xr.DataArray(globe_irr_land_1981.T, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "ax = plt.subplot(223, projection=proj, frameon=False)\n",
    "ax.text(0.005, 0.95, 'c', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "plot_land_single(ax, data_xr_1, '1981', 'irrigated fraction')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "\n",
    "data_xr_1 = xr.DataArray(globe_irr_land_2014.T, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "ax = plt.subplot(224, projection=proj, frameon=False)\n",
    "ax.text(0.005, 0.95, 'd', color='dimgrey', fontsize=16, transform=ax.transAxes, weight='bold')\n",
    "plot_land_single(ax, data_xr_1, '2014', 'irrigated fraction')\n",
    "ax.add_feature(cfeature.OCEAN, color='lightgrey')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}