Converting GPCC gridded rainfall 1901-2010@0.5° to monthly Geotiffs

GPCC rainfall data contains the largest database worldwide, with over 90 000 weather-stations interpolated to a 0.5° grid. It includes S. Nicholsons dataset for Africa and is thus a good source for gridded monthly precipitation. It covers the period 1901-2010, however, it comes in a strange dataformat and it is a long way until we get such diagrams for the desired location:

GPCC annual rainfall averaged over my study area in northern Senegal

GPCC annual rainfall averaged over my study area in northern Senegal

As I couldn’t handle with the format provided at the DWD download site, I used the climate explorer to generate a ncdf-file. Use this link to get the global dataset: GPCC-V6.

Then this file is loaded into R. After loading the packages, the working directory which contains the file is set:

library(raster)
library(rgdal)
library(ncdf)
setwd("~/Dissertation/GPCC/")

Then the variable “prcp” is loaded into R and my study area clipped. To find out the correct variable, use gdalinfo e.g.

gpcc = brick("gpcc_V6_05.nc", varname="prcp")
gpcc
plot(gpcc)
e <- extent(-20, 5, 9, 26)
gpccclip <- crop(gpcc, e)

Then it’s exported as a multilayer Geotiff:

writeRaster(gpccclip, filename="gpcc_v6_westafr.tif", format="GTiff")

now I extract the individual rasters and create 1320 single files. I exit R for this:

cd ~/Dissertation/GPCC/
gdalinfo gpcc_v6_westafr.tif
for ((i = 1; i<= 1320; i++)) do gdal_translate -b $i -of GTiff gpcc_v6_westafr.tif `basename $i`.tif; done

Finaly the files are renamed. This may not be an elegant way, but it works:

for ((i=1;$i<=1320;i++)) do echo "mv $i.tif gpcc_v6_$(( ($i-1) / 12 + 1901 ))_$(( ($i-1) % 12 + 1 )).tif"; done >Umbenennungsscript.sh
sed -i 's/_1.tif/_01.tif/g' Umbenennungsscript.sh
sed -i 's/_2.tif/_02.tif/g' Umbenennungsscript.sh
sed -i 's/_3.tif/_03.tif/g' Umbenennungsscript.sh
sed -i 's/_4.tif/_04.tif/g' Umbenennungsscript.sh
sed -i 's/_5.tif/_05.tif/g' Umbenennungsscript.sh
sed -i 's/_6.tif/_06.tif/g' Umbenennungsscript.sh
sed -i 's/_7.tif/_07.tif/g' Umbenennungsscript.sh
sed -i 's/_8.tif/_08.tif/g' Umbenennungsscript.sh
sed -i 's/_9.tif/_09.tif/g' Umbenennungsscript.sh
chmod 0755 Umbenennungsscript.sh
./Umbenennungsscript.sh

The final result are 1320 Geotiffs, named gpcc_v6_1901_01.tif until gpcc_v6_2010_12.tif with monthly rainfall for West Africa, ready to be used in any GIS. Now we can easily create beatiful maps like this:

Mean annual rainfall 1950-2010 (source: GPCC v6)

Mean annual rainfall 1950-2010 (source: GPCC v6)

 

Advertisements

Automatically downloading and processing TRMM rainfall data

TRMM rainfall data maybe is the most accurate rainfall data derived from satellite measurements and a valuable source in regions with scarse weather-stations. It has a good spatial (0.25°) and temporal (daily) resolution and is available since 1998. However, downloading and processing may be a lot of work, if not scripted. The following script may be badly coded, but it works. All you need is the open source software R. This script is written for a Linux environment.

First, load some packages:

library(raster)
 library(rgdal)
 library(RCurl)

Then, download the data from the FTP server: ftp://disc2.nascom.nasa.gov/data/TRMM/Gridded/. This example downloads 3B43 monthly data version 7 for 2012. The first line has to be changed for the desired data:

url = "ftp://disc2.nascom.nasa.gov/data/TRMM/Gridded/3B43_V7/2012/"
filenames = getURL(url, ftp.use.epsv = FALSE, ftplistonly=TRUE, crlf=TRUE)
 filenames = paste(url, strsplit(filenames, "\r*\n") [[1]], sep="")
filenames = filenames[grep(filenames, pattern="*precipitation.accum")] # this line allows to download only files with a certain pattern, e.g. only certain months or days. "*precipitation.accum" means monthly accumulated rainfall here.
filenames # list all files
mapply(download.file, filenames, basename(filenames)) # download files

Now I create a virtual file as the downloaded TRMM data come as binarys. The VRT-file contains all the downloaded binary files with the appropriate geo-informations. To automate the process, the following script generates the VRT-file (TRMM.vrt) for all 2012 data. Change “3B43.12” accordingly to the downloaded data. Save this example as trmm2012.sh.

echo '<VRTDataset rasterXSize="1440" rasterYSize="400"> 
 <Geotransform>-180,0.25,0,50, 0, -0.25</Geotransform> 
 <SRS>WGS84</SRS>' >TRMM.vrt;
for i in 3B43.12*
do echo '<VRTRasterBand dataType="Float32" band="1" subClass="VRTRawRasterBand"> 
 <SourceFilename relativeToVRT="1">'$i'</SourceFilename> 
 <ByteOrder>MSB</ByteOrder>
 <ImageOffset>0</ImageOffset>
 <PixelOffset>4</PixelOffset>
 <LineOffset>5760</LineOffset>
 </VRTRasterBand>' >>TRMM.vrt
done;
echo '</VRTDataset>' >>TRMM.vrt

Within R, the script (trmm2012.sh) is executed and the virtual-file (TRMM.vrt) loaded as a rasterbrick. This is flipped in y direction and the files written as multilayer Geotiff. This Geotiff contains all the layers for 2012 and can be opened in every GIS software.

system("./trmm2012.sh")
 b=brick("TRMM.vrt")
 b
 trmm = flip(b, direction='y')
 writeRaster(trmm, filename="trmm_acc_2012.tif", format="GTiff", overwrite=TRUE)