Land Use Planning - Arumeru District -Tanzania

20090618

WORLDCLIM Data

All WORLDCLIM data for Tanzania (Zone 36 and 37) downloaded and imported in GRASS.

WorldClim is a set of global climate layers (climate grids) with a spatial resolution of a square kilometer. The current version is Version 1.4 (release 3).

Some notes: tiles are in Generic grid (raster) format. EPSG:4326. Spatial resolutions: 30 seconds (0.93 x 0.93 = 0.86 km²at the equator). Data files (.BIL) are sequential binary files in which values are stored line by line from the upper-left to the lower-right corner. Each cell (pixel) is an integer value (2 bytes; -32768..32767). If you simply open the data as images or raster data sets, there is a problem. ESRI software assumes that BIL files do not have negative values. These values (x) are replaced by (65535 - x); E.g., -10 becomes 65525. The nodata value of -9999 is not recognized... a BIG mess results... I got the same problem with GRASS. For a workaround:
Use DIVA-GIS (free download)
Data\Import to gridfile\Multiple Files (BIL/BIP/BSQ) and then
Data\Export to gridfile\Multiple Files (ESRI Ascii)
Now GRASS works, with module r.in.arc. Or Qgis GRASS plugin, r.in.arc module.

Informations are:
tmean = average monthly mean temperature (°C * 10)
tmin = average monthly minimum temperature (°C * 10)
tmax = average monthly maximum temperature (°C * 10)
prec = average monthly precipitation (mm)
bio = bioclimatic variables derived from the tmean, tmin, tmax and prec
alt = altitude (elevation above sea level) (m) (from SRTM)

Note that temperature is expressed as degree C * 10.

BIOCLIM
Bioclimatic variables are derived from the monthly temperature and rainfall values in order to generate more biologically meaningful variables. These are often used in ecological niche modeling (e.g., BIOCLIM, GARP). The bioclimatic variables represent annual trends (e.g., mean annual temperature, annual precipitation) seasonality (e.g., annual range in temperature and precipitation) and extreme or limiting environmental factors (e.g., temperature of the coldest and warmest month, and precipitation of the wet and dry quarters). A quarter is a period of three months (1/4 of the year).
They are coded as follows:
BIO1 = Annual Mean Temperature
BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
BIO3 = Isothermality (P2/P7) (* 100)
BIO4 = Temperature Seasonality (standard deviation *100)
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (P5-P6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO10 = Mean Temperature of Warmest Quarter
BIO11 = Mean Temperature of Coldest Quarter
BIO12 = Annual Precipitation
BIO13 = Precipitation of Wettest Month
BIO14 = Precipitation of Driest Month
BIO15 = Precipitation Seasonality (Coefficient of Variation)
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quarter

This scheme follows that of ANUCLIM, except that for temperature seasonality the standard deviation was used because a coefficient of variation does not make sense with temperatures between -1 and 1).

Links:
WORLDCLIM Docs:
http://www.worldclim.org/methods.htm
http://www.worldclim.org/format.htm
http://www.worldclim.org/bioclim.htm
DIVA-GIS
WORLDCLIM Data:
http://www.worldclim.org/tiles.php

AFRICOVER Data

All AFRICOVER data for Tanzania downloaded and imported in GRASS:

    •    tz_boundaries
    ⁃    tz_politicalboundary
    ⁃    tz_adm

    •    tz_towns
    ⁃    tz_majortowns
    ⁃    tz_othertowns

    •    tz_roads
    •    tz_rv
    •    tz_landform_agg
    •    tz_lithology_agg
    •    tz_spatial_agg
    •    tz_cultiv_agg
    •    tz_grass_agg
    •    tz_woody_agg

All vector in shape format, EPSG:4326. No problem with GRASS Qgis plugin, Module v.in.ogr.qgis.
Database uses sqlite driver, set up by GRASS Qgis plugin, Module db.connect.
AFRICOVER.sqlite created BEFORE importing data, with sqlite-manager extension for firefox.

Links:
AFRICOVER: http://www.africover.org/
Sqlite-Manager: http://code.google.com/p/sqlite-manager/
Mozilla Firefox: http://www.mozilla.com/en-US/

SRTM DEM

All SRTM tiles for Tanzania downloaded and imported in GRASS (EPSG:4326):

42_13 | 43_13 | 44_13
42_14 | 43_14 | 44_14
42_15 | 43_15 | 44_15 |45_15

All tiles patched by:
r.patch input="srtm_42_13,srtm_42_14,srtm_42_15,srtm_43_13,srtm_43_14,srtm_43_15,srtm_44_13,srtm_44_14,srtm_44_15,srtm_45_15" output="srtm_full"

