{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Spatial Operations and `rio`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Today we are going to use the spatial background we learned yesterday and apply it to dataset exploration of an AVIRIS image. We are going to use the `rio` module of `xarray`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Libraries for image visualization\n", "from IPython.display import Image\n", "from IPython.core.display import HTML " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import xarray as xr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First let's open the datafile. We are going to use the argument `engine='rasterio'` because we are using an ENVI file type. This is also the recommended method for opening TIF files." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "aviris_ds = xr.open_dataset('./data/subset_f131205t01p00r10rdn_e_sc01_ort_img', engine='rasterio')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Extract just the band DataArray from the dataset\n", "aviris = aviris_ds['band_data']" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'band_data' (band: 224, y: 300, x: 100)>\n",
       "[6720000 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * band         (band) int64 1 2 3 4 5 6 7 8 ... 218 219 220 221 222 223 224\n",
       "    xc           (y, x) float64 5.777e+05 5.777e+05 ... 5.814e+05 5.814e+05\n",
       "    yc           (y, x) float64 4.151e+06 4.151e+06 ... 4.148e+06 4.148e+06\n",
       "    spatial_ref  int64 0\n",
       "Dimensions without coordinates: y, x\n",
       "Attributes:\n",
       "    long_name:                 ('Band 1', 'Band 2', 'Band 3', 'Band 4', 'Band...\n",
       "    bands:                     224\n",
       "    band_names:                Band 1,Band 2,Band 3,Band 4,Band 5,Band 6,Band...\n",
       "    byte_order:                0\n",
       "    coordinate_system_string:  PROJCS["unnamed",GEOGCS["GCS_WGS_1984",DATUM["...\n",
       "    data_type:                 2\n",
       "    description:               ./output_data/subset_f131205t01p00r10rdn_e_sc0...\n",
       "    file_type:                 ENVI Standard\n",
       "    header_offset:             0\n",
       "    interleave:                bsq\n",
       "    lines:                     300\n",
       "    samples:                   100
" ], "text/plain": [ "\n", "[6720000 values with dtype=float32]\n", "Coordinates:\n", " * band (band) int64 1 2 3 4 5 6 7 8 ... 218 219 220 221 222 223 224\n", " xc (y, x) float64 5.777e+05 5.777e+05 ... 5.814e+05 5.814e+05\n", " yc (y, x) float64 4.151e+06 4.151e+06 ... 4.148e+06 4.148e+06\n", " spatial_ref int64 0\n", "Dimensions without coordinates: y, x\n", "Attributes:\n", " long_name: ('Band 1', 'Band 2', 'Band 3', 'Band 4', 'Band...\n", " bands: 224\n", " band_names: Band 1,Band 2,Band 3,Band 4,Band 5,Band 6,Band...\n", " byte_order: 0\n", " coordinate_system_string: PROJCS[\"unnamed\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"...\n", " data_type: 2\n", " description: ./output_data/subset_f131205t01p00r10rdn_e_sc0...\n", " file_type: ENVI Standard\n", " header_offset: 0\n", " interleave: bsq\n", " lines: 300\n", " samples: 100" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "aviris" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we are opening TIF or ENVI files using the `rasterio` engine is that we get an additional set of attributes and methods that give us information about the image. We access this using `.rio`" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(577666.150342955, 4155829.70969225, 579276.150342955, 4150999.70969225)\n", "EPSG:32610\n", "300\n", "100\n", "None\n", "| 14.08, 7.81, 577666.15|\n", "| 7.81,-14.08, 4150999.71|\n", "| 0.00, 0.00, 1.00|\n" ] } ], "source": [ "print(aviris.rio.bounds()) # bounding box\n", "print(aviris.rio.crs) # coordinate reference system\n", "print(aviris.rio.height) # number of pixels tall\n", "print(aviris.rio.width) # number of pixels wide\n", "print(aviris.rio.nodata) # a nodata fill value, if one has been set\n", "print(aviris.rio.transform()) # A CRS property called an affine transformation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These metadata can be really helpful for orienting you to your data.\n", "\n", "One field I want to make a note about is the nodata field. The nodata value is sometimes not set, but that doesn't mean that there isn't a nodata value in the dataset. There still may be a nodata value present in the dataset even if one is not set." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{admonition} 📝 Check your understanding\n", ":class: tip\n", "\n", "Look at the EPSG code of this dataset. Reference the figure below from the spatial data lecture to figure out approximately where in the world the dataset is located. Is the image in the northern or southern hemisphere?\n", "\n", ":::" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Image(url= \"http://www.dmap.co.uk/utmworld.gif\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although the EPSG code is useful for general orientation and an overall sanity check, we usually want to know more specifically where our data is located. We can do this by extracting the bounds and converting the coordinates to EPSG 4326.\n", "\n", "The general syntax for that is:\n", "\n", "```\n", "t = Transformer.from_crs(INPUT_PROJ, OUTPUT_PROJ, always_xy=True).transform\n", "transform(t, SHAPELY_POINT)\n", "```\n", "Where everything in ALL_CAPS is a variable you will be inserting yourself.\n", "\n", "Here's an example" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from pyproj import Transformer\n", "from shapely.ops import transform" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# Get the bounding box\n", "bbox = aviris.rio.bounds()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "from shapely.geometry import box" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# Create a shapely object with the bounding box\n", "bbox_geom = box(bbox[0], bbox[3], bbox[2], bbox[1])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# Consider this as a copy-and-paste able code snippet that you can use if you ever need to\n", "# transform coordinates\n", "t = Transformer.from_crs('epsg:32610', 'epsg:4326', always_xy=True).transform\n", "bbox_4326 = transform(t, bbox_geom)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'POLYGON ((-122.10308958299767 37.50252918542879, -122.10256826050976 37.5460596647619, -122.12079199682005 37.54619677057729, -122.12130273640575 37.50266607712935, -122.10308958299767 37.50252918542879))'" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# View the output as a wkt\n", "bbox_4326.wkt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3.10.4 ('sarp')", "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.10.4" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "c6a7549e0a21eb75b694ff081490f6e8db8aea94086d29d6a7e085c556b251de" } } }, "nbformat": 4, "nbformat_minor": 2 }