1 #example how to make a sub-sampling of an array
   2 demsamp=10 #takes e.g. every etnth pixel
   3 ELEV=ELEV[::demsamp,::demsamp]

   1 def read_gtopo30(filename):
   2     '''
   3     routine to read GTOPO 30 data
   4     
   5     INPUT:
   6     filename without extension
   7     
   8     OUTPUT:
   9     DEM and LAT/LON GRID
  10     '''
  11     
  12     print 'Reading GTOPO30 data: ', filename
  13     
  14     #/// read header ///
  15     hdrname = filename + '.HDR'
  16     data = [line.split() for line in open(hdrname).readlines()]
  17     hdr = dict(data)
  18 
  19     ulx,uly = float(hdr['ULXMAP']),float(hdr['ULYMAP'])
  20     nx, ny = int(hdr['NCOLS']), int(hdr['NROWS'])
  21     dx, dy = float(hdr['XDIM']), float(hdr['YDIM'])
  22     nodata = float(hdr['NODATA'])
  23 
  24     #/// read data ///
  25     elev = None
  26     file = open(filename + '.DEM','rb')
  27     a=array.array('h')
  28     a.read(file,nx*ny)
  29     a.byteswap()
  30     elev=reshape(a,(ny,nx))
  31     elev=flipud(elev)
  32     a=None
  33     
  34     #/// geometry ///
  35     dem_lon = arange(ulx,ulx+nx*dx,dx)
  36     dem_lat = arange(uly-ny*dy,uly,dy)
  37     LON,LAT=meshgrid(dem_lon,dem_lat)
  38     #elev[elev==nodata]=
  39     
  40     return LON,LAT,elev