Comparing CHIRPS with observed data

Compare satellite derived data with field observations

In previous exercises we have used Python and CHIRPS rainfall data to calculate several statistics about rainfall in Vietnam. The python scripts we used to do this always followed to following steps:

1. Define which data we want to use/open
2. Create a new variable to store results
3. Open the data from step 1
4. Do a calculation
5. Save the result in the variable from step 2

Using these 5 steps, we are now going to compare observed rainfall data from a meteo-station in Hanoi to data from CHIRPS. The observed rainfall data can be downloaded here (or get it from the USB stick).

In order to compare the observed values with CHIRPS data, we need to follow the five steps described above. To do this, you can use some of these functions:

Function 1 This function requires 1 input called “folder”. When you input a string of the location of your folder, it returns a list of all the tif-files in the folder.

def ListFilesInFolder(folder):
    list_of_files = [os.path.join(folder,fn) for fn in next(os.walk(folder))[2] if fn[-4:] == '.tif']
    return list_of_files

Function 2 This functions open a tif-file and returns the values of the file as an python array

def OpenAsArray(filename, Bandnumber = 1):
    DataSet = gdal.Open(filename, gdal.GA_ReadOnly)
    Band = DataSet.GetRasterBand(Bandnumber)
    Array = Band.ReadAsArray()
    return Array

Function 3 This function converts WGS84 coordinates to pixel coordinates.

def PixelCoordinates(lon,lat,filename):
    SourceDS = gdal.Open(filename, gdal.GA_ReadOnly)
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    assert (lon >= GeoT[0]) & (lon <= GeoT[0] + xsize * GeoT[1]), 'longitude is not on the map'
    assert (lat <= GeoT[3]) & (lat >= GeoT[3] + ysize * GeoT[5]), 'latitude is not on the map'
    location = GeoT[0]
    xpixel = -1
    while location <= lon:         location += GeoT[1]         xpixel += 1     location = GeoT[3]     ypixel = -1     while location >= lat:
        location += GeoT[5]
        ypixel += 1
    return xpixel, ypixel