{{{#!python
#example how to make a sub-sampling of an array
demsamp=10 #takes e.g. every etnth pixel
ELEV=ELEV[::demsamp,::demsamp]
}}}


== Code to read GTOPO30 data ==
{{{#!python
def read_gtopo30(filename):
    '''
    routine to read GTOPO 30 data
    
    INPUT:
    filename without extension
    
    OUTPUT:
    DEM and LAT/LON GRID
    '''
    
    print 'Reading GTOPO30 data: ', filename
    
    #/// read header ///
    hdrname = filename + '.HDR'
    data = [line.split() for line in open(hdrname).readlines()]
    hdr = dict(data)

    ulx,uly = float(hdr['ULXMAP']),float(hdr['ULYMAP'])
    nx, ny = int(hdr['NCOLS']), int(hdr['NROWS'])
    dx, dy = float(hdr['XDIM']), float(hdr['YDIM'])
    nodata = float(hdr['NODATA'])

    #/// read data ///
    elev = None
    file = open(filename + '.DEM','rb')
    a=array.array('h')
    a.read(file,nx*ny)
    a.byteswap()
    elev=reshape(a,(ny,nx))
    elev=flipud(elev)
    a=None
    
    #/// geometry ///
    dem_lon = arange(ulx,ulx+nx*dx,dx)
    dem_lat = arange(uly-ny*dy,uly,dy)
    LON,LAT=meshgrid(dem_lon,dem_lat)
    #elev[elev==nodata]=
    
    return LON,LAT,elev

}}}