Links:
SRTM (Original dataset Docs): http://www2.jpl.nasa.gov/srtm/index.html
SRTM (V4): http://srtm.csi.cgiar.org/

Updates

Switched to:

GRASS 6.4 RC5
Qgis 1.1.0

QGIS/GRASS and Frameworks Mac binaries: http://www.kyngchaos.com/software:unixport

20090331

Reprojecting data

Another map has been were georeferenced in WGS84LL using Georeferencing Qgis Plugin:

"Agroecological zones of Tanzania", Ministry Of Agriculture, Food and Cooperative

Some reprojection has been done on raster maps using gdal tools:

No Datum change:

gdalwarp -t_srs EPSG:32737 4326_HRImage.tif 32737_HRImage.tif

move the source WGS84LL source GeoTiff file to the WGS84UTM37 huge, uncompressed target file, so you need to reduce it with

gdal_translate -co "COMPRESS=JPEG" -co "PHOTOMETRIC=YCBCR" -co "JPEG_QUALITY=100" 32737_HRImage.tif 32737_HRImage2.tif

and add pyramid overviews (load very fast in Qgis) with

gdaladdo --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL 32737_HRImage2.tif 2 4 8 16

Datum change:

gdalwarp -t_srs EPSG:32737 -s_srs "1285_TZ.prf" 21037_Topo.tif 32737_Topo.tif

move the source ARC60UTM37 source GeoTiff file to the WGS84UTM37 file, using the right transformation parameters written in:

1285_TZ.prf
+proj=utm +zone=37 +south +a=6378249.145 +b=6356514.96582849 +units=m +towgs84=-175,-23,-303

The process works backward: from standard EPSG_code --> to custom_CRS

In this case you can add right info to the header of the result with:

gdal_translate -a_srs "21037_WK.txt" source.tif 21037_target.tif

with 21037_WK.txt:

