SRCCON 2014 Raster Data

srccon.md

Telling Stories with Raster Data

SRCCON 2014 | Al Shaw | @A_L | almshaw@gmail.com

This file: http://j.mp/srccon-raster

Announcement: SimpleTiles/SimplerTiles raster

• Shameless plug: SimpleTiles + SimplerTiles now support raster data.

• SimplerTiles slippy map example:

    map = SimplerTiles::Map.new do |m|
      m.slippy params[:x].to_i, params[:y].to_i, params[:z].to_i
      m.raster_layer "path/to/geo.tif"
    end
  

• SimplerTiles Raster DEMO <--- Image is scene LC80130322014100LGN00 available on the thumb drive, hastily pansharpened in Photoshop

Vector maps are everywhere, thanks to Google Maps and Census data

• NYT 2010 Census Explorer
• WaPo Super ZIPs
• Fusion tables, Google Maps API, Mapbox/OpenStreetMap, d3, Landline

(via http://flowingdata.com/2014/07/07/19-maps-that-will-blow-your-mind/)

Raster data is images associated with geographic data

• Crafting the Blue Marble
• NYPL MapWarper: Manhattan, 1891 Georectified
• Raster to Vector! NYPL Perris Atlas
• PublicLab: Grand Isle, Louisiana mapping
• NYT: Battle for the Baiji Oil Refinery
• NYT: A Rogue State Along Two Rivers
• MapBox: Monitoring Oil Reserves from Space
• Drone pilot finds “river of blood” outside Dallas meatpacking plant
• NASA Earth Observatory: Glacier Retreat near Indonesia's Grasberg Mine
• Mapbox: Cloudless Atlas
• Google Earth Engine + Timelapses
• Raster imagery after a disaster: Hurricane Sandy

When reporting ruins a good story: Georeferencing FEMA's paper flood maps after Weld County, Colo. flooding

Our Bands Could Be Your Life

• Three kinds of resolution: spatial, temporal, wavelength
• Comparisons:

(via: http://www.geog.ucsb.edu/~jeff/projects/thesis/)

• Landsat 7/8 bands compared

• Common Landsat 8 band combinations
• NDVI: Measuring vegetation. Quick 'n dirty NDVI Python script.
• Pansharpening in Photoshop, Orfeo

Where do I get it?

• NOAA Historical Charts
• USGS Earth Explorer
• Landsat via Google EarthEngine
• NASA LPDAAC
• Blue Marble (Derived from Modis)
• Emerging private companies: PlanetLab, Skybox, etc.
• Georeference any image you want, using QGIS!

How do I use it?

• Georeference in QGIS

• Check the geo headers of a file with gdalinfo

    Driver: GTiff/GeoTIFF
    Files: LC80130322014100LGN00_B1.TIF
    Size is 7661, 7791
    Coordinate System is:
    PROJCS["WGS 84 / UTM zone 18N",
        GEOGCS["WGS 84",
            DATUM["WGS_1984",
                SPHEROID["WGS 84",6378137,298.257223563,
                    AUTHORITY["EPSG","7030"]],
                AUTHORITY["EPSG","6326"]],
            PRIMEM["Greenwich",0],
            UNIT["degree",0.0174532925199433],
            AUTHORITY["EPSG","4326"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",-75],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",500000],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AUTHORITY["EPSG","32618"]]
    Origin = (524085.000000000000000,4582215.000000000000000)
    Pixel Size = (30.000000000000000,-30.000000000000000)
    Metadata:
      AREA_OR_POINT=Point
    Image Structure Metadata:
      INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  (  524085.000, 4582215.000) ( 74d42'42.88"W, 41d23'27.98"N)
    Lower Left  (  524085.000, 4348485.000) ( 74d43'14.63"W, 39d17' 7.42"N)
    Upper Right (  753915.000, 4582215.000) ( 71d57'53.48"W, 41d21' 5.12"N)
    Lower Right (  753915.000, 4348485.000) ( 72d 3'27.41"W, 39d14'54.75"N)
    Center      (  639000.000, 4465350.000) ( 73d21'49.68"W, 40d19'37.63"N)
    Band 1 Block=7661x1 Type=UInt16, ColorInterp=Gray
  

• Hint: find the EPSG SRS code from Proj string with prj2epsg

• Hint: reproject coordinates together with cs2cs, also available online

• Reproject images with gdalwarp

    for band in {4,3,2}
    do
      gdalwarp -t_srs EPSG:3857 img_$band.tif img_$band_projected.tif
    done
  

• Stitch scenes together with gdalwarp

    for scenes in $(ls)
    do
      gdalwarp --config GDAL_CACHEMAX 3000 -wm 3000 $scenes stitched.tif
    done
  

• Merge bands together with gdal_merge.py

    # preserve nodata
    gdal_merge.py -n 0 -a_nodata 0 -o out.tiff image1_band1.tiff image1_band2.tif
  

• Add geographic headers to any TIF with gdal_edit.py

    gdal_edit.py -a_srs EPSG:3857 scene.tif
  

• Manipulate and combine images using convert, then put geo headers back in with listgeo

     # generate a tfw file
     listgeo -tfw scene.tif
  

• Do calculations between bands with gdal_calc.py

• Putting it all together: charlie-loyd.rake based on Charlie Loyd's post.

• Mapbox guide to processing satellite imagery

• The old fashioned way: Rob Simmon's "How to Make a True-Color Landsat 8 Image"

• GDAL cheat sheet from Derek Watkins

• GDAL scripts from Dan Stahlke

What we're working on

• True color (4,3,2) Landsat 8 image of southeast Louisiana

• False color (7,5,3) image of southeast Louisiana to accentuate coastal erosion

• Land loss in Lafitte, La. since 1956

Big problems/discussion

• How to acquire and process data programmatically, to create "real time" raster apps?
• How do we stitch and blend together large amounts of data, especially with cloud cover?
• How can we store and serve massive GeoTIFFs?
• Design/presentation issues: How do we show false color images to readers without being misleading?
• Deciding between vector and raster maps
• Processing your own imagery vs. using Google or Mapbox imagery
• The future of raster - privacy concerns with increased spatial and temporal resolution
• The uncanny valley of sharpening and post-processing: when does an image become an illustration?

On the thumb drive

• Landsat 8 scenes LC80130322014100LGN00 and LC80140322014187LGN00 (NYC meets Philly <3 <3 <3)