Create a gist now

Instantly share code, notes, and snippets.

@ashaw /srccon.md Secret
Last active Aug 29, 2015

Embed
What would you like to do?
SRCCON 2014 Raster Data

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

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

Raster data is images associated with geographic data

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

Where do I get it?

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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment