Copernicus - retrieve elevation for any point in Europe (Pandas version).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieve elevation for any point in Europe\n",
    "- The dataset is the [European Digital Elevation Model (EU-DEM), version 1.1](https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1).\n",
    "- The source is the European Union's Copernicus Programme. More info [here](https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1?tab=metadata).\n",
    "- The dataset is in [Raster format](https://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/what-is-raster-data.htm) (geotiff), divided into 1000x1000km tiles, 25m resolution, accuracy of [(~7m RMSE)](https://ec.europa.eu/eurostat/documents/7116161/7172326/Report-EU-DEM-statistical-validation-August2014.pdf).\n",
    "- The only dependencies of this script are Numpy, Pandas, Pyproj and Rasterio.\n",
    "- The process was designed to minimize memory consumption. This means that pieces, tiles or chunks of the rasters are read sequentially from storage.\n",
    "- A distributed version of this script can be easily created with Dask for Big Data use cases.\n",
    "- Source: https://github.com/carlosg-m/goose-lab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import rasterio\n",
    "import pyproj\n",
    "\n",
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Function to create regular grid, and function to intersect points with regular grid\n",
    "- Efficient and simple, built in Numpy.\n",
    "- Compatible with Rasterio windows.\n",
    "- No need for spatial index or brute force queries.\n",
    "- Shapely, Pygeos or handling geometries is also not required.\n",
    "- More information: [Grid creation and intersection notebooks.](https://github.com/carlosg-m/goose-lab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_grid(xmin, ymin, xmax, ymax, xstep, ystep):\n",
    "    \"\"\"\n",
    "    Create bounding coordinates for a regular 2d grid with numpy.\n",
    "    These coordinates can then be used with shapely, pygeos and geopandas.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Generate two evenly spaced ranges\n",
    "    x = np.linspace(xmin, xmax, int(xstep+1))\n",
    "    y = np.linspace(ymin, ymax, int(ystep+1))\n",
    "    \n",
    "    # Get center points for each range\n",
    "    x = np.mean((x[1:], x[:-1]), axis=0)\n",
    "    y = np.mean((y[1:], y[:-1]), axis=0)\n",
    "    \n",
    "    # Cartesian product of both ranges \n",
    "    x, y = np.meshgrid(x, y)\n",
    "    x, y = x.ravel(), y.ravel()\n",
    "\n",
    "    # Calculate delta x and delta y\n",
    "    xd = (xmax - xmin) / (xstep * 2)\n",
    "    yd = (ymax - ymin) / (ystep * 2)\n",
    "    \n",
    "    # Repeat center points and add respective delta\n",
    "    grid = np.column_stack((x,y))\n",
    "    grid = np.tile(grid, 2) + np.array([[-xd, -yd, xd, yd]])\n",
    "    \n",
    "    # Generate index\n",
    "    grid_index = np.arange(0, xstep * ystep)\n",
    "    \n",
    "    # Output coordinates: xmin, ymin, xmax, ymax\n",
    "    return grid, grid_index\n",
    "\n",
    "def intersect_points_grid(xy, xmin, ymin, xmax, ymax, xstep, ystep, include_xmax=False, include_ymax=False):\n",
    "    \"\"\"\n",
    "    Intersect 2d Points with a regular 2d Grid using NumPy.\n",
    "    Returns the index of the Grid rectangle that each Point belongs to.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Calculate delta x and delta y\n",
    "    xd = (xmax - xmin) / xstep\n",
    "    yd = (ymax - ymin) / ystep\n",
    "    \n",
    "    # Center points at zero\n",
    "    xy = xy - np.array([[xmin, ymin]])\n",
    "    \n",
    "    # Get index of x and y coordinates\n",
    "    xi = np.floor(xy[:,0] / xd) \n",
    "    yi = np.floor(xy[:,1] / yd)\n",
    "    \n",
    "    # Cast x and y index to 64 bit integer\n",
    "    xi = xi.astype(np.int64)\n",
    "    yi = yi.astype(np.int64)\n",
    "    \n",
    "    # Opened interval for x\n",
    "    if include_xmax == True:\n",
    "        xi = xi - (xy[:,0] % xd == 0)\n",
    "    \n",
    "    # Opened interval for y\n",
    "    if include_ymax == True:\n",
    "        yi = yi - (xy[:,1] % yd == 0)\n",
    "        \n",
    "    # Mask out-of-bound indexes on the x and y axis\n",
    "    xi = np.ma.array(xi, mask=((xi < 0) | (xi >= xstep)))\n",
    "    yi = np.ma.array(yi, mask=((yi < 0) | (yi >= ystep)))\n",
    "        \n",
    "    # Calculate index\n",
    "    index = xi + (yi * xstep)\n",
    "    \n",
    "    # Unmask NumPy array\n",
    "    index = index.filled(np.iinfo(np.int64).min)\n",
    "    \n",
    "    return index\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Logic to manage the EUR-DEM v1.1 dataset and to retrieve elevation for a dataframe of coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CopernicusDEM:\n",
    "    \"\"\"\n",
    "    Extract elevation values for a Pandas dataframe of WGS84 coordinates,\n",
    "    through raster files in Geotiff format provided by Copernicus programme.\n",
    "    \n",
    "    Since the rasters are of considerable size, 1000x1000km at 25 meter resolution,\n",
    "    Rasterio's reading windows are used to process chunks individually, \n",
    "    limiting memory consuption.\n",
    "    \n",
    "    The process is completely implemented in Pandas, Numpy, Pyproj and Rasterio.\n",
    "    \n",
    "    - Source: Copernicus programme\n",
    "    - Resource: European Digital Elevation Model (EU-DEM), version 1.1\n",
    "    - URL: https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, raster_paths, grid_xstep=50, grid_ystep=50, elevation_alias='elevation'):\n",
    "        \"\"\"\n",
    "        raster_paths: list of paths to the EU-DEM Geotiff raster files.\n",
    "                      ['eu_dem_v11_E00N00.TIF', 'eu_dem_v11_E00N00.TIF'].\n",
    "        \n",
    "        grid_xstep: number of chunks in the x-axis used to create the reading windows. \n",
    "                    Must be diviseble by x-range. Affects performance and memory.\n",
    "                    Default value is 50.\n",
    "        \n",
    "        grid_ystep: number of chunks in the y-axis used to create the reading windows. \n",
    "                    Must be diviseble by y-range. Affects performance and memory.\n",
    "                    Default value is 50.\n",
    "                    \n",
    "        elevation_alias: name of output column with elevation values appended to the dataframe.         \n",
    "        \"\"\"\n",
    "        \n",
    "        # ETRS-LAEA grid parameters\n",
    "        self.laea_xmin = 0\n",
    "        self.laea_ymin = 0\n",
    "        self.laea_xmax = 8000000\n",
    "        self.laea_ymax = 6000000\n",
    "        self.laea_xstep = 8\n",
    "        self.laea_ystep = 6\n",
    "        self.laea_include_xmax = False\n",
    "        self.laea_include_ymax = True\n",
    "        \n",
    "        # User parameters and properties\n",
    "        self.raster_paths = raster_paths\n",
    "        self.grid_xstep = grid_xstep\n",
    "        self.grid_ystep = grid_ystep\n",
    "        \n",
    "        # Input validations\n",
    "        self._validate_grid_xstep()\n",
    "        self._validate_grid_ystep()\n",
    "        \n",
    "        # Properties\n",
    "        self.raster_crs = self._get_raster_crs()\n",
    "        \n",
    "        # Columns\n",
    "        self.col_path = '_path'\n",
    "        self.col_laea_tile = '_etrs_laea_tile'\n",
    "        self.col_grid_index = '_grid_index'\n",
    "        self.col_temp_key = '_temp_key'\n",
    "        self.col_grid_xmin = '_grid_xmin'\n",
    "        self.col_grid_ymin = '_grid_ymin'\n",
    "        self.col_grid_xmax = '_grid_xmax'\n",
    "        self.col_grid_ymax = '_grid_ymax'\n",
    "        self.col_proj_x = '_x'\n",
    "        self.col_proj_y = '_y'\n",
    "        self.col_elevation = elevation_alias\n",
    "                    \n",
    "    def _validate_grid_xstep(self):\n",
    "        \"\"\"\n",
    "        Check if xstep parameter is valid.\n",
    "        \"\"\" \n",
    "        \n",
    "        x_range = self.laea_xmax - self.laea_xmin\n",
    "        x_step = self.laea_xstep * self.grid_xstep\n",
    "        \n",
    "        is_valid = x_range % x_step == 0\n",
    "        assert_message = 'Parameter \"grid_xstep\" must be divisible by xrange ({}).'.format(x_range)\n",
    "        \n",
    "        assert is_valid, assert_message\n",
    "        \n",
    "    def _validate_grid_ystep(self):\n",
    "        \"\"\"\n",
    "        Check if ystep parameter is valid.\n",
    "        \"\"\" \n",
    "        \n",
    "        y_range = self.laea_ymax - self.laea_ymin\n",
    "        y_step = self.laea_ystep * self.grid_ystep\n",
    "        \n",
    "        is_valid = y_range % y_step == 0\n",
    "        assert_message = 'Parameter \"grid_ystep\" must be divisible by yrange ({}).'.format(y_range)\n",
    "        \n",
    "        assert is_valid, assert_message\n",
    "    \n",
    "    def _get_raster_crs(self):\n",
    "        \"\"\"\n",
    "        Read first geotiff raster file and retrieve Coordinate Reference System (CRS).\n",
    "        \"\"\"\n",
    "        \n",
    "        first_raster = self.raster_paths[0]\n",
    "        \n",
    "        with rasterio.open(first_raster) as src:\n",
    "            crs = src.crs\n",
    "        \n",
    "        return crs\n",
    "    \n",
    "    def _get_etrs_laea_tile(self, df, xmin, ymin):\n",
    "        \"\"\"\n",
    "        Convert rectangle bounds into ETRS-LAEA tile denomination (\"E30N30\").\n",
    "        The first digit, at the kilometer level, of the xmin bound is East and of the ymin bound is North.\n",
    "        \"\"\"\n",
    "        \n",
    "        etrs_e = np.floor(df[xmin] / 10**6)\n",
    "        etrs_e = (etrs_e * 10).astype(np.int64)\n",
    "        etrs_e = 'E' + etrs_e.astype('str').str.zfill(2)\n",
    "        \n",
    "        etrs_n = np.floor(df[ymin] / 10**6)\n",
    "        etrs_n = (etrs_n * 10).astype(np.int64)\n",
    "        etrs_n = 'N' + etrs_n.astype('str').str.zfill(2)\n",
    "        \n",
    "        df[self.col_laea_tile] = etrs_e + etrs_n\n",
    "        \n",
    "        return df\n",
    "        \n",
    "    def _read_raster_metadata(self):\n",
    "        \"\"\"\n",
    "        Read the metadata of every raster in the \"raster_paths\" parameter.\n",
    "        Only the bounds of the raster and filename is kept.\n",
    "        The bounds are then converted to the ETRS-LAEA tile denomination and dropped.\n",
    "        \"\"\"\n",
    "        \n",
    "        df = []\n",
    "        for path in self.raster_paths:\n",
    "            with rasterio.open(path) as src:\n",
    "                bounds = dict(zip(('_raster_xmin', '_raster_ymin', '_raster_xmax', '_raster_ymax'), src.bounds))\n",
    "                bounds[self.col_path] = path\n",
    "                df.append(bounds)\n",
    "\n",
    "        df = pd.DataFrame(df)\n",
    "        df = self._get_etrs_laea_tile(df, xmin='_raster_xmin', ymin='_raster_ymin')\n",
    "        df = df.drop(['_raster_xmin', '_raster_ymin', '_raster_xmax', '_raster_ymax'], axis=1)\n",
    "        \n",
    "        return df\n",
    "    \n",
    "    def _create_reading_grid(self):\n",
    "        \"\"\"\n",
    "        Create the reading grid that will be the base for retrieving raster tiles.\n",
    "        This grid must coincide and contain the original 1000x1000km ETRS-LAEA grid, and can \n",
    "        be subdivided into finer tiles, according to the \"grid_xstep\" and \"grid_ystep\" parameters.\n",
    "        The ETRS-LAEA tile denomination is also added.\n",
    "        Each tile of the created grid will have an index, represented by an intenger.\n",
    "        \"\"\"\n",
    "        \n",
    "        raster_grid, raster_grid_index = create_grid(xmin=self.laea_xmin, \n",
    "                                                     ymin=self.laea_ymin, \n",
    "                                                     xmax=self.laea_xmax, \n",
    "                                                     ymax=self.laea_ymax, \n",
    "                                                     xstep=self.laea_xstep*self.grid_xstep, \n",
    "                                                     ystep=self.laea_ystep*self.grid_ystep)\n",
    "        \n",
    "        df = pd.DataFrame(raster_grid, columns=[self.col_grid_xmin,self.col_grid_ymin,self.col_grid_xmax,self.col_grid_ymax])\n",
    "        df[self.col_grid_index] = raster_grid_index\n",
    "        df = self._get_etrs_laea_tile(df, xmin=self.col_grid_xmin, ymin=self.col_grid_ymin)\n",
    "        \n",
    "        return df\n",
    "    \n",
    "    def _get_grid_scope(self):\n",
    "        \"\"\"\n",
    "        Grid scope is the intersection of the rasters included in input parameters and the reading grid.\n",
    "        The ETRS-LAEA tile denomination is used as key in a database join.\n",
    "        \"\"\"\n",
    "        \n",
    "        df = self._create_reading_grid().merge(self._read_raster_metadata(), 'inner', self.col_laea_tile)\n",
    "        \n",
    "        df = df.drop(self.col_laea_tile, axis=1)\n",
    "        \n",
    "        return df\n",
    "    \n",
    "    def _project_points(self, df, lat_col, lon_col, input_crs):\n",
    "        \"\"\"\n",
    "        Projects WGS84 coordinates into the raster CRS (EPSG:3035).\n",
    "        The original coordinates are droped.\n",
    "        \"\"\"\n",
    "        \n",
    "        projection = pyproj.Transformer.from_crs(input_crs, self.raster_crs)\n",
    "        \n",
    "        df[[self.col_proj_x, self.col_proj_y]] = np.column_stack(projection.transform(df[lat_col].values, df[lon_col].values))\n",
    "        \n",
    "        df = df.drop([lat_col, lon_col], axis=1)\n",
    "        \n",
    "        return df\n",
    "    \n",
    "    def _get_grid_index(self, df):\n",
    "        \"\"\"\n",
    "        Converts the projected points (x, y), into the respective reading grid index they belong to.\n",
    "        \"\"\"\n",
    "        \n",
    "        df[self.col_grid_index] = intersect_points_grid(xy=df[[self.col_proj_x, self.col_proj_y]].values, \n",
    "                                                        xmin=self.laea_xmin, \n",
    "                                                        ymin=self.laea_ymin, \n",
    "                                                        xmax=self.laea_xmax, \n",
    "                                                        ymax=self.laea_ymax, \n",
    "                                                        xstep=self.laea_xstep*self.grid_xstep, \n",
    "                                                        ystep=self.laea_ystep*self.grid_ystep,\n",
    "                                                        include_xmax=self.laea_include_xmax,\n",
    "                                                        include_ymax=self.laea_include_ymax)\n",
    "        \n",
    "        return df\n",
    "    \n",
    "    def _generate_temp_key(self, df):\n",
    "        \"\"\"\n",
    "        Generates a temporary key that serves as reference \n",
    "        for the input dataset through the whole process.\n",
    "        \"\"\"\n",
    "        \n",
    "        df[self.col_temp_key] = range(len(df))\n",
    "        \n",
    "        return df\n",
    "        \n",
    "    def _get_elevation_pandas(self, df):\n",
    "        \"\"\"\n",
    "        Function to be applied to partitions of the Pandas dataframe.\n",
    "        It uses rasterio and all the predefined parameters to retrieve the elevation \n",
    "        of all points contained in the partition.\n",
    "        The partitions are defined by grid index. This means that each window is only read once.\n",
    "        The auxiliary columns that are no longer required are dropped to save memory.\n",
    "        \"\"\"\n",
    "        \n",
    "        path = df[self.col_path].iloc[0]\n",
    "        xmin = df[self.col_grid_xmin].iloc[0]\n",
    "        ymin = df[self.col_grid_ymin].iloc[0]\n",
    "        xmax = df[self.col_grid_xmax].iloc[0]\n",
    "        ymax = df[self.col_grid_ymax].iloc[0]\n",
    "        \n",
    "        cols_drop = [self.col_grid_index, \n",
    "                     self.col_path, \n",
    "                     self.col_grid_xmin, \n",
    "                     self.col_grid_ymin, \n",
    "                     self.col_grid_xmax, \n",
    "                     self.col_grid_ymax]\n",
    "        \n",
    "        df = df.drop(cols_drop, axis=1)\n",
    "        \n",
    "        with rasterio.open(path) as src:\n",
    "\n",
    "            raster_affine = src.transform\n",
    "            rasterio_window = rasterio.windows.from_bounds(left=xmin, \n",
    "                                                           bottom=ymin, \n",
    "                                                           right=xmax, \n",
    "                                                           top=ymax, \n",
    "                                                           transform=raster_affine)\n",
    "\n",
    "            window_affine = src.window_transform(rasterio_window)\n",
    "            index = rasterio.transform.rowcol(window_affine, df[self.col_proj_x], df[self.col_proj_y])\n",
    "            df = df.drop([self.col_proj_x, self.col_proj_y], axis=1)\n",
    "            \n",
    "            elevation = src.read(1, window=rasterio_window)\n",
    "            elevation = elevation[index]\n",
    "            df[self.col_elevation] = elevation\n",
    "\n",
    "            return df\n",
    "    \n",
    "    def _get_df_elevation(self, df, lat_col, lon_col, input_crs):\n",
    "        \"\"\"\n",
    "        Base process to retrieve elevation for valid points.\n",
    "        The reading grid with only tiles in scope is retrieved.\n",
    "        The original df is copied, maintaining only the required columns (temp_id, lat, lon).\n",
    "        Null latitude and longitude coordinates are filtered.\n",
    "        The coordinates are then projected, and the grid index is obtained.\n",
    "        The dataframe is filtered only to include points that fall within the bounds of the input rasters.\n",
    "        The reading metadata, path and window bounds, are appended to the dataframe.\n",
    "        Finaly the dataset is grouped or partioned through grid index and the elvation is obtained.\n",
    "        \"\"\"\n",
    "        \n",
    "        # Get grid tiles that are within scope of available input raster files.\n",
    "        grid = self._get_grid_scope()\n",
    "        \n",
    "        # Copy original dataframe and retrieve necessary columns.\n",
    "        # Drop null coordinates.\n",
    "        df = df[[self.col_temp_key, lat_col, lon_col]].copy()\n",
    "        df = df.dropna(subset=[lat_col, lon_col])\n",
    "        \n",
    "        # Project coordinates.\n",
    "        # Get grid index.\n",
    "        df = self._project_points(df, lat_col, lon_col, input_crs)\n",
    "        df = self._get_grid_index(df)\n",
    "        \n",
    "        # Filter points that fall ouside of available grid bounds.\n",
    "        # Join reading metadata.\n",
    "        df = df.merge(grid[[self.col_grid_index]], 'inner', self.col_grid_index)\n",
    "        df = df.merge(grid, 'inner', self.col_grid_index)\n",
    "        \n",
    "        # Group by grid index to obtain partitions.\n",
    "        # Apply rasterio based function to get elevation.\n",
    "        df = df.groupby(self.col_grid_index)\n",
    "        df = df.apply(self._get_elevation_pandas)\n",
    "        \n",
    "        return df\n",
    "    \n",
    "    def get_elevation(self, df, lat_col, lon_col, input_crs='epsg:4326'):\n",
    "        \"\"\"\n",
    "        User interface. Allows to retrieve elevation for a dataframe that contains WGS84 coordinates.\n",
    "        The elevation column is appended to the input dataframe.\n",
    "        All original columns of the dataframe are also returned (e.g. keys).\n",
    "        Null input values values are supported and kept. \n",
    "        Elevation values for points that are out-of-bounds are treated as null.\n",
    "        The raster uses a mask, a large negative value, around water bodies and the ocean.\n",
    "        \n",
    "        df: Pandas dataframe with coordinates and miscellaneous columns.\n",
    "        lat_col: name of column that contains latitude.\n",
    "        lon_col: name of columns that contains longitude.\n",
    "        input_crs: input CRS, by default is WGS84 (EPSG:4326).\n",
    "        \"\"\"\n",
    "        \n",
    "        df = self._generate_temp_key(df)\n",
    "        df = df.merge(self._get_df_elevation(df, lat_col, lon_col, input_crs), 'left', self.col_temp_key)\n",
    "        df = df.drop(self.col_temp_key, axis=1)\n",
    "        \n",
    "        return df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example 1: Generate WGS84 test points over Portugal and Spain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_test_points(sample=3000):\n",
    "    \"\"\"\n",
    "    Generate test points over Portugal and Spain.\n",
    "    \"\"\"\n",
    "    \n",
    "    lat = np.linspace(36, 43, sample+1)\n",
    "    lon = np.linspace(-11, -4, sample*2+1)\n",
    "    \n",
    "    lat, lon = np.meshgrid(lat, lon)\n",
    "\n",
    "    lat = lat.ravel()\n",
    "    lon = lon.ravel()\n",
    "\n",
    "    df = np.column_stack((lat, lon))\n",
    "    df = pd.DataFrame(df, columns=['lat','lon'])\n",
    "    df['some_key'] = np.arange(len(df))\n",
    "    \n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of examples: 18009001\n",
      "Dataframe memory usage (MB): 432.216152\n"
     ]
    }
   ],
   "source": [
    "test_points = get_test_points()\n",
    "\n",
    "print('Number of examples:', len(test_points))\n",
    "print('Dataframe memory usage (MB):', test_points.memory_usage(index=True, deep=True).sum() / 10**6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>lat</th>\n",
       "      <th>lon</th>\n",
       "      <th>some_key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>36.000000</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>36.002333</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>36.004667</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>36.007000</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>36.009333</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         lat   lon  some_key\n",
       "0  36.000000 -11.0         0\n",
       "1  36.002333 -11.0         1\n",
       "2  36.004667 -11.0         2\n",
       "3  36.007000 -11.0         3\n",
       "4  36.009333 -11.0         4"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_points.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Instantiate object with EUR-DEM raster paths\n",
    "- The dataset can be downloaded from the official source: https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1\n",
    "- The raster tiles don't have to be consecutive. You can download just the ones you need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "copernicus = CopernicusDEM(raster_paths=['eu_dem_v11_E20N10.TIF', 'eu_dem_v11_E20N20.TIF'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Retrieve elevation values for test dataframe\n",
    "- The elevation column is appended to the input dataframe.\n",
    "- All original columns of the dataframe are also returned (e.g. keys).\n",
    "- Null input values values are supported and kept. \n",
    "- Elevation values for points that are out-of-bounds are treated as null.\n",
    "- The raster uses a mask, a large negative value, around water bodies and the ocean.\n",
    "- Execution time and memory consumption displayed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 3971.30 MiB, increment: 3429.87 MiB\n",
      "CPU times: user 1min 19s, sys: 11.3 s, total: 1min 30s\n",
      "Wall time: 1min 30s\n"
     ]
    }
   ],
   "source": [
    "%time %memit test_points = copernicus.get_elevation(test_points, lat_col='lat', lon_col='lon')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of points: 18009001\n",
      "Number of points with elevation: 14379213\n"
     ]
    }
   ],
   "source": [
    "print('Number of points:', len(test_points))\n",
    "print('Number of points with elevation:', test_points['elevation'].count())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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>lat</th>\n",
       "      <th>lon</th>\n",
       "      <th>some_key</th>\n",
       "      <th>elevation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>36.000000</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>0</td>\n",
       "      <td>-3.402823e+38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>36.002333</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>1</td>\n",
       "      <td>-3.402823e+38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>36.004667</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>2</td>\n",
       "      <td>-3.402823e+38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>36.007000</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>3</td>\n",
       "      <td>-3.402823e+38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>36.009333</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>4</td>\n",
       "      <td>-3.402823e+38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18008996</th>\n",
       "      <td>42.990667</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>18008996</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18008997</th>\n",
       "      <td>42.993000</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>18008997</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18008998</th>\n",
       "      <td>42.995333</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>18008998</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18008999</th>\n",
       "      <td>42.997667</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>18008999</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18009000</th>\n",
       "      <td>43.000000</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>18009000</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>18009001 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                lat   lon  some_key     elevation\n",
       "0         36.000000 -11.0         0 -3.402823e+38\n",
       "1         36.002333 -11.0         1 -3.402823e+38\n",
       "2         36.004667 -11.0         2 -3.402823e+38\n",
       "3         36.007000 -11.0         3 -3.402823e+38\n",
       "4         36.009333 -11.0         4 -3.402823e+38\n",
       "...             ...   ...       ...           ...\n",
       "18008996  42.990667  -4.0  18008996           NaN\n",
       "18008997  42.993000  -4.0  18008997           NaN\n",
       "18008998  42.995333  -4.0  18008998           NaN\n",
       "18008999  42.997667  -4.0  18008999           NaN\n",
       "18009000  43.000000  -4.0  18009000           NaN\n",
       "\n",
       "[18009001 rows x 4 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example 2: Retrieve elevation values for smaller dataframe with airport information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 1485.55 MiB, increment: 24.62 MiB\n",
      "CPU times: user 802 ms, sys: 72.2 ms, total: 874 ms\n",
      "Wall time: 990 ms\n"
     ]
    }
   ],
   "source": [
    "airports = pd.DataFrame([['LPPT', 38.775600, -9.135400],\n",
    "                         ['LPPR', 41.242100, -8.678600],\n",
    "                         ['LPFR', 37.017600, -7.969700],\n",
    "                         ['LPBJ', 38.063700, -7.939200]], columns=['ICAO', 'Latitude', 'Longitude'])\n",
    "\n",
    "%time %memit airports = copernicus.get_elevation(airports, lat_col='Latitude', lon_col='Longitude')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>ICAO</th>\n",
       "      <th>Latitude</th>\n",
       "      <th>Longitude</th>\n",
       "      <th>elevation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>LPPT</td>\n",
       "      <td>38.7756</td>\n",
       "      <td>-9.1354</td>\n",
       "      <td>101.134453</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>LPPR</td>\n",
       "      <td>41.2421</td>\n",
       "      <td>-8.6786</td>\n",
       "      <td>60.038509</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LPFR</td>\n",
       "      <td>37.0176</td>\n",
       "      <td>-7.9697</td>\n",
       "      <td>5.508499</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>LPBJ</td>\n",
       "      <td>38.0637</td>\n",
       "      <td>-7.9392</td>\n",
       "      <td>186.556244</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   ICAO  Latitude  Longitude   elevation\n",
       "0  LPPT   38.7756    -9.1354  101.134453\n",
       "1  LPPR   41.2421    -8.6786   60.038509\n",
       "2  LPFR   37.0176    -7.9697    5.508499\n",
       "3  LPBJ   38.0637    -7.9392  186.556244"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "airports"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}