notebooks/00_create_naip_on_aws_gpkg.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NAIP On AWS\n",
    "\n",
    "This Jupyter notebook provides code to scrape the [NAIP on AWS](https://registry.opendata.aws/naip/) manifest and create a GPKG file that provides a geospatial footprint for all available imagery.\n",
    "\n",
    "Then, if we are interested in a given OSM ID, it's easy to look up which NAIP images intersect this ID and download them directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1150/1150 [00:00<00:00, 125232.36it/s]\n"
     ]
    }
   ],
   "source": [
    "import pathlib\n",
    "\n",
    "from tqdm import tqdm\n",
    "from osgeo import ogr, osr\n",
    "ogr.UseExceptions()\n",
    "\n",
    "from rsc.common import aws_naip\n",
    "from rsc.common.aws_naip import AWS_PATH\n",
    "\n",
    "# Get NAIP manifest\n",
    "manifest = aws_naip.get_naip_manifest()\n",
    "\n",
    "# Filter out shapefiles to download\n",
    "shp = [e for e in manifest if e.split('.')[-1].lower() in \\\n",
    "    ('shp', 'dbf', 'shx', 'prj', 'sbn')]\n",
    "\n",
    "# Fetch all the shapefiles\n",
    "for object_name in tqdm(shp):\n",
    "    aws_naip.get_naip_file(object_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conv(s: str) -> str:\n",
    "    \"\"\" Convert object paths as seen in manifest to those that might be seen in the shapefiles.\n",
    "        It's silly they don't match.\"\"\"\n",
    "    return '%s.tif' % '_'.join(s.split('_')[:6])\n",
    "\n",
    "# Read the manifest, and convert the TIF files to those that might be seen in the shapefile metadata\n",
    "with open(AWS_PATH / 'manifest.txt', 'r') as f:\n",
    "    mani = {conv(pathlib.Path(p).stem): p for p in (e.strip() for e in f.readlines()) if p.endswith('.tif')}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do all the work!\n",
    "\n",
    "# Create SRS (EPSG:4326: WGS-84 decimal degrees)\n",
    "srs = osr.SpatialReference()\n",
    "srs.ImportFromEPSG(4326)\n",
    "srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)\n",
    "\n",
    "# Create GPKG file for writing\n",
    "driver: ogr.Driver = ogr.GetDriverByName('GPKG')\n",
    "ds_w: ogr.DataSource = driver.CreateDataSource(str(AWS_PATH / 'naip_on_aws.gpkg'))\n",
    "layer_w: ogr.Layer = ds_w.CreateLayer('footprints', srs=srs, geom_type=ogr.wkbPolygon)\n",
    "\n",
    "# Define output fields\n",
    "state_field = ogr.FieldDefn('STATE', ogr.OFTString)\n",
    "band_field = ogr.FieldDefn('BAND', ogr.OFTString)\n",
    "usgs_id_field = ogr.FieldDefn('USGSID', ogr.OFTString)\n",
    "src_img_date_field = ogr.FieldDefn('SRCIMGDATE', ogr.OFTString)\n",
    "filename_field = ogr.FieldDefn('FILENAME', ogr.OFTString)\n",
    "object_field = ogr.FieldDefn('OBJECT', ogr.OFTString)\n",
    "\n",
    "# Create output fields in layer\n",
    "layer_w.CreateField(state_field)\n",
    "layer_w.CreateField(band_field)\n",
    "layer_w.CreateField(usgs_id_field)\n",
    "layer_w.CreateField(src_img_date_field)\n",
    "layer_w.CreateField(filename_field)\n",
    "layer_w.CreateField(object_field)\n",
    "\n",
    "# Get layer feature definition to load in features\n",
    "feat_defn = layer_w.GetLayerDefn()\n",
    "\n",
    "# Loop through all fetched shapefiles\n",
    "for p in AWS_PATH.rglob('*.shp'):\n",
    "\n",
    "    # Load them in OGR, get layer and spatial reference\n",
    "    ds_r: ogr.DataSource = ogr.Open(str(p))\n",
    "    layer_r: ogr.Layer = ds_r.GetLayer()\n",
    "    srs_r: osr.SpatialReference = layer_r.GetSpatialRef()\n",
    "    srs_r.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)\n",
    "\n",
    "    # Loop throught the features in the layer\n",
    "    for _ in range(layer_r.GetFeatureCount()):\n",
    "        feat_r: ogr.Feature = layer_r.GetNextFeature()\n",
    "\n",
    "        # Quickly crosscheck with manifest. Skip if not in there\n",
    "        filename = feat_r.GetFieldAsString('FileName')\n",
    "        filename_conv = conv(filename.split('.')[0])\n",
    "        if not filename_conv in mani:\n",
    "            continue\n",
    "\n",
    "        # Parse remaining metadata\n",
    "        try:\n",
    "            state = feat_r.GetFieldAsString('ST')\n",
    "        except RuntimeError:\n",
    "            state = feat_r.GetFieldAsString('QUADST')\n",
    "        band = feat_r.GetFieldAsString('Band')\n",
    "        usgs_id = feat_r.GetFieldAsString('USGSID')\n",
    "        src_img_date = feat_r.GetFieldAsString('SrcImgDate')\n",
    "\n",
    "        # Fetch geometry and convert to desired spatial reference\n",
    "        trans = osr.CoordinateTransformation(srs_r, srs)\n",
    "        geom = ogr.CreateGeometryFromWkt(feat_r.GetGeometryRef().ExportToWkt())\n",
    "        geom.Transform(trans)\n",
    "\n",
    "        # Create our new feature\n",
    "        feat_w = ogr.Feature(feat_defn)\n",
    "        feat_w.SetGeometry(geom)\n",
    "        feat_w.SetField('STATE', state)\n",
    "        feat_w.SetField('BAND', band)\n",
    "        feat_w.SetField('USGSID', usgs_id)\n",
    "        feat_w.SetField('SRCIMGDATE', src_img_date)\n",
    "        feat_w.SetField('FILENAME', filename)\n",
    "        feat_w.SetField('OBJECT', mani[filename_conv])\n",
    "\n",
    "        # Save!\n",
    "        layer_w.CreateFeature(feat_w)\n",
    "\n",
    "        # Cleanup features\n",
    "        feat_w = None\n",
    "        feat_r = None\n",
    "\n",
    "    # Cleanup read dataset\n",
    "    layer_r = None\n",
    "    ds_r = None\n",
    "\n",
    "# Cleanup write dataset\n",
    "layer_w = None\n",
    "ds_w = None"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.8 64-bit",
   "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.8"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "4d102384ded633c24f1031e288c2ecf1ababc4ef37e402995ad37064232eefd1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}