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/

Transcript

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

   View Slide

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

   View 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 Slide

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

   View 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 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 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 Slide

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

   View Slide

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

   View 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 Slide

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

    View Slide

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

    View 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 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 Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View 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 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 Slide

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

    View 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 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 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 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 Slide

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

    View Slide

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

    View 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 Slide

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

    View Slide

32. Inspect  and  download:  
    32  
    

    View 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 Slide

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

    View 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 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 Slide

37. PYTHON  MAPPING  IN  DETAIL  
    A  worked  example  
    37  
    

    View Slide

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

    View Slide

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

    View 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 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 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 Slide

43. View Slide

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  
    

    View Slide

45. View Slide

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  
    

    View Slide

47. View Slide

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  
    

    View Slide

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

    View Slide