{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Using GeoPandas with Rasterio to sample point data\n", "\n", "This example shows how to use GeoPandas with Rasterio. [Rasterio](https://rasterio.readthedocs.io/en/latest/index.html) is a package for reading and writing raster data.\n", "\n", "In this example a set of vector points is used to sample raster data at those points.\n", "\n", "The raster data used is Copernicus Sentinel data 2018 for Sentinel data.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import geopandas\n", "import rasterio\n", "import matplotlib.pyplot as plt\n", "from shapely.geometry import Point" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create example vector data\n", "=============================\n", "\n", "Generate a geodataframe from a set of points\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create sampling points\n", "points = [\n", " Point(625466, 5621289),\n", " Point(626082, 5621627),\n", " Point(627116, 5621680),\n", " Point(625095, 5622358),\n", "]\n", "gdf = geopandas.GeoDataFrame([1, 2, 3, 4], geometry=points, crs=32630)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ``GeoDataFrame`` looks like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gdf.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Open the raster data\n", "=============================\n", "\n", "Use ``rasterio`` to open the raster data to be sampled" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "src = rasterio.open(\"s2a_l2a_fishbourne.tif\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see the raster data with the point data overlaid.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "nbsphinx-thumbnail" ] }, "outputs": [], "source": [ "from rasterio.plot import show\n", "\n", "fig, ax = plt.subplots()\n", "\n", "# transform rasterio plot to real world coords\n", "extent = [src.bounds[0], src.bounds[2], src.bounds[1], src.bounds[3]]\n", "ax = rasterio.plot.show(src, extent=extent, ax=ax, cmap=\"pink\")\n", "\n", "gdf.plot(ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sampling the data\n", "===============\n", "Rasterio requires a list of the coordinates in x,y format rather than as the points that are in the geomentry column.\n", "\n", "This can be achieved using the code below" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "coord_list = [(x, y) for x, y in zip(gdf[\"geometry\"].x, gdf[\"geometry\"].y)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Carry out the sampling of the data and store the results in a new column called `value`. Note that if the image has more than one band, a value is returned for each band." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gdf[\"value\"] = [x for x in src.sample(coord_list)]\n", "gdf.head()" ] } ], "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.2" } }, "nbformat": 4, "nbformat_minor": 4 }