PROJCS[" Projection Name = Arc_1960_UTM_Zone_37S Units = meters GeoTIFF Units = meters",
    GEOGCS["Arc 1960",
        DATUM["Arc_1960",
            SPHEROID["Clarke 1880 (RGS)",6378249.145,293.4649999999983,
                AUTHORITY["EPSG","7012"]],
            AUTHORITY["EPSG","6210"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4210"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",39],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

that can be extract with gdalinfo>info.txt on the original toposheets (see Defining the Area of Intervention)

In this way we obtained:

./Raster/21037:
21037_Agro_Eco_Zones.tif
21037_Arusha_soil.tif
21037_EcoClimatic_Zones.tif
21037_HRImage.tif
21037_Topo.tif

The same structure is repeated for 32737 and 4326.

With

gdaltindex index.shp *.tif

we obtained the coverage of all the rasters in every folder.

The same transformation procedures work with vector data:

ogr2ogr 4326_Index.shp 21037_Index.shp -t_srs EPSG:4326 -s_srs "1285.prf"

ogr2ogr mosaic2.shp mosaic.shp -t_srs EPSG:4326 -s_srs "+init=EPSG:21037 +wgs84=-175,-23,-303"

Links:
http://www.gdal.org/gdalwarp.html

20090327

Georeferencing raster data

Two maps:

"Soils of Arusha and Monduli Area 1:250.000" (Soil Research Institute, Ottawa, Canada, 1974)
"Tanzania Eco-Climatic Zones" (Institute of Resource Assestment).

were georeferenced in WGS84LL using Georeferencing Qgis Plugin.

Other sources can be georeferenced using gdal utilities.
For testing purposes is possible to download high resolution imagery covering part of Area of Intervention from Google dataset using GoogleMV. Level 17 has enough resolution.
The result is a big *.jpg file and a *.map file (Ozi Explorer format, see here info on the format).
The *.map is simple text, and coordinates of the corners of the image in latlongwgs84 can be found near the end of the file:

MMPLL,1, 36.754167, -3.076833
MMPLL,2, 36.903752, -3.076833
MMPLL,3, 36.903752, -3.233316
MMPLL,4, 36.754167, -3.233316

(in red are marked the Upper Left and Lower Right corners coordinates).

gdal_translate can be used now:

gdal_translate -a_srs EPSG:4326 -a_ullr 36.754167 -3.076833 36.903752 -3.233316 -co "COMPRESS=JPEG" -co "PHOTOMETRIC=YCBCR" -co "JPEG_QUALITY=100" source.jpg target.tif

the -a_srs parameter assign an EPSG code to the image
the -a_ullr parameter define upper left and lower right corners coordinates (taken from the *.map file)
the -co parameters define the compression algoritm used

The result is a GeoTIFF with this header (read it with gdalinfo):

Driver: GTiff/GeoTIFF
Files: target.tif
Size is 13943, 14608
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235630016,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (36.754167000000002,-3.076833000000000)
Pixel Size = (0.000010728322456,-0.000010712144031)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
SOURCE_COLOR_SPACE=YCbCr
COMPRESSION=JPEG
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 36.7541670, -3.0768330) ( 36d45'15.00"E, 3d 4'36.60"S)
Lower Left ( 36.7541670, -3.2333160) ( 36d45'15.00"E, 3d13'59.94"S)
Upper Right ( 36.9037520, -3.0768330) ( 36d54'13.51"E, 3d 4'36.60"S)
Lower Right ( 36.9037520, -3.2333160) ( 36d54'13.51"E, 3d13'59.94"S)
Center      ( 36.8289595, -3.1550745) ( 36d49'44.25"E, 3d 9'18.27"S)
Band 1 Block=13943x16 Type=Byte, ColorInterp=Red
Band 2 Block=13943x16 Type=Byte, ColorInterp=Green
Band 3 Block=13943x16 Type=Byte, ColorInterp=Blue

The file is huge, so pyramid layers can be added to fasten the opening in GIS with gdaladdo:

gdaladdo --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL target.tif 2 4 8 16

In this way the last lines of the header of GeoTIFF are substituted with:

Band 1 Block=13943x16 Type=Byte, ColorInterp=Red
Overviews: 6972x7304, 3486x3652, 1743x1826, 872x913
Band 2 Block=13943x16 Type=Byte, ColorInterp=Green
Overviews: 6972x7304, 3486x3652, 1743x1826, 872x913
Band 3 Block=13943x16 Type=Byte, ColorInterp=Blue
Overviews: 6972x7304, 3486x3652, 1743x1826, 872x913

You can obtain a shapefile with the coverage of the GeoTiff with:

gdaltindex target.shp target.tif

Links:

GoogleMV:

http://www.silber2004.narod.ru/trevel.htm (this link is in Russian, you can download the program here)

GDAL Utilities:

http://www.gdal.org/gdal_utilities.html

OZI *.map format:

http://www.rus-roads.ru/gps/help_ozi/map_file_format.html

20090325

Defining the Area of Intervention

The Project AI (Area of Intervention) is defined by four of the Standard Topo Sheets "East Africa 1:50.000 (United Republic of Tanzania)":

55/1 Oldonyo Sambu
55/2 Engare Nanyuki
55/3 Arusha
55/4 Usa River.

These maps (geotiff) are georeferenced in Standard Arc 1960 UTM (EPSG:21037).

To obtain a single raster mosaic covering all the area:

ls -1 *.tif > tiff_list.txt
gdal_merge.py -v -o mosaic.tif --optfile tiff_list.txt

Then with:

gdalinfo mosaic.tif

You can obtain:

Driver: GTiff/GeoTIFF
Files: mosaic.tif
Size is 13109, 13042
Coordinate System is:
PROJCS[" Projection Name = Arc_1960_UTM_Zone_37S Units = meters GeoTIFF Units = meters",
    GEOGCS["Arc 1960",
        DATUM["Arc_1960",
            SPHEROID["Clarke 1880 (RGS)",6378249.145,293.4649999999983,
                AUTHORITY["EPSG","7012"]],
            AUTHORITY["EPSG","6210"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4210"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",39],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (222098.049999999988358,9668234.250000000000000)
Pixel Size = (4.250000000000000,-4.250000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 222098.050, 9668234.250) ( 36d30'0.01"E, 2d59'56.28"S)
Lower Left ( 222098.050, 9612805.750) ( 36d29'55.57"E, 3d29'59.98"S)
Upper Right ( 277811.300, 9668234.250) ( 37d 0'3.48"E, 3d 0'0.01"S)
Lower Right ( 277811.300, 9612805.750) ( 36d59'59.92"E, 3d30'4.33"S)
Center      ( 249954.675, 9640520.000) ( 36d44'59.74"E, 3d15'0.26"S)
Band 1 Block=13109x1 Type=Byte, ColorInterp=Red
Band 2 Block=13109x1 Type=Byte, ColorInterp=Green
Band 3 Block=13109x1 Type=Byte, ColorInterp=Blue

These are the limits (and other useful information) of the AI in Arc60, UTM and LatLong.

You can obtain a shapefile covering AI (EPSG:21037)with:

gdaltindex mosaic.shp mosaic.tif

Now you can add a *.prj file to the shape with:

ogr2ogr mosaic2.shp mosaic.shp -a_srs EPSG:21037

And read parameters with:

ogrinfo -al -so mosaic2.shp

Links:

GDAL:

OGR: