20090331

Reprojecting data

Another map has been were georeferenced in WGS84LL using Georeferencing Qgis Plugin:
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"
or
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



Powered by ScribeFire.

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:
GDAL Utilities:
OZI *.map format:




Powered by ScribeFire.

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:
  1. ls -1 *.tif > tiff_list.txt
  2. 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:




Powered by ScribeFire.

20090324

EPSG Codes and GRASS LOCATIONS

Software used:

Mac OS X 10.5.6 on Intel Macbook (late 2007)
Quantum GIS (1.0.1-Kore)
GRASS GIS (6.4.0RC3)


CRS used:
  • EPSG:4326 - WGS84LL
  • EPSG:32737 - WGS84UTM37S
  • EPSG:21037 - ARC60UTM37S
  • EPSG:4210 - ARC60LL
  • Custom ARC60UTM37S (Kenya_Tanzania, based on EPSG:1122): +proj=utm +zone=37 +south +ellps=clrk80 +units=m +towgs84=-160,-6,-302
  • Custom ARC60UTM37S (Tanzania Onshore and Offshore, based on EPSG:1285): +proj=utm +zone=37 +south +a=6378249.145 +b=6356514.96582849 +units=m +towgs84=-175,-23,-303
Projection Parameters (as written in the Standard Tanzanian Topo Sheets):

  • Grid: UTM Zone 37
  • Projection: Transverse Mercator
  • Spheroid: Clarke 1880 (Modified)
  • Unit of Measurement: metre
  • Meridian of Origin: 39°00' East of Greenwich
  • Latitude of Origin: Equator
  • Scale Factor at Origin: 0.9996
  • False Coords of Origin: 500.000 Easting; 10.000.000 Northing
  • Datum: New (1960) Arc
GRASS Locations

Standard EPSG 4326, 21037, 32737

  1. In QGIS/GRASS Plugin Create new mapset
  2. Choose the database (path_to_GIS_data)
  3. Create New Location: WGS84LL ("Projected": EPSG:4326)*
  4. Default Wind: Tanzania
  5. Choose a name for the Mapset (ARUMERU)
* Repeat for WGS84UTM37S (EPSG:32737) and for ARC60UTM37S (EPSG:21037)

Links:
QGIS: http://www.qgis.org/
GRASS GIS: http://grass.itc.it/
QGIS/GRASS and Frameworks Mac binaries: http://www.kyngchaos.com/software:unixport
EPSG Codes: http://www.epsg-registry.org/

Powered by ScribeFire.



Village Land Demarcation

Village land demarcation activity was done to define the village administrative boundaries. This is important step before starting the village land use planning exercise. The challenges of village Land demarcation activities were immense. The land issue in the Meru and Maasai communities is very sensitive due to complexity and historical issues related to land tenure. The villagers decided the boundaries of the Villages to follow the colonial settler farms that were in the old map projections (TTM), so we had to re-project the geographical data using free GIS software to UTM ARC 60. We had to seek for data related to the boundaries of the settler farms, West Kilimanjaro Ranch, Arusha National Park and the Meru District Council boundaries as stipulated in the Government Notice (GN). Getting all this old data was a big challenge and it takes a lot of time and other resources. After all the challenges and with right projections of the coordinates, the beacons were set on the village boundary points. The village leaders, village land adjudication committees, the district surveyor and representatives from institutions such as West Kilimanjaro Ranch were involved. out of this exercise, Village maps were produced.

Legal and Institutional issues related to Village Land Demarcation.
- Village Councils
- Village land adjudication committees
- Village Assembly
- District Council
- Ministry of Land and human Settlements Development
Important documents
- National Land Policy, 1995
- Village Land Act, 1999
- Village Land Act Regulations, 2002
- Village Land Act Guidelines developed by ministry of land.
Important trainings
- Stakeholders workshop - target (all stakeholders)
- Village land Act, 1999 - Target (district staff, ward, villages)
- Village boundary agreements procedures and processes -Target ( Village councils, Village Land Committees)
- Basic GIS training - Target (district council staff)
- Compass use and basic Map reading - Target (village game scouts [VGS]), village leaders

20090323

Helo world

Here I (and the rest of the MAE Team, I hope) will publish progresses, constraints, challenges on project implementation.

At least, I'll try to do it...

Rock
and
Roll