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)) 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) & (lon <= GeoT + xsize * GeoT), 'longitude is not on the map' assert (lat <= GeoT) & (lat >= GeoT + ysize * GeoT), 'latitude is not on the map' location = GeoT xpixel = -1 while location <= lon: location += GeoT xpixel += 1 location = GeoT ypixel = -1 while location >= lat: location += GeoT ypixel += 1 return xpixel, ypixel