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

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

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  

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

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  

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  

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  

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

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

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  

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

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

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

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

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

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

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

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

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

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

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   …  

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

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

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

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  

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  

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  

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

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

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  

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

Inspect  and  download:   32  

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  

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

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  

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.  

PYTHON  MAPPING  IN  DETAIL   A  worked  example   37  

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

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

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  

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  

…  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  

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  

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  

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  

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