run_DEM_EUDEMv11_ipy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''Python for Processing Digital Elevation Models (DEMs).\n",
    "\n",
    "Author: He Zhang @ University of Exeter\n",
    "Date: 16th March 2019 (Update: 20th April 2019)\n",
    "Contact: hz298@exeter.ac.uk zhangheupc@126.com\n",
    "\n",
    "Copyright (c) 2019 He Zhang\n",
    "'''\n",
    "\n",
    "'''\n",
    "DEM Data:\n",
    "    EUDEMv11  -  EU-DEMv1.1\n",
    "    Download Link: https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1\n",
    "\n",
    "Terms and Abbreviations:\n",
    "    EPSG  -  European Petroleum Survey Group\n",
    "    GCS   -  Geographic Coordinate System (Identified by an unique EPSG code)\n",
    "    PCS   -  Projected Coordinate System (Identified by an unique EPSG code)\n",
    "    WGS-84 [EPSG 4326]   -  1984 World Geodetic System [GCS]\n",
    "    Merc [EPSG 3857]     -  Mercator -> Web Mercator -> Pseudo Mercator [PCS of WGS-84]\n",
    "    OSGB-36 [EPSG 4277]  -  1936 Ordnance Survey Great Britain [GCS]\n",
    "    BNG [EPSG 27700]     -  British National Grid [PCS of OSGB-36]\n",
    "    ETRS-89 [EPSG 4258]  -  1989 European Terrestrial Reference System [GCS]\n",
    "    LAEA [EPSG 3035]     -  Lambert Azimuthal Equal-Area [PCS of ETRS-89]\n",
    "    lat  -  Latitude\n",
    "    lng  -  Longitude\n",
    "    Transform  -  GCS to GCS\n",
    "    Project    -  GCS to PCS, PCS to PCS\n",
    "    Convert    -  PCS to PCS\n",
    "\n",
    "Functions:\n",
    "    Convert DEM from LAEA PCS to ETRS-89 GCS.\n",
    "\n",
    "    Transform DEM from ETRS-89 GCS to WGS-84 GCS.\n",
    "    Transform DEM from ETRS-89 GCS to OSGB-36 GCS.\n",
    "\n",
    "    Get the elevation from DEM in WGS-84 GCS.\n",
    "    Get the elevation from DEM in OSGB-36 GCS.\n",
    "    Get the elevation from DEM in ETRS-89 GCS.\n",
    "\n",
    "******************** Important Information of Code Usage ********************\n",
    "- Use 'GDAL.GetProjection()' to check the GCS/PCS information of DEM (in TIF format).\n",
    "- Use 'GDAL.GetGeoTransform()' to check the resolution of DEMs.\n",
    "- Use 'GDAL.GetGeoTransform()' to check (lat, lng) of top-left corner of DEM (in GCS).\n",
    "- The (lat, lng) of other locations in DEM can therefore be calculated.\n",
    "- Each GCS has one related PCS (GCS/PCS is identified by an unique EPSG code).\n",
    "- You can transform DEM between different GCSs (e.g., WGS-84 <-> OSGB-36 <-> ETRS-89).\n",
    "- You can project DEM in GCS to the related PCS and then display its 2D image.\n",
    "- You can not display DEM in GCS as 2D image (e.g., WGS-84 -> 2D image is wrong).\n",
    "- You can not project DEM in GCS to the unrelated PCS (e.g., WGS-84 -> BNG is wrong).\n",
    "- You can not project DEM between different PCSs (e.g., Pseudo Mercator <-> BNG is wrong).\n",
    "* Use 'gdalwarp' command to transform/project DEM to different GCSs/PCSs might be correct.\n",
    "'''\n",
    "\n",
    "# Python 3.7\n",
    "\n",
    "# import os\n",
    "# import re\n",
    "# import shutil\n",
    "# import subprocess\n",
    "\n",
    "# import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from osgeo import gdal\n",
    "from pandas import read_csv\n",
    "\n",
    "from pyDEM_function import get_dem_info\n",
    "from pyDEM_function import get_elevation\n",
    "# from pyDEM_function import get_file_names\n",
    "# from pyDEM_function import show_2d_dem\n",
    "from pyDEM_function import transprojcnvt_dem\n",
    "# from pyDEM_function import write_dem\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify user settings.\n",
    "\n",
    "# Set the EPSG code.\n",
    "EPSG_WGS84 = 4326\n",
    "EPSG_OSGB36 = 4277\n",
    "EPSG_ETRS89 = 4258\n",
    "EPSG_LAEA = 3035\n",
    "\n",
    "# Set the format of DEMs.\n",
    "DEM_FORMAT = '.tif'\n",
    "\n",
    "# Set the path of DEMs.\n",
    "PATH_EUDEM = 'DATA/DATA_EUDEMv11/'\n",
    "# PATH_EUDEM_SOURCE = 'DATA/DATA_EUDEMv11/EPSG3035_s/'  # The folder of source DEMs must exist.\n",
    "\n",
    "# Set the name of DEMs in GCSs.\n",
    "EUDEM_GCS_WD = 'EUDEMv11_EPSG4326.tif'\n",
    "EUDEM_GCS_UK = 'EUDEMv11_EPSG4277.tif'\n",
    "EUDEM_GCS_EU = 'EUDEMv11_EPSG4258.tif'\n",
    "\n",
    "# Set the name of DEMs in PCSs.\n",
    "EUDEM_PCS_WD = 'EUDEMv11_EPSG3857.tif'\n",
    "EUDEM_PCS_UK = 'EUDEMv11_EPSG27700.tif'\n",
    "EUDEM_PCS_EU = 'EUDEMv11_EPSG3035.tif'\n",
    "\n",
    "# Set the path of location data file.\n",
    "PATH_LD_STATION_DATA = 'DATA/DATA_LD_AirQuality/London_AirQuality_Stations.csv'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> <EUDEMv11> Convert DEM from LAEA PCS to ETRS-89 GCS.\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 40000\n",
      "The number of column (width) is: 40000\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (3000000.0, 25.0, 0.0, 4000000.0, 0.0, -25.0)\n",
      "The GCS/PCS information is:\n",
      " PROJCS[\"ETRS89_ETRS_LAEA\",GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4258\"]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",52],PARAMETER[\"longitude_of_center\",10],PARAMETER[\"false_easting\",4321000],PARAMETER[\"false_northing\",3210000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30505\n",
      "The number of column (width) is: 52428\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4258\"]]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# <EUDEMv11> Convert DEM from LAEA PCS to ETRS-89 GCS.\n",
    "\n",
    "print('\\n>>> <EUDEMv11> Convert DEM from LAEA PCS to ETRS-89 GCS.')\n",
    "\n",
    "path = PATH_EUDEM\n",
    "dem_in = path + EUDEM_PCS_EU\n",
    "dem_out = path + EUDEM_GCS_EU\n",
    "\n",
    "epsg_in = EPSG_LAEA\n",
    "epsg_out = EPSG_ETRS89\n",
    "\n",
    "transprojcnvt_dem(dem_in, epsg_in, dem_out, epsg_out)\n",
    "\n",
    "data = gdal.Open(dem_in)\n",
    "get_dem_info(data, if_print=True)\n",
    "data = gdal.Open(dem_out)\n",
    "get_dem_info(data, if_print=True)\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> <EUDEMv11> Transform DEM from ETRS-89 GCS to WGS-84 GCS.\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30505\n",
      "The number of column (width) is: 52428\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4258\"]]\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30505\n",
      "The number of column (width) is: 52428\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# <EUDEMv11> Transform DEM from ETRS-89 GCS to WGS-84 GCS.\n",
    "\n",
    "print('\\n>>> <EUDEMv11> Transform DEM from ETRS-89 GCS to WGS-84 GCS.')\n",
    "\n",
    "path = PATH_EUDEM\n",
    "dem_in = path + EUDEM_GCS_EU\n",
    "dem_out = path + EUDEM_GCS_WD\n",
    "\n",
    "epsg_in = EPSG_ETRS89\n",
    "epsg_out = EPSG_WGS84\n",
    "\n",
    "transprojcnvt_dem(dem_in, epsg_in, dem_out, epsg_out)\n",
    "\n",
    "data = gdal.Open(dem_in)\n",
    "get_dem_info(data, if_print=True)\n",
    "data = gdal.Open(dem_out)\n",
    "get_dem_info(data, if_print=True)\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> <EUDEMv11> Transform DEM from ETRS-89 GCS to OSGB-36 GCS.\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30505\n",
      "The number of column (width) is: 52428\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4258\"]]\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30506\n",
      "The number of column (width) is: 52430\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266102928602475, 0.0003392475335051166, 0.0, 58.98812292836723, 0.0, -0.0003392475335051166)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"OSGB 1936\",DATUM[\"OSGB_1936\",SPHEROID[\"Airy 1830\",6377563.396,299.3249646000044,AUTHORITY[\"EPSG\",\"7001\"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY[\"EPSG\",\"6277\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4277\"]]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# <EUDEMv11> Transform DEM from ETRS-89 GCS to OSGB-36 GCS.\n",
    "\n",
    "print('\\n>>> <EUDEMv11> Transform DEM from ETRS-89 GCS to OSGB-36 GCS.')\n",
    "\n",
    "path = PATH_EUDEM\n",
    "dem_in = path + EUDEM_GCS_EU\n",
    "dem_out = path + EUDEM_GCS_UK\n",
    "\n",
    "epsg_in = EPSG_ETRS89\n",
    "epsg_out = EPSG_OSGB36\n",
    "\n",
    "transprojcnvt_dem(dem_in, epsg_in, dem_out, epsg_out)\n",
    "\n",
    "data = gdal.Open(dem_in)\n",
    "get_dem_info(data, if_print=True)\n",
    "data = gdal.Open(dem_out)\n",
    "get_dem_info(data, if_print=True)\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> Read London air quality monitoring station data file.\n",
      "  Unnamed: 0 api_data need_prediction  historical_data   Latitude  Longitude  \\\n",
      "0        BX9     True             NaN             True  51.465983   0.184877   \n",
      "1        BX1     True             NaN             True  51.465983   0.184877   \n",
      "2        BL0     True            True             True  51.522287  -0.125848   \n",
      "\n",
      "           SiteType                   SiteName  \n",
      "0          Suburban  Bexley - Slade Green FDMS  \n",
      "1          Suburban       Bexley - Slade Green  \n",
      "2  Urban Background        Camden - Bloomsbury  \n",
      "\n",
      "*==> The number of stations is: 24\n",
      "\n",
      "*==> The location (lat, lng) of stations are:\n",
      " [[51.46598327  0.18487713]\n",
      " [51.46598327  0.18487713]\n",
      " [51.522287   -0.125848  ]\n",
      " [51.52770662 -0.12905321]\n",
      " [51.544219   -0.175284  ]\n",
      " [51.51452534 -0.10451563]\n",
      " [51.51384718 -0.07776568]\n",
      " [51.410039   -0.127523  ]\n",
      " [51.490532    0.074003  ]\n",
      " [51.45258     0.070766  ]\n",
      " [51.486957    0.095111  ]\n",
      " [51.456357    0.040725  ]\n",
      " [51.4563      0.085606  ]\n",
      " [51.617327   -0.298775  ]\n",
      " [51.52078746  0.20546071]\n",
      " [51.48878    -0.441627  ]\n",
      " [51.52104675 -0.21349214]\n",
      " [51.52104675 -0.21349214]\n",
      " [51.474954   -0.039641  ]\n",
      " [51.56948433  0.08290747]\n",
      " [51.42525604 -0.34560829]\n",
      " [51.3892869  -0.14166153]\n",
      " [51.51504617 -0.00841849]\n",
      " [51.52254    -0.15459   ]]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Read London air quality monitoring station data file.\n",
    "\n",
    "print('\\n>>> Read London air quality monitoring station data file.')\n",
    "\n",
    "site_data = read_csv(PATH_LD_STATION_DATA)\n",
    "print(site_data.head(3))\n",
    "\n",
    "site_num = site_data['SiteName'].count()\n",
    "print('\\n*==> The number of stations is: %d' % site_num)\n",
    "\n",
    "# Get the latitude and longitude of stations.\n",
    "site_latlng = np.zeros((site_num, 2))\n",
    "site_latlng[:, 0] = site_data['Latitude']  # The 0-th column - Latitude.\n",
    "site_latlng[:, 1] = site_data['Longitude']  # The 1-th column - Longitude.\n",
    "\n",
    "np.set_printoptions(suppress=True)  # Print numbers without scientific notation.\n",
    "print('\\n*==> The location (lat, lng) of stations are:\\n', site_latlng)\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> <EUDEMv11> Get the elevation from DEM in WGS-84 GCS.\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30505\n",
      "The number of column (width) is: 52428\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]\n",
      "\n",
      "The 6 GeoTransform parameters of DEM are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "\n",
      "*==> The elevation information of stations is:\n",
      " [[    0.            51.46598327     0.18487713 22174.\n",
      "  36706.            12.63551998]\n",
      " [    1.            51.46598327     0.18487713 22174.\n",
      "  36706.            12.63551998]\n",
      " [    2.            51.522287      -0.125848   22008.\n",
      "  35790.            37.60099411]\n",
      " [    3.            51.52770662    -0.12905321 21992.\n",
      "  35781.            30.68862724]\n",
      " [    4.            51.544219      -0.175284   21943.\n",
      "  35645.            61.01919174]\n",
      " [    5.            51.51452534    -0.10451563 22031.\n",
      "  35853.            25.66654587]\n",
      " [    6.            51.51384718    -0.07776568 22033.\n",
      "  35932.            29.40543556]\n",
      " [    7.            51.410039      -0.127523   22339.\n",
      "  35785.            36.34508514]\n",
      " [    8.            51.490532       0.074003   22102.\n",
      "  36379.            11.27258968]\n",
      " [    9.            51.45258        0.070766   22214.\n",
      "  36370.            66.39667511]\n",
      " [   10.            51.486957       0.095111   22112.\n",
      "  36442.            11.88549519]\n",
      " [   11.            51.456357       0.040725   22202.\n",
      "  36281.            31.31890869]\n",
      " [   12.            51.4563         0.085606   22203.\n",
      "  36414.            64.95016479]\n",
      " [   13.            51.617327      -0.298775   21728.\n",
      "  35281.            79.23999786]\n",
      " [   14.            51.52078746     0.20546071 22012.\n",
      "  36767.             7.36576653]\n",
      " [   15.            51.48878       -0.441627   22107.\n",
      "  34859.            27.24001312]\n",
      " [   16.            51.52104675    -0.21349214 22012.\n",
      "  35532.            25.52720642]\n",
      " [   17.            51.52104675    -0.21349214 22012.\n",
      "  35532.            25.52720642]\n",
      " [   18.            51.474954      -0.039641   22148.\n",
      "  36044.            14.32332325]\n",
      " [   19.            51.56948433     0.08290747 21869.\n",
      "  36406.            15.99652195]\n",
      " [   20.            51.42525604    -0.34560829 22294.\n",
      "  35142.            13.35388565]\n",
      " [   21.            51.3892869     -0.14166153 22400.\n",
      "  35744.            32.15060806]\n",
      " [   22.            51.51504617    -0.00841849 22029.\n",
      "  36136.             5.19105673]\n",
      " [   23.            51.52254       -0.15459    22007.\n",
      "  35706.            36.85784912]]\n",
      "\n",
      "*==> The elevation value of stations is:\n",
      " [12 12 37 30 61 25 29 36 11 66 11 31 64 79  7 27 25 25 14 15 13 32  5 36]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# <EUDEMv11> Get the elevation from DEM in WGS-84 GCS.\n",
    "\n",
    "print('\\n>>> <EUDEMv11> Get the elevation from DEM in WGS-84 GCS.')\n",
    "\n",
    "dem_gcs = gdal.Open(PATH_EUDEM + EUDEM_GCS_WD)\n",
    "get_dem_info(dem_gcs, if_print=True)\n",
    "\n",
    "site_ele_eudem_wgs = get_elevation(dem_gcs, site_latlng)\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "print('\\n*==> The elevation information of stations is:\\n', site_ele_eudem_wgs)\n",
    "print('\\n*==> The elevation value of stations is:\\n', site_ele_eudem_wgs[:, 5].astype(int))\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> <EUDEMv11> Get the elevation from DEM in OSGB-36 GCS.\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30506\n",
      "The number of column (width) is: 52430\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266102928602475, 0.0003392475335051166, 0.0, 58.98812292836723, 0.0, -0.0003392475335051166)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"OSGB 1936\",DATUM[\"OSGB_1936\",SPHEROID[\"Airy 1830\",6377563.396,299.3249646000044,AUTHORITY[\"EPSG\",\"7001\"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY[\"EPSG\",\"6277\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4277\"]]\n",
      "\n",
      "The 6 GeoTransform parameters of DEM are:\n",
      " (-12.266102928602475, 0.0003392475335051166, 0.0, 58.98812292836723, 0.0, -0.0003392475335051166)\n",
      "\n",
      "*==> The elevation information of stations is:\n",
      " [[    0.            51.46598327     0.18487713 22173.\n",
      "  36702.            12.47727203]\n",
      " [    1.            51.46598327     0.18487713 22173.\n",
      "  36702.            12.47727203]\n",
      " [    2.            51.522287      -0.125848   22007.\n",
      "  35786.            37.75251389]\n",
      " [    3.            51.52770662    -0.12905321 21991.\n",
      "  35776.            26.8823452 ]\n",
      " [    4.            51.544219      -0.175284   21942.\n",
      "  35640.            58.50819397]\n",
      " [    5.            51.51452534    -0.10451563 22030.\n",
      "  35849.            30.98666763]\n",
      " [    6.            51.51384718    -0.07776568 22032.\n",
      "  35928.            32.92979813]\n",
      " [    7.            51.410039      -0.127523   22338.\n",
      "  35781.            35.5172081 ]\n",
      " [    8.            51.490532       0.074003   22101.\n",
      "  36375.            11.85224247]\n",
      " [    9.            51.45258        0.070766   22213.\n",
      "  36365.            71.22161865]\n",
      " [   10.            51.486957       0.095111   22111.\n",
      "  36437.             8.77665329]\n",
      " [   11.            51.456357       0.040725   22201.\n",
      "  36277.            30.45423889]\n",
      " [   12.            51.4563         0.085606   22202.\n",
      "  36409.            66.69743347]\n",
      " [   13.            51.617327      -0.298775   21727.\n",
      "  35276.            78.91999817]\n",
      " [   14.            51.52078746     0.20546071 22011.\n",
      "  36762.             8.10325527]\n",
      " [   15.            51.48878       -0.441627   22106.\n",
      "  34855.            26.04798698]\n",
      " [   16.            51.52104675    -0.21349214 22011.\n",
      "  35527.            24.9171505 ]\n",
      " [   17.            51.52104675    -0.21349214 22011.\n",
      "  35527.            24.9171505 ]\n",
      " [   18.            51.474954      -0.039641   22147.\n",
      "  36040.             8.97394276]\n",
      " [   19.            51.56948433     0.08290747 21868.\n",
      "  36401.            16.96926498]\n",
      " [   20.            51.42525604    -0.34560829 22293.\n",
      "  35138.            13.82033825]\n",
      " [   21.            51.3892869     -0.14166153 22399.\n",
      "  35739.            32.94750977]\n",
      " [   22.            51.51504617    -0.00841849 22028.\n",
      "  36132.             9.25135612]\n",
      " [   23.            51.52254       -0.15459    22006.\n",
      "  35701.            35.54499435]]\n",
      "\n",
      "*==> The elevation value of stations is:\n",
      " [12 12 37 26 58 30 32 35 11 71  8 30 66 78  8 26 24 24  8 16 13 32  9 35]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# <EUDEMv11> Get the elevation from DEM in OSGB-36 GCS.\n",
    "\n",
    "print('\\n>>> <EUDEMv11> Get the elevation from DEM in OSGB-36 GCS.')\n",
    "\n",
    "dem_gcs = gdal.Open(PATH_EUDEM + EUDEM_GCS_UK)\n",
    "get_dem_info(dem_gcs, if_print=True)\n",
    "\n",
    "site_ele_eudem_osgb = get_elevation(dem_gcs, site_latlng)\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "print('\\n*==> The elevation information of stations is:\\n', site_ele_eudem_osgb)\n",
    "print('\\n*==> The elevation value of stations is:\\n', site_ele_eudem_osgb[:, 5].astype(int))\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> <EUDEMv11> Get the elevation from DEM in ETRS-89 GCS.\n",
      "\n",
      "The information of DEM:\n",
      "The number of row (height) is: 30505\n",
      "The number of column (width) is: 52428\n",
      "The number of band is: 1\n",
      "The 6 GeoTransform parameters are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "The GCS/PCS information is:\n",
      " GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4258\"]]\n",
      "\n",
      "The 6 GeoTransform parameters of DEM are:\n",
      " (-12.266271034206424, 0.00033920980218715875, 0.0, 58.98761768429964, 0.0, -0.00033920980218715875)\n",
      "\n",
      "*==> The elevation information of stations is:\n",
      " [[    0.            51.46598327     0.18487713 22174.\n",
      "  36706.            12.63551998]\n",
      " [    1.            51.46598327     0.18487713 22174.\n",
      "  36706.            12.63551998]\n",
      " [    2.            51.522287      -0.125848   22008.\n",
      "  35790.            37.60099411]\n",
      " [    3.            51.52770662    -0.12905321 21992.\n",
      "  35781.            30.68862724]\n",
      " [    4.            51.544219      -0.175284   21943.\n",
      "  35645.            61.01919174]\n",
      " [    5.            51.51452534    -0.10451563 22031.\n",
      "  35853.            25.66654587]\n",
      " [    6.            51.51384718    -0.07776568 22033.\n",
      "  35932.            29.40543556]\n",
      " [    7.            51.410039      -0.127523   22339.\n",
      "  35785.            36.34508514]\n",
      " [    8.            51.490532       0.074003   22102.\n",
      "  36379.            11.27258968]\n",
      " [    9.            51.45258        0.070766   22214.\n",
      "  36370.            66.39667511]\n",
      " [   10.            51.486957       0.095111   22112.\n",
      "  36442.            11.88549519]\n",
      " [   11.            51.456357       0.040725   22202.\n",
      "  36281.            31.31890869]\n",
      " [   12.            51.4563         0.085606   22203.\n",
      "  36414.            64.95016479]\n",
      " [   13.            51.617327      -0.298775   21728.\n",
      "  35281.            79.23999786]\n",
      " [   14.            51.52078746     0.20546071 22012.\n",
      "  36767.             7.36576653]\n",
      " [   15.            51.48878       -0.441627   22107.\n",
      "  34859.            27.24001312]\n",
      " [   16.            51.52104675    -0.21349214 22012.\n",
      "  35532.            25.52720642]\n",
      " [   17.            51.52104675    -0.21349214 22012.\n",
      "  35532.            25.52720642]\n",
      " [   18.            51.474954      -0.039641   22148.\n",
      "  36044.            14.32332325]\n",
      " [   19.            51.56948433     0.08290747 21869.\n",
      "  36406.            15.99652195]\n",
      " [   20.            51.42525604    -0.34560829 22294.\n",
      "  35142.            13.35388565]\n",
      " [   21.            51.3892869     -0.14166153 22400.\n",
      "  35744.            32.15060806]\n",
      " [   22.            51.51504617    -0.00841849 22029.\n",
      "  36136.             5.19105673]\n",
      " [   23.            51.52254       -0.15459    22007.\n",
      "  35706.            36.85784912]]\n",
      "\n",
      "*==> The elevation value of stations is:\n",
      " [12 12 37 30 61 25 29 36 11 66 11 31 64 79  7 27 25 25 14 15 13 32  5 36]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# <EUDEMv11> Get the elevation from DEM in ETRS-89 GCS.\n",
    "\n",
    "print('\\n>>> <EUDEMv11> Get the elevation from DEM in ETRS-89 GCS.')\n",
    "\n",
    "dem_gcs = gdal.Open(PATH_EUDEM + EUDEM_GCS_EU)\n",
    "get_dem_info(dem_gcs, if_print=True)\n",
    "\n",
    "site_ele_eudem_etrs = get_elevation(dem_gcs, site_latlng)\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "print('\\n*==> The elevation information of stations is:\\n', site_ele_eudem_etrs)\n",
    "print('\\n*==> The elevation value of stations is:\\n', site_ele_eudem_etrs[:, 5].astype(int))\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      ">>> Compare the elevation obtained from different DEMs.\n",
      "\n",
      "EUDEMv11_WD [12 12 37 30 61 25 29 36 11 66 11 31 64 79  7 27 25 25 14 15 13 32  5 36]\n",
      "\n",
      "EUDEMv11_UK [12 12 37 26 58 30 32 35 11 71  8 30 66 78  8 26 24 24  8 16 13 32  9 35]\n",
      "\n",
      "EUDEMv11_EU [12 12 37 30 61 25 29 36 11 66 11 31 64 79  7 27 25 25 14 15 13 32  5 36]\n",
      "\n",
      ">>> Complete!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Compare the elevation obtained from different DEMs.\n",
    "\n",
    "print('\\n>>> Compare the elevation obtained from different DEMs.')\n",
    "\n",
    "print('\\nEUDEMv11_WD', site_ele_eudem_wgs[:, 5].astype(int))\n",
    "print('\\nEUDEMv11_UK', site_ele_eudem_osgb[:, 5].astype(int))\n",
    "print('\\nEUDEMv11_EU', site_ele_eudem_etrs[:, 5].astype(int))\n",
    "\n",
    "print('\\n>>> Complete!\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}