In [1]:

    

import numpy as np
import matplotlib.pyplot as plt

from osgeo import gdal
from osgeo import gdal_array

import rasterio
import cartopy.crs as ccrs

## DEM plot via rasterio, pyplot imshow
##  Live 8.5  *  darkblue-b
##


    

In [2]:

    

## Rasterio
##   Clean and fast and geospatial raster I/O for Python programmers who use Numpy

with rasterio.open('SanMateo_CA.tif') as src:
    data = src.read()
    data = np.transpose( data, [1,2,0])


## show selected data attributes
## -------------------------------------
# type(data)  numpy.ma.core.MaskedArray
# data.ndim   3
# data.shape  (1080, 864, 1)
# data.dtype  dtype('float32')

## NOTE: when using rasterio and matplotlib only,
##  the column order for the trivial case of lat/lon/measure
##  is NOT handled automatically..  
##  column reordering is REQUIRED   np.transpose()


    

In [3]:

    

## Rasterio supplies a simple dictionary of important GeoTIFF metadata
src.meta


    

    
Out[3]:





{'affine': Affine(0.0002777777777780012, 0.0, -122.44000000003291,
       0.0, -0.0002777777777779992, 37.69999999999724),
 'count': 1,
 'crs': {'init': u'epsg:4269'},
 'driver': u'GTiff',
 'dtype': numpy.float32,
 'height': 1080,
 'nodata': -3.4028234663852886e+38,
 'transform': (-122.44000000003291,
  0.0002777777777780012,
  0.0,
  37.69999999999724,
  0.0,
  -0.0002777777777779992),
 'width': 864}

In [4]:

    

## numpy ndarray indexing example
##  show the measure values for a 2x2 area
##  no-data options are available using a masked array
##    http://docs.scipy.org/doc/numpy/reference/maskedarray.generic.html#what-is-a-masked-array

data[0:2,0:2,0]


    

    
Out[4]:





masked_array(data =
 [[ 203.41319275  203.94624329]
 [ 198.21229553  197.35855103]],
             mask =
 False,
       fill_value = -3.40282e+38)

In [6]:

    

## pyplot Image-Show 
##  http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.imshow
##
##  Example of using matplotlib to directly display a GeoTIFF
##
##  Cartopy supplies a library of mapping transformations
##   use an idiom to calculate the correct display bounds
##
##  data[:,:,0] refers to a numpy ndarray: all-x, all-y, 0th measure
##

xmin = src.transform[0]
xmax = src.transform[0] + src.transform[1]*src.width
ymin = src.transform[3] + src.transform[5]*src.height
ymax = src.transform[3]

##   Mercator etc..
crs = ccrs.PlateCarree()

ax = plt.axes(projection=crs)
plt.imshow(  data[:,:,0], origin='upper', 
             extent=[xmin, xmax, ymin, ymax], 
             cmap=plt.get_cmap('gist_earth'),
             transform=crs, interpolation='nearest')


    

    
Out[6]:





<matplotlib.image.AxesImage at 0xacd99d0c>






    

In [7]:

    

## Show the same content, but with a colorbar legend for the data
##
plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())

#elev, crs, extent = srtm_composite(12, 47, 1, 1)
plt.imshow( data[:,:,0], transform=ccrs.PlateCarree(),
           cmap='gist_earth')
cb = plt.colorbar(orientation='vertical')
cb.set_label('Altitude')
plt.title("DEM")


    

    
Out[7]:





<matplotlib.text.Text at 0xac62568c>






    

In [8]:

    

%whos


    

    




Variable     Type                            Data/Info
------------------------------------------------------
ax           GeoAxesSubplot                  < GeoAxes: <cartopy.crs.P<...>e object at 0xad545f2c> >
cb           matplotlib.colorbar.Colorbar    <matplotlib.colorbar.Colo<...>r instance at 0xac6049ec>
ccrs         module                          <module 'cartopy.crs' fro<...>ackages/cartopy/crs.pyc'>
crs          PlateCarree                     <cartopy.crs.PlateCarree object at 0xadf643bc>
data         MaskedArray                     [[[  2.03413193e+02]\n  [<...>]\n  [  9.26997986e+01]]]
gdal         module                          <module 'osgeo.gdal' from<...>packages/osgeo/gdal.pyc'>
gdal_array   module                          <module 'osgeo.gdal_array<...>es/osgeo/gdal_array.pyc'>
np           module                          <module 'numpy' from '/us<...>ages/numpy/__init__.pyc'>
plt          module                          <module 'matplotlib.pyplo<...>7/matplotlib/pyplot.pyc'>
rasterio     module                          <module 'rasterio' from '<...>s/rasterio/__init__.pyc'>
src          RasterReader                    <closed RasterReader name<...>anMateo_CA.tif' mode='r'>
xmax         float                           -122.2
xmin         float                           -122.44
ymax         float                           37.7
ymin         float                           37.4

In [ ]: