{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Geospatial analysts often have a set of vector data (points, lines, and polygons) that label areas where they want to look at satellite imagery.\n", "\n", "The code below works through some common steps for this kind of analysis:\n", "\n", "1. Reading the vector data from common formats\n", "1. Reading the Earth OnDemand catalog using the vector geometry\n", "1. Filtering and reasoning on the vector data and EO metadata\n", "1. Reading the imagery as a Spark DataFrame\n", "\n", "After this is complete, in our analysis we will compute the normalized difference built-up index (NDBI) for large US urban areas in 2018 in two analysis steps.\n", "\n", "Analysis steps:\n", "\n", "1. \"Burning-in\" or rasterizing the geometry alongside the image tiles\n", "1. Computing mean built-up index\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "from earthai.init import *\n", "\n", "import pyspark.sql.functions as F\n", "import geopandas\n", "import pandas as pd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "We will download some vector data from the US Census Bureau, delineating different urban areas.\n", "\n", "## Read Vector Data\n", "\n", "The data is in a shapefile. We will read it as a GeoDataFrame, which is a Pandas DataFrame with extras for handling vector data!\n", "\n", "You can also read shapefiles with Spark using `spark.read.shapefile`.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "! wget -nc -nv https://www2.census.gov/geo/tiger/TIGER2019/UAC/tl_2019_us_uac10.zip" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "! unzip -n tl_2019_us_uac10.zip" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "uac = geopandas.read_file('tl_2019_us_uac10.shp')\n", "print(\"Rows: \", len(uac))\n", "_ = uac.centroid.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Let's reduce the scope to the top 50 urban areas by the US Census reported land area.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "uac_50 = uac.nlargest(50, 'ALAND10')\n", "print(\"Rows: \", len(uac_50))\n", "_ = uac_50.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "uac_50.NAME10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Query the Earth OnDemand Catalog\n", "\n", "Just pass the GeoDataFrame into the `geo` parameter of the `earth_ondemand.read_catalog`. Let's look at which imagery sources we can use.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "earth_ondemand.collections()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Let's get some low cloud MODIS imagery for the first quarter of 2018. Note here the use of `simplify` to reduce the complexity of the city boundaries. The MODIS scenes are very, very large relative to the city boundaries so we still get all the relevant scenes but the query is much faster.\n", "\n", "Another interesting point here is we can pass in almost any kind of vector data. It can be a file-like python object such as an open vector file, the str path to a vector file, a `shapely` geometry object, etc. See `help(earth_ondemand.read_catalog)` for more.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "catalog = earth_ondemand.read_catalog(\n", " uac_50.geometry.simplify(0.2),\n", " '2018-01-01',\n", " '2018-03-31',\n", " max_cloud_cover=0,\n", " collections='mcd43a4'\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Here is the magic bit, we will do a @ref:[_spatial join_](https://astraeahelp.zendesk.com/hc/en-us/articles/360043896371) between our catalog and our city DataFrame.\n", "\n", "This attaches the city information to the scene information!\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "uac_catalog = geopandas.sjoin(uac_50, catalog.to_crs(uac_50.crs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Summarize the number of images for each city. The `id` column contains the unique scene identifier for the raster data product.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "uac_catalog.groupby('NAME10').count().id" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Filtering the Catalog\n", "\n", "Now that we know which city each image belongs to, let's find the lowest cloud image for each city.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "# have to reset the index because the left DF in the join has had its index repeated across rows\n", "uac_catalog.reset_index(inplace=True)\n", "least_cloud_img_per_city = uac_catalog.groupby('NAME10').eo_cloud_cover.idxmin()\n", "\n", "uac_least_cloudy = uac_catalog.loc[least_cloud_img_per_city]\n", "assert len(uac_least_cloudy) == len(uac_50)\n", "\n", "uac_least_cloudy[['index', 'NAME10', 'id', 'eo_cloud_cover', 'datetime']]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "uac_least_cloudy.columns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Read the Raster Data\n", "\n", "We will read the short-wave infrared (~1.6μm) and near-infrared from MODIS MCD43A4. These are bands __6__ and __2__ respectively. We can find this out from the `earth_ondemand.item_assets` function.\n", "\n", "Note the renaming we do. This allows us to forget about the specific band designations for the rest of the code. This is a good idea if we want to revisit a similar analysis with a different collection.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "earth_ondemand.item_assets('mcd43a4')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "rf = spark.read.raster(uac_least_cloudy, catalog_col_names=['B06', 'B02']) \\\n", " .withColumnRenamed('B06', 'swir16') \\\n", " .withColumnRenamed('B02', 'nir')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Create a new column reprojecting the city polygon to the same CRS as the imagery. Filter so we remove any rows where the image does not cover the city. We will also apply a very simple filter to remove `fill` values from the Landsat imagery. Note that our cloud cover is low, so we expect otherwise quite clean images.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "rf = rf.withColumn('geo_native',\n", " st_reproject('geometry', rf_mk_crs(str(uac_50.crs)), rf_crs('swir16'))) \\\n", " .filter(st_intersects(rf_geometry('swir16'), 'geo_native'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Take a look at a sample of the data.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "rf.select('swir16', 'nir', 'NAME10', st_centroid('geometry'), st_centroid(rf_geometry('swir16')))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Analysis\n", "\n", "### Burn-in the City Geometry\n", "\n", "We will create a new Tile aligned to the existing raster data that \"burns in\" the city polygon. This is discussed more under [zonal stats](https://rasterframes.io/zonal-algebra.html) on the RasterFrames docs. The cell below is going to work through the steps leaving the explanation to those docs.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "rf = rf.withColumn('city_tile',\n", " rf_rasterize('geo_native', rf_geometry('swir16'), F.lit(1), rf_dimensions('swir16').cols, rf_dimensions('swir16').rows)) \\\n", " .filter(rf_data_cells('city_tile') > 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Now we will mask the imagery with the rasterized city polygons.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "rf = rf.withColumn('swir16', rf_mask('swir16', 'city_tile')) \\\n", " .withColumn('nir', rf_mask('nir', 'city_tile'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Look at a sample\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": "peek", "option_string": "name = \"peek\", evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "rf.select('swir16', 'nir', 'city_tile', 'NAME10') \\\n", " .filter(rf_data_cells('city_tile') > 2000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "### Compute Mean Built-up Index\n", "\n", "More built-up / dense areas have higher values in this index.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "city_ndbi = rf.select('NAME10', rf_normalized_difference('swir16', 'nir').alias('ndbi')) \\\n", " .groupBy('NAME10') \\\n", " .agg(rf_agg_stats('ndbi').alias('ndbi_stats')) \\\n", " .select(F.col('NAME10').alias('city'),\n", " F.col('ndbi_stats').mean.alias('ndbi_mean'),\n", " F.col('ndbi_stats').data_cells.alias('ndbi_count'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "city_ndbi_pd = city_ndbi.toPandas()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "city_ndbi_pd.sort_values('ndbi_mean')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "auto", "collapsed": false, "jupyter": { "outputs_hidden": false }, "options": { "caption": false, "complete": true, "display_data": true, "display_stream": true, "dpi": 200, "echo": true, "evaluate": false, "f_env": null, "f_pos": "htpb", "f_size": [ 6, 4 ], "f_spines": true, "fig": true, "include": true, "name": null, "option_string": "evaluate=False", "results": "verbatim", "term": false, "wrap": "output" } }, "outputs": [], "source": [ "_ = city_ndbi_pd.plot.scatter('ndbi_count', 'ndbi_mean')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "This simple analysis could be improved by choosing more than just 1 scene per city. In some of the images above we see areas of missing data in the imagery that overlap with the city polygon.\n", "\n", "But we now have the basic steps for joining the vector data and the Earth OnDemand catalog, then working with the geometry and imagery together.\n" ] } ], "metadata": { "kernel_info": { "name": "echo" }, "kernelspec": { "display_name": "EarthAI Environment", "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.9" }, "zendesk": { "draft": true, "id": 360051411512, "section_id": 360008877912, "title": "Mean Built-up Index for USA Cities Using RasterFrames" } }, "nbformat": 4, "nbformat_minor": 4 }