Adding rainfall

In the previous post, we produced a polygonal area of Central America within a range of elevation and latitude.

Now we want to limit to areas with roughly the same annual rainfall.  I've hunted around for a public global annual precipitation dataset. It's actually pretty hard to find.  Found:

1. World Bank Average precipitation is annual, but only per country.
2. WorldClim - Global Climate Data, which is broken down by month, not annual.
3. NASA NEO Rainfall TRMM is also only available in monthly form.
4. ArcGIS Total Annual Precipitation is annual, which is what we need, but unfortunately doesn't let you download it.  It's based on WorldClim.

Well, that's frustrating. I also tried asking on Twitter but the gistribe had no leads.

Then, it occurred to me that I could try adding up all 12 months myself, using "raster arithmetic".  I did the following:

1. From WorldClim 1.4: Current conditions: download by tile, I clicked on tiles "22" and "23" and downloaded the GeoTiff precipitation zipfile for each.
2. For a tile, unzipping the zip gives 12 files of the form precN_T.tif, where N is the month 1-12, and T is the tile number.
3. I dropped all 12 months into QGIS, and did some examination of the values.  As expected, each pixel in each layer has an integer number of mm of rainfall that fall there in the given month of the year:
   
4. QGIS command: Raster: Raster Calculator. I entered an expression that is the sum of all 12 layers:
   
5. That only takes a few seconds to produce a new layer, I called it "precip_annual_tile22.tif", which contains the rainfall for the whole year!  I used "Identify" to manually check a few spots, and the amount of rain looked believable, ranging from less than 10cm in the desert of Baja California, to more than 4 meters in the mountains of northern Chiapas.
6. I did the same thing for tile 23.  Now all of Central America is summed up:
   

Now that we have total annual rainfall, I'd like to identify the areas which have rainfall similar to my part of Hawaii, which is roughly 120-160 inches (3-4 meters).  We can do the same thing we did for elevation, using GRASS (in QGIS: Processing Toolbox>Grass>Raster>r.recode).  This time the "recode rules" are a tiny file containing "3000:4000:1:1".  The region extent is the same as before: -120,-60,15,25

• 

The result is... surprisingly small again.  Only a few remote mountainous areas have as much rainfall as my part of Hawaii!  I decided to relax the range a little, from 2500 to 4500 mm of rainfall, and ran it on both tiles.  I made them purple (Layer Properties: Style: Singleband pseudocolor: Add value for 1) to make them easier to see against the greyscale rainfall:

That's a little better.  I want to be able to visit people in these places, so I hope it's not all remote wilderness areas!  It's basically some mountains of Puebla, Veracruz, Oaxaca, Chiapas, parts of northern Guatemala and Honduras, and small spots in the Caribbean islands of Jamaica, Haiti, and Puerto Rico.

The next step turns these rasters into polygonal vector layers (again with gdal_polygonize) and join them with Union.  Now we can see the result along with the elevation range polygons we computed before, set to 50% transparent, and zooming into part of Mexico so you can see where they overlap:

Finally, we want the intersection of these two polygonal layers.

1. I try "Vector: Geoprocessing: Intersection" and give my two layers as input.  Unfortunately, I get: Algorithm Intersection starting... Input layer A contains invalid geometries (feature 43). Unable to complete intersection algorithm.
2. I press "Close" and try again, with the layers in the other order.  This time, Intersection freezes, which makes QGIS lock up, forcing me to kill and restart the application.

Well, that's less than friendly.  It could at least have told us what kind of "invalid geometries" it found?  As far as I know, these are simply polygonal layers produced by gdal_polygonize and Union.  They shouldn't have any bad geometry in them... and applications shouldn't hang.

Since Intersection seems buggy, perhaps we can at least view these two polygonal layers on top of a real world map, as separate layers.  I'll tackle that in my next post.