Where did all the trees go?

I had a conversation with the VSO programme director the other day about the environment in Jeneponto district where I'll be working. His best guess was that there isn't any forest left in Jeneponto - which may make my job as a forestry mapping advisor quite difficult, or ridiculously easy, depending on which way you look at it.

I'm travelling to Jeneponto at the end of next week, but it would be useful to have some background into what the area looks like before I get there. I've been playing around with MODIS imagery (http://rapidfire.sci.gsfc.nasa.gov/) with little success, and this week decided to have a look at Landsat data instead.

A quick not-too-techy guide to the differences. MODIS (Moderate Resolution Imaging Spectrometer) has 30 something spectral bands, a 16 day repeat pass (if you use both satellites), but a best resolution of 250m. The most recent Landsat has7 bands at a resolution of 30m, plus a panchromatic band at 15m. The number of bands is essentially the number of "colours" the satellite can see, so MODIS is good for big things where you want to know the colour precisely, such as for ocean temperature, atmospheric aerosols, and sea ice. Landsat is better at seeing smaller things - and as some of the plantations and forest clearings may be only a few hundred metres across, this may be useful.

Landsat data can be downloaded here http://glcfapp.glcf.umd.edu:8080/esdi/index.jsp - this is the Global Landcover Survey, so only archives a few cloud free scenes rather than everything acquired by the satellite.

Each band is downloaded separately as a geotiff, so a bit of processing is required to turn these into a multiband image. I've done this in python using the gdal bindings:

from osgeo import osr, gdal 

src1=gdal.Open("L71114064_06420060822_B10.TIF") 

src2=gdal.Open("L71114064_06420060822_B20.TIF") 

src3=gdal.Open("L71114064_06420060822_B30.TIF") 

 

band1=src1.GetRasterBand(1).ReadAsArray() 

band2=src2.GetRasterBand(1).ReadAsArray() 

band3=src3.GetRasterBand(1).ReadAsArray() 

drv=gdal.GetDriverByName("GTiff") 

dst=drv.Create("Bands321.tiff",src1.RasterXSize,src1.RasterYSize,3,gdal.GDT_Byte) 

srs=osr.SpatialReference() 

srs.SetWellKnownGeogCS("WGS84") 

srs.SetUTM(50,0) 

dst.SetProjection(srs.ExportToWkt()) 

dst.SetGeoTransform(src1.GetGeoTransform()) 

dst.GetRasterBand(1).WriteArray(band3) 

dst.GetRasterBand(2).WriteArray(band2) 

dst.GetRasterBand(3).WriteArray(band1) 

dst.FlushCache()

As in other things, it's very important to flush before you finish.

[I'm going to sort out code highlighting soon, honest...]

Finally I wanted a vector land-sea boundary. This is picked out well in the infrared Landsat ETM+ band 5. Using command line gdal:

gdal_contour -fl 30 L71114064_06420060822_B50.TIF contour.shp

This generates a vector shoreline, after a bit of manual editing.

The results look OK. You can see the city of Makasar in the top left, and the scar left by the collapse of the Mount Bawakaraeng crater (centre right of image). I'm now playing with different infrared band combinations to see which best picks out forest and other vegetation. I'm also trying to merge the panchromatic band (15m resolution) to sharpen up the colour images (30m resolution). More to follow...