{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example demonstrates how to write out a RasterFrame to a GeoTIFF. RasterFrames provides a specialized Spark DataFrame writer for rendering a RasterFrame to a GeoTIFF. It is an expensive operation since creating a GeoTIFF requires that all of the data be in the memory of one computer. \n",
    "\n",
    "In this example, we use a lower resolution collection and limit the area of interest to create a manageable data size. We run through the steps of 1) acquiring imagery scenes from the EarthAI Catalog, 2) using RasterFrames to read imagery, and 3) writing a RasterFrame to GeoTIFF.\n",
    "\n",
    "# Import Libraries\n",
    "\n",
    "We will start by importing all of the Python libraries used in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from earthai.init import *\n",
    "import earthai.chipping.strategy\n",
    "\n",
    "import pyspark.sql.functions as F\n",
    "import ipyleaflet\n",
    "import geopandas "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Query the EarthAI Catalog\n",
    "\n",
    "We read in a GeoJSON file containing U.S. state boundaries and filter the GeoDataFrame to North Carolina. \n",
    "\n",
    "We use the __geometry__ column in the GeoDataFrame to query the EarthAI catalog for MODIS surface reflectance data from September 1, 2020 covering North Carolina. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_url ='https://raw.githubusercontent.com/datasets/geo-admin1-us/master/data/admin1-us.geojson'\n",
    "states_gdf = geopandas.read_file(states_url)\n",
    "\n",
    "nc_gdf = states_gdf[states_gdf[\"name\"] == \"North Carolina\"]\n",
    "\n",
    "cat = earth_ondemand.read_catalog(nc_gdf.geometry, \n",
    "                                  start_datetime='2020-09-01', \n",
    "                                  end_datetime='2020-09-01', \n",
    "                                  collections='mcd43a4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in MODIS Imagery\n",
    "\n",
    "We join the catalog back to the GeoDataFrame containing North Carolina state boundaries in order to match the state boundary to the intersecting image scene. This step is critical for use of the chip reader since the chipping strategy needs the state boundary polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat = geopandas.sjoin(cat, nc_gdf, how='right').rename(columns={\"geometry\":\"nc_bounds\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since MODIS scenes are very large, we use `spark.read.chip` to read only the imagery intersecting North Carolina state boundaries. We read in the __B01__ (red), __B04__ (green), and __B03__ (blue) bands. The feature-aligned grid strategy creates a grid across North Carolina using the specified `tile_dimensions`, and returns the generated chips. \n",
    "\n",
    "_To view all of the available bands for the MODIS collection, you can run `earth_ondemand.item_assets('mcd43a4')`._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf = spark.read.chip(cat,\n",
    "                     catalog_col_names=['B01', 'B04', 'B03'],\n",
    "                     geometry_col_name='nc_bounds',\n",
    "                     chipping_strategy=earthai.chipping.strategy.FeatureAlignedGrid(256),\n",
    "                     tile_dimensions=(256,256))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To view the North Carolina state boundary and chip outlines on a map using [ipyleaflet](https://ipyleaflet.readthedocs.io/en/latest/), we need to reproject the chip geometries to the standard \"EPSG:4326\" projection before plotting the outlines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf = rf.select(F.col('B01').alias('red'), \n",
    "               F.col('B04').alias('green'),\n",
    "               F.col('B03').alias('blue')) \\\n",
    "       .withColumn('chip_geom_wgs84', st_reproject(rf_geometry('red'), rf_crs('red'), F.lit('EPSG:4326')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = ipyleaflet.Map(center=(nc_gdf.geometry.centroid.y.mean(), nc_gdf.geometry.centroid.x.mean()), zoom=6) \n",
    "\n",
    "nc_layer = ipyleaflet.GeoData(geo_dataframe=nc_gdf, \n",
    "                                      name='North Carolina', \n",
    "                                      style={'fillColor': '#f003fc', 'color': '#f003fc'}) \n",
    "\n",
    "chips_gdf = geopandas.GeoDataFrame(rf.select('chip_geom_wgs84').toPandas(),\n",
    "                                   geometry='chip_geom_wgs84', crs='EPSG:4326')\n",
    "\n",
    "chips_layer = ipyleaflet.GeoData(geo_dataframe=chips_gdf,\n",
    "                                  name='Chips', \n",
    "                                  style={'fillColor': '#32a852', 'color': '#32a852'}) \n",
    "\n",
    "m.add_layer(nc_layer) \n",
    "m.add_layer(chips_layer)\n",
    "m.add_control(ipyleaflet.LayersControl())\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write RasterFrame to GeoTIFF\n",
    "\n",
    "As mentioned earlier, it's important to limit the imagery size or downsample the resolution to create a single GeoTIFF of your imagery since all of the data to be encoded has to be in the memory of one computer. You can provide a `raster_dimensions` parameter, which will downsample (or upsample) your pixel resolution depending on the dimensions using bilinear resampling. If no `raster_dimensions` parameter is specified, the RasterFrame contents are written at full resolution. \n",
    "\n",
    "Since we used a lower resolution collection and limited our area of interest to North Carolina, we can write out the contents at full resolution. \n",
    "\n",
    "The RasterFrame writer will interpret a RasterFrame with three or four tile columns as red, green, blue, and optionally alpha bands. You can select any combination of three or four tile columns to create your intended color composite. Any other number of tile columns will result in a greyscale interpretation.\n",
    "\n",
    "We select the __red__, __green__, and __blue__ bands to write out a natural color image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf.select('red', 'green', 'blue').write.geotiff('geotiff-overview.tif', crs='EPSG:4326')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__geotiff-overview.tif__ will be written out to the same directory where your notebook resides. You can right click on the file in the left menu and select ___Download___ to save the file to your local machine for viewing in an external GIS software."
   ]
  }
 ],
 "metadata": {
  "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": 360051136412,
   "section_id": 360008732711,
   "title": "Writing a RasterFrame to a GeoTIFF File"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}