0uncertainties estimation_regional error propagation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e3e53333-d220-4ec3-bcd6-6f32c2d0519f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geoutils as gu\n",
    "import numpy as np\n",
    "import xdem\n",
    "import os\n",
    "import shutil\n",
    "import xarray as xr\n",
    "import rioxarray as rioxr\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import proplot as pplt\n",
    "from xdem.spatialstats import nmad\n",
    "import pandas as pd\n",
    "from matplotlib.ticker import FuncFormatter\n",
    "import gc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e330f9c4-8f3c-42fb-a892-b82a450ba985",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "09f610fa-27e8-44a1-9d82-7f4063204da5",
   "metadata": {},
   "outputs": [],
   "source": [
    "pth_vgm='D:\\\\Key_works\\\\20231212dh_calculate\\\\xnhl_test\\\\spatial_corr_analysis\\\\vgm_params'\n",
    "pth_err='D:\\\\Key_works\\\\20231212dh_calculate\\\\xnhl_test\\\\group_err'\n",
    "pth_gla='D:\\\\Key_works\\\\20231212dh_calculate\\\\gla\\\\CGI1'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f7b1f69-25b7-4d61-9862-4d30a294b535",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### test for one group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "26d7ec98-e83a-4c90-a4c5-46eceb60c008",
   "metadata": {},
   "outputs": [],
   "source": [
    "pth1=pth_err+'\\\\g0001_heteroscedasticity.tif'\n",
    "pth2=pth_vgm+'\\\\g0001_corr.csv'\n",
    "pth3=pth_gla+'\\\\g0001_gla.shp'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "e5bcb586-88ba-4e60-9ac1-79945cea6d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "params_variogram_model=pd.read_csv(pth2,index_col=0)\n",
    "errors=xdem.DEM(pth1)\n",
    "glacier_outlines = gpd.read_file(pth3)\n",
    "\n",
    "merged_gdf = gpd.GeoDataFrame(geometry=[glacier_outlines.geometry.unary_union], crs=glacier_outlines.crs)\n",
    "#mask_glacier = glacier_outlines.create_mask(dh)\n",
    "\n",
    "stderr_glaciers = xdem.spatialstats.spatial_error_propagation(\n",
    "    areas=[merged_gdf], errors=errors, params_variogram_model=params_variogram_model\n",
    ")\n",
    "\n",
    "stderr_glaciers1 = xdem.spatialstats.spatial_error_propagation(\n",
    "    areas=merged_gdf.geometry.area.values, errors=errors, params_variogram_model=params_variogram_model\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "a3f485c4-ed93-48e0-ae26-8a3f8c20aabb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([3.5554939606581177], [3.5115590323268253])"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stderr_glaciers,stderr_glaciers1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df215073-2b6a-4fe6-90d4-6652a3fb4588",
   "metadata": {},
   "source": [
    "### all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2ece2255-45d5-4c8d-a309-a7e420eaf079",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "261"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glst=[i[0:5] for i in os.listdir(pth_vgm)]\n",
    "len(glst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "db8601d8-d6ce-4fbf-8dac-5ad5986dfd93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g0423 done!\n",
      "g0424 done!\n"
     ]
    }
   ],
   "source": [
    "#out_lst=[]\n",
    "for gpi in glst:\n",
    "    out_lst=[]\n",
    "    if os.path.exists('tmp_file\\\\'+gpi+'.csv'):\n",
    "        continue\n",
    "    print(gpi,end='')\n",
    "    parai=pd.read_csv(pth_vgm+'\\\\'+gpi+'_corr.csv',index_col=0)\n",
    "    erri=xdem.DEM(pth_err+'\\\\'+gpi+'_heteroscedasticity.tif')\n",
    "    gla_gpi=gpd.read_file(pth_gla+'\\\\'+gpi+'_gla.shp',index_col=0)\n",
    "    merged_gdf = gpd.GeoDataFrame(geometry=[gla_gpi.geometry.unary_union], crs=gla_gpi.crs)\n",
    "\n",
    "    try:\n",
    "        sign='poly'\n",
    "        stderr_glaciers = xdem.spatialstats.spatial_error_propagation(\n",
    "            areas=[merged_gdf], errors=erri, params_variogram_model=parai)\n",
    "    except:\n",
    "        sign='value'\n",
    "        stderr_glaciers = xdem.spatialstats.spatial_error_propagation(\n",
    "            areas=merged_gdf.geometry.area.values, errors=erri, params_variogram_model=parai)\n",
    "    out_lst.append([gpi,stderr_glaciers[0],sign])\n",
    "    dfi=pd.DataFrame(out_lst,columns=['group','error_1sigma','calculation_form'])\n",
    "    dfi.to_csv('tmp_file\\\\'+gpi+'.csv')\n",
    "    print(' done!')\n",
    "    gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0b740268-a48e-4ad5-908c-0b371d9ae704",
   "metadata": {},
   "outputs": [],
   "source": [
    "flst=[i for i in os.listdir('tmp_file') if i[-4:]=='.csv']\n",
    "dfflst=[]\n",
    "for nm in flst:\n",
    "    dffi=pd.read_csv('tmp_file\\\\'+nm)\n",
    "    dfflst.append(dffi)\n",
    "df_all=pd.concat(dfflst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "142f2e61-e02f-4096-8e43-60fb62e5ea45",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all=df_all.drop_duplicates()[['group','error_1sigma','calculation_form']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0c698954-ad66-4f14-90bb-be5db82aa490",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all.to_csv('D:\\\\Key_works\\\\20231212dh_calculate\\\\xnhl_test\\\\dh_error_propagation.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcf8f315-7a39-4d24-9806-8555bb99c9f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all=pd.DataFrame(out_lst,columns=['group','error_1sigma','calculation_form'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9f246546-800d-4f57-9452-bf0f1c4cac77",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "203"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(df_all[df_all['calculation_form']=='poly'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "723ff7e8-0590-4042-86ff-110c4707a1ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>group</th>\n",
       "      <th>error_1sigma</th>\n",
       "      <th>calculation_form</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>g0001</td>\n",
       "      <td>3.551806</td>\n",
       "      <td>poly</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>g0002</td>\n",
       "      <td>1.934397</td>\n",
       "      <td>poly</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>g0004</td>\n",
       "      <td>4.056870</td>\n",
       "      <td>poly</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>g0005</td>\n",
       "      <td>2.301977</td>\n",
       "      <td>poly</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>g0006</td>\n",
       "      <td>7.864847</td>\n",
       "      <td>poly</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   group  error_1sigma calculation_form\n",
       "0  g0001      3.551806             poly\n",
       "1  g0002      1.934397             poly\n",
       "2  g0004      4.056870             poly\n",
       "3  g0005      2.301977             poly\n",
       "4  g0006      7.864847             poly"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_all.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e09014f8-1a1b-4425-87c3-032054b9c63c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py39",
   "language": "python",
   "name": "py39"
  },
  "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.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}