Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Chris Waigl - Satellite mapping for everyone

PyCon 2015
April 18, 2015

Chris Waigl - Satellite mapping for everyone

Concerned about urban sprawl, landscape change or ecosystem recovery? Wildfire, drought or flooding? A vast amount of satellite data, collected since the 1970s, is freely available for your next mapping project. I will demonstrate how Python helps to make sense of odd scientific data and metadata formats and produce beautiful visualization and map products.

https://us.pycon.org/2015/schedule/presentation/425/

PyCon 2015

April 18, 2015
Tweet

More Decks by PyCon 2015

Other Decks in Technology

Transcript

  1. Satellite  mapping  for  everyone  
    Chris  Waigl  
    @chrys  
    Talk  presented  at  PyCon  2015,  Montréal  
    1  

    View full-size slide

  2. SATELLITE  MAPPING  FOR  PYTHON  
    PROGRAMMERS  
    …  and  Python  programming  for  satellite  mappers  
    2  

    View full-size slide

  3. There  have  been  many  great  talks  about  
    geospatial  Python  or  GIS  (geographic  
    information  systems)  at  previous  PyCons…  
    •  Mele  Sax-­‐Barnett  (PyCon  2014)  
    •  Jason  Scheirer  (PyCon  2012)  
    •  Zain  Memon  (PyCon  2012)  
    •  Sara  Safavi  (PyTexas  2014)    
    •  …  
     
    (They  are  all  on  pyvideo.org.)  
    3  

    View full-size slide

  4. …  but  they  have  all  focused  on  vector  data,  
    not  raster  data.  
    4  
    Prakash,  Waigl,  McNabb,  UAF  

    View full-size slide

  5. Both  are…  just  data.  We  can  deal  with  data.  
     vector  data  
     
    5  
    tabular/key-­‐value  data  
    geolocation  data  
    metadata  
    2D  arrays  (stacked)  
    geolocation  data  
    metadata  
    raster  data  

    View full-size slide

  6. Using  raster  data  is  not  (really)  harder.  
    •  OK,  fewer  people  do  it.    
    •  BUT  a  lot  of  data  is  freely  available.    
    •  Landsat,  for  example:  
    –  30  m  spatial  resolution  (mostly)  
    –  since  the  80s,  and  early  version  70s  
    –  any  point  on  the  globe  (except  the  poles)  
    –  public  domain  
    •  more  users  ⇒  better  tools  
    6  

    View full-size slide

  7. Satellite  mapping  for  Python  programmers  
    Why  would  you  want  to  do  this?  
    Where  would  you  get  data  and  what  does  it  
    look  like?  
    How  can  you  make  maps  from  it  most  
    efficiently?    
     
     
    7  

    View full-size slide

  8. WHY  SATELLITE  MAPPING?  
    …  and  first:  what  is  satellite  mapping?    
    8  

    View full-size slide

  9. Polar  orbiting  (=  sun-­‐synchronous)  satellite  
    (EUMETSAT  =  European  met  satellite  agency)  
    9  
    Video  at  https://youtu.be/y_jM_BxQGvE    

    View full-size slide

  10. Let’s  get  used  to  some  basic  concepts.  
    •  swath  width  (in  km)    
    •  orbital  inclination  
    •  repeat  interval  
    •  trade-­‐off  between  spatial  and  temporal  
    resolution  
    •  coordinate  reference  system  (and  map  
    projection)  
    10  

    View full-size slide

  11. Example:  Maybe  you  live  in  a  place  like  this.  
    2012  vs.  2013,  May.  
    11  
    upenn.edu  
    earthobservatory.nasa.gov  

    View full-size slide

  12. Example:  Maybe  you  live  in  a  place  like  this.  
    2012  vs.  2013,  May.  
    12  
    upenn.edu  
    earthobservatory.nasa.gov  

    View full-size slide

  13. 30  ×  30  m  imagery:  visualize  changes  on  a  city  or  
    regional  scale.  Example:  Las  Vegas,  1984  vs  2009  
    13  
    earthobservatory.nasa.gov  

    View full-size slide

  14. 30  ×  30  m  imagery:  visualize  changes  on  a  city  or  
    regional  scale.  Example:  Las  Vegas,  1984  vs  2009  
    14  
    earthobservatory.nasa.gov  

    View full-size slide

  15. Higher  resolution  images  are  available.  
    15  
    Alaska  SDMI,  
    SPOT5,  source:  
    GINA,  UAF  

    View full-size slide

  16. Higher  resolution  images  are  available.  
    16  
    Alaska  SDMI,  
    SPOT5,  source:  
    GINA,  UAF  

    View full-size slide

  17. SPOT5  2.5  m  resolution,  village  of  Togiak  
    17  
    Alaska  SDMI,  
    SPOT5,  source:  
    GINA,  UAF  

    View full-size slide

  18. Visual  data  journalism  :  The  ProPublica  project  
    on  Louisiana’s  disappearing  coastline  
    18  

    View full-size slide

  19. A  LOOK  AT  SOME  DATA  
    Where  to  get  it,  what  it  looks  like  
    19  

    View full-size slide

  20. Back  to  raster  and  vector  data.  How  do  we  
    process  this?    
    20  
    Prakash,  Waigl,  McNabb,  UAF  

    View full-size slide

  21. matplotlib  
    We  have  a  hierarchy  of  libraries  at  our  
    disposal  to  deal  with  the  data.      
    UPPER  CASE:  C  libraries  
    lower  case:  Python  libraries  
    21  
    HDF5  
    (h5py)  
    pygaarst  
    HDF4  
    (python-­‐hdf4)  
    GDAL/OGR  (osgeo)  
    shapely   fiona  
    GEOS  
    PROJ4  
    (pyproj)  
    numpy  
    rasterio   …  

    View full-size slide

  22. Simple  case:  highly  pre-­‐processed  file  (NASA).  
    Imagine  you  live  close  to  a  mine  like  this:    
    22  
    mtpolley_oli_2014210_geo.tif  

    View full-size slide

  23. Simple  case:  highly  pre-­‐processed  file  (NASA).  
    …  and  then  one  day  this  happens:    
    23  
    mtpolley_oli_2014217_geo.tif  

    View full-size slide

  24. Exploring  a  dataset  with  GDAL  
    >>>  from  osgeo  import  gdal  
    >>>  gdal.UseExceptions()  
    >>>  
    >>>  ds  =  gdal.Open(fpath)  
    >>>  print(ds.GetDescription())  
    mtpolley_oli_2014210_clip.tif  
    24  

    View full-size slide

  25. GDAL:  Dataset  attributes  and  methods  
    >>>  print(ds.RasterCount)  
    3  
    >>>  print(ds.RasterXSize)  
    1200  
    >>>  print(ds.RasterYSize)  
    900  
    >>>  print(ds.GetProjection())  
    PROJCS["WGS  84  /  UTM  zone  
    10N",GEOGCS["WGS  
    84",DATUM["WGS_1984",SPHEROID["WGS  84",
    6378137,298.257223563,  …]   25  

    View full-size slide

  26. Exploring  a  dataset  with  rasterio  
    >>>  import  rasterio  
    >>>  
    >>>  with  rasterio.open(fname)  as  src:  
    >>>          print(src.width,  src.height)  
    >>>          print(src.crs)  
    >>>          print(src.crs_wkt)  
    >>>          print(src.bounds)  
    >>>          print(src.count)  
    >>>          print(src.indexes)  
    26  

    View full-size slide

  27. Exploring  a  dataset  with  rasterio  
    1200  900  
    {'init':  u'epsg:32610'}  
    PROJCS["WGS  84  /  UTM  zone  10N"....]  
    BoundingBox(left=586762.5,  
    bottom=5816392.5,  right=604762.5,  
    top=5829892.5)  
    3  
    [1,  2,  3]  
    27  

    View full-size slide

  28. Accessing  the  data  with  GDAL  
    >>>  ds  =  gdal.Open(fpath)  
    >>>  band1  =  ds.GetRasterBand(1)  
    >>>  data  =  band1.ReadAsArray()  
    >>>  print(type(data))  
     
    28  

    View full-size slide

  29. Accessing  the  data  with  rasterio  
    >>>  with  rasterio.open(fname)  as  src:  
           b1,  b2,  b3  =  src.read()  
    >>>  print(type(b1))  
     
     
    29  

    View full-size slide

  30. Both  GDAL  and  rasterio  have  command-­‐line  
    interfaces  
    $  gdalinfo  mtpolley_oli_2014210_geo.tif  
    [a  whole  slew  of  information]  
    $  rio  info  mtpolley_oli_2014210_geo.tif  
    [a  concise  summary  as  in  JSON  format]    
     
    BUT:  GDAL  scripts  are  very  powerful,  especially  
    gdal_transform  and  gdal_merge.py:    
    •  cropping  (“subsetting”)  and  RGB  composites  
    •  read  the  source  code  of  the  Python  scripts  for  
    more  complicated  tasks  
    30  

    View full-size slide

  31. Full-­‐scale  Landsat  data  recommended  portal  
    (and  many  others):  earthexplorer.usgs.gov  
    31  

    View full-size slide

  32. Inspect  and  download:  
    32  

    View full-size slide

  33. Landsat  full  scene  data  
    Downloaded  file,  for  example:  
    LC80680152014261LGN00.tar.gz  (980  MB)  
     
     
     
     
    33  
    LC80680152013258LGN00_B1.TIF    (each  138  MB)  
    LC80680152013258LGN00_B10.TIF  
    LC80680152013258LGN00_B11.TIF  
    …  
    LC80680152013258LGN00_B8.TIF  
    LC80680152013258LGN00_B9.TIF  
    LC80680152013258LGN00_BQA.TIF  
    LC80680152013258LGN00_MTL.txt  

    View full-size slide

  34. Introduction  to  remote  sensing  in  1  min  
    34  
    band  =  channel  =  layer      
    radiometric  resolution  =    bit  depth  (8  bit,  12  bit,  16  bit)  

    View full-size slide

  35. Landsat  full  scene  data  basics:  
    Landsat  5,  7:  
    •  8  bit  ,  approx.  200-­‐500  MB  per  zip  file  
    •  8  bands:  use  7/5/3  for  false-­‐natural  color  RGB  
    Landsat  8:  
    •  16  bit,  approx.  1  GB  per  zip  file  
    •  11  bands:  use  7/6/4  for  false-­‐natural  color  RGB  
    35  

    View full-size slide

  36. Life  beyond  Landsat:  Datasets  may  have  
    subdatasets  and  more  complicated  geolocation  
    36  
    You  may  encounter  data  in  
    HDF4  or  HDF5  format.  
     
    This  data  has  subdatasets.  It  
    may  have  groups  of  
    subdatasets.  
     
    >>>  sd  =  ds.GetSubDatasets()    
    >>>  gdal.Open(sd[2][0])  
     
    Use  h5py  or  python-­‐hdf4  
    libraries  instead  of  GDAL.  

    View full-size slide

  37. PYTHON  MAPPING  IN  DETAIL  
    A  worked  example  
    37  

    View full-size slide

  38. Example:  Let’s  look  at  the  coast  right  here  in  
    Eastern  Canada  
    38  

    View full-size slide

  39. National  Snow  and  Ice  Data  Center  (NSIDC):  How  
    much  ice  is  there  in  the  winter  in  St.  Laurence  Bay  
    39  

    View full-size slide

  40. Step  1:  Read  the  documentation  
    •  GeoTiff  files  (one  of  the  options)  
    •  25x25  km  pixel  size  
    •  polar  stereographic  projection  
    •  8  bit  data  
    •  0-­‐250  maps  to  0-­‐100%  sea  ice  cover  
    •  251-­‐255:  used  for  special  data  values:  
    coastline,  land,  no  data,  bad  data…  
    40  

    View full-size slide

  41. Use  pygaarst  to  manipulte  the  data    
    (https://github.com/chryss/pygaarst)  
    >>>  import  numpy  as  np  
    >>>  from  pygaarst  import  raster  
    >>>  nsdic  =  raster.GeoTIFF(fn)  
    >>>  print(nsdic.nrow,  nsdic.ncol)  
    448  304  
    >>>  icedat  =  np.empty_like(nsdic.data)  
    >>>  icedat[:]  =  nsdic.data/2.5  
    >>>  maskdat  =  np.empty_like(nsdic.data)  
    >>>  maskdat  =  np.ma.MaskedArray(nsdic.data,  
    mask=nsdic.data  <=  250)  
    41  

    View full-size slide

  42. …  and  matplotlib  …  
    >>>  f,  ax  =  plt.subplots()  
    >>>  plt1  =  ax.imshow(icedat,  
    cmap="Blues_r",  vmin=0,  vmax=100,  
    interpolation='nearest')  
    >>>  plt2  =  ax.imshow(maskdat,  vmin=251,  
    vmax=255,  cmap=customcmap,  
    interpolation='nearest')  
    42  

    View full-size slide

  43. Find  the  Gulf  of  St.  Laurence  and  Anticosti  
    Island  …    
    >>>  lat_center,  lon_center  =  49.5,  -­‐63  
    >>>  geog2stereo  =  nsdic.coordtrans  
    >>>  x_center,  y_center  =  
    geog2stereo(lon_center,  lat_center)    
    >>>  i_center,  j_center  =  
    nsdic.xy2ij(x_center,  y_center,  
    precise=True)  
    >>>  print(i_center,  j_center)  
    407.930553598  97.4816627248  
     
    44  

    View full-size slide

  44. Use  the  basemap  toolkit  to  map  the  data  in  
    Albers  Equal  Area  projection  
    >>>  from  mpl_toolkits.basemap  import  Basemap  
     
    >>>  fig  =  plt.figure())  
    >>>  bmap  =  Basemap(width=…,  height=…,  
    resolution='i',  projection='aea',  lat_1=40.,  
    lat_2=60.,  lat_0=50.,  lon_0=-­‐63.)  
    >>>  x,  y  =  bmap(nsdic.Lon,  nsdic.Lat)  
    >>>  implt  =  bmap.pcolormesh(x,  y,  
    np.flipud(icedat),  …)  
     
    >>>  x,  y  =  bmap(lon_center,  lat_center)  
    >>>  bmap.scatter(x,  y,  marker='o',  c='tomato',  
    s=150,  zorder=20,  label='Anticosti',  alpha=0.7)  
    >>>  bmap.drawcoastlines(zorder=10)  
    >>>  bmap.drawrivers(color=water)  
    >>>  …    
    46  

    View full-size slide

  45. To  sum  up  
    •  Satellite  data  is  available  for  your  use  
    •  Tools  are  powerful  and  usability  is  
    improving  
    •  Find  the  last  example  as  an  IPython  
    notebook  on  
    http://nbviewer.ipython.org/gist/chryss/
    0e7894d9dc78a0c31708    
    •  If  you  would  like  to  get  together,  find  me  
    and  let’s  talk  
    48  

    View full-size slide

  46. COLLEGE OF NATURAL
    SCIENCE & MATHEMATICS
    University of Alaska Fairbanks
    NASA  Earth  and  Space  
    Science  Fellowship  
    Thank  you!    
    49  

    View full-size slide