Chris Waigl - Satellite mapping for everyone

D5710b3bca38f1233274b4cbc523dc4b?s=47 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/

D5710b3bca38f1233274b4cbc523dc4b?s=128

PyCon 2015

April 18, 2015
Tweet

Transcript

  1. Satellite  mapping  for  everyone   Chris  Waigl   @chrys  

    Talk  presented  at  PyCon  2015,  Montréal   1  
  2. SATELLITE  MAPPING  FOR  PYTHON   PROGRAMMERS   …  and  Python

     programming  for  satellite  mappers   2  
  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  
  4. …  but  they  have  all  focused  on  vector  data,  

    not  raster  data.   4   Prakash,  Waigl,  McNabb,  UAF  
  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  
  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  
  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  
  8. WHY  SATELLITE  MAPPING?   …  and  first:  what  is  satellite

     mapping?     8  
  9. Polar  orbiting  (=  sun-­‐synchronous)  satellite   (EUMETSAT  =  European  met

     satellite  agency)   9   Video  at  https://youtu.be/y_jM_BxQGvE    
  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  
  11. Example:  Maybe  you  live  in  a  place  like  this.  

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

    2012  vs.  2013,  May.   12   upenn.edu   earthobservatory.nasa.gov  
  13. 30  ×  30  m  imagery:  visualize  changes  on  a  city

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

     or   regional  scale.  Example:  Las  Vegas,  1984  vs  2009   14   earthobservatory.nasa.gov  
  15. Higher  resolution  images  are  available.   15   Alaska  SDMI,

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

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

    Alaska  SDMI,   SPOT5,  source:   GINA,  UAF  
  18. Visual  data  journalism  :  The  ProPublica  project   on  Louisiana’s

     disappearing  coastline   18  
  19. A  LOOK  AT  SOME  DATA   Where  to  get  it,

     what  it  looks  like   19  
  20. Back  to  raster  and  vector  data.  How  do  we  

    process  this?     20   Prakash,  Waigl,  McNabb,  UAF  
  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   …  
  22. Simple  case:  highly  pre-­‐processed  file  (NASA).   Imagine  you  live

     close  to  a  mine  like  this:     22   mtpolley_oli_2014210_geo.tif  
  23. Simple  case:  highly  pre-­‐processed  file  (NASA).   …  and  then

     one  day  this  happens:     23   mtpolley_oli_2014217_geo.tif  
  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  
  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  
  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  
  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  
  28. Accessing  the  data  with  GDAL   >>>  ds  =  gdal.Open(fpath)

      >>>  band1  =  ds.GetRasterBand(1)   >>>  data  =  band1.ReadAsArray()   >>>  print(type(data))   <type  'numpy.ndarray'>   28  
  29. Accessing  the  data  with  rasterio   >>>  with  rasterio.open(fname)  as

     src:          b1,  b2,  b3  =  src.read()   >>>  print(type(b1))   <type  'numpy.ndarray'>     29  
  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  
  31. Full-­‐scale  Landsat  data  recommended  portal   (and  many  others):  earthexplorer.usgs.gov

      31  
  32. Inspect  and  download:   32  

  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  
  34. Introduction  to  remote  sensing  in  1  min   34  

    band  =  channel  =  layer       radiometric  resolution  =    bit  depth  (8  bit,  12  bit,  16  bit)  
  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  
  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.  
  37. PYTHON  MAPPING  IN  DETAIL   A  worked  example   37

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

    Eastern  Canada   38  
  39. National  Snow  and  Ice  Data  Center  (NSIDC):  How   much

     ice  is  there  in  the  winter  in  St.  Laurence  Bay   39  
  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  
  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  
  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  
  43. 43  

  44. 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  
  45. 45  

  46. 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  
  47. 47  

  48. 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  
  49. COLLEGE OF NATURAL SCIENCE & MATHEMATICS University of Alaska Fairbanks

    NASA  Earth  and  Space   Science  Fellowship   Thank  you!     49