Landsat 8 example using pan-sharpening in Orfeo Toolbox

February this year USGS has sent another satellite up into orbit to continue observing the earth. This satellite, named Landsat 8, provides the same spatial resolution in all bands as Landsat 7. That is one pixel equals 15m pan-chromatic and 30m for all other bands (excluding thermal). In addition, it has a couple extra bands for water, cloud and surface temperature (this link explains the other similarities and differences in more detail).

Here is a sample that I have pan-sharpened using the open source Orfeo Toolbox. Bands displayed are 6, 5 & 2, which is comparable to 5, 4 & 1 in Landsat 5 & 7. This image was capture 27th of May 2013. Most of the bright green paddocks are canola and you can see some of the earlier cereals coming through.

Bands 6, 5 & 2. Extract from scene Path 91, Row 80 on  27-May-2013
Landsat 8: Bands 6, 5 & 2. Extract from scene Path 91, Row 80 on 27-May-2013

More info: As suggested in comments, Spectral Transformer for Landsat-8 imagery by GeoSage is another free option for pan-sharpening Landsat 8 imagery. I have tried it and it works well if you don’t mind following some simple command line instructions. Check it out at this link.

Even more info: To get an idea of the best Landsat 8 band combinations and comparison with Landsat 7 bands check out this blog article from ESRI.

 
http://rcm-na.amazon-adsystem.com/e/cm?t=agma-20&o=1&p=8&l=as1&asins=1439845379&ref=tf_til&fc1=000000&IS2=1&lt1=_blank&m=amazon&lc1=0000FF&bc1=FFFFFF&bg1=FFFFFF&f=ifr

Using MultiSpec to composite Landsat images

As most GIS users will know the Landsat archive of satellite images is available for free download from USGS. This is a valuable resourse of historic and recent course resolution images covering most of the earth. I find the USGS Global Visualization Viewer the best way to sort through and download images. You need to add images to the basket. Often you can download them directly from the basket, but sometimes they need processing time which can take a day or two. Note that images from L7 2003 onwards have lines of no data caused by a faulty sensor. Up until recently L5 was still working producing quality images but is now out of service. Shortly we will have Landsat 8 – it has been launched and more free imagery is not far off!

Once the image file is downloaded you will have a compressed .tar.gz file which can extracted using built in Windows functions or 7Zip. Once extracted you have a TIFF (.tif) file for each of the bands of the image. For L7, that is 9 TIFF files about 60mb each except band 8 which is higher resolution band for pan-sharpening usually 200mb+.

The next step is to composite these bands to make a coloured image to load into your GIS. If you are a Windows or Mac user I suggest you try out MultiSpec. Download and extract to a folder of you choice. I keep it at C:MultiSpec for simplicity.

A comprehensive guide on how to composite images has been written up by the creators of MultiSpec and is available here. The process I give below is a simplified explanation based on the MultiSpec tutorial.

1. Launch MultiSpec.

2. Go to File > Open and select all bands you want to composite into one image.

Open multiple files
Open multiple files

3. For the ‘Set Display Specifications’ dialogue press OK (Try setting Bands to Red: 5, Green: 4 and Blue: 1. This is not critical for this exercise but is useful to visually inspect the image inside MultiSpec).

4. Press Ctrl + R to open ‘Set Image File Format Change Specifications’

Set Image File Format Change Sepecifications
Set Image File Format Change Specifications

5. You can select Channels and reorder by choosing ‘Subset…’ (I usually include bands 1-5) and also select file format in the ‘Header’ drop box (I always stick with GeoTIFF format). Press OK and choose where you want your composite image to reside and give it a useful name (do not forget to add the .tif extension).

6. Load up your favorite GIS desktop software, QGIS for me, and add your new TIFF file.

7. Edit layer properties and try a few different band combinations and contrast enhancements.
For good visibility of vegetation try Red: Band 5, Green: Band 4 and Blue: Band 1.
For a more natural look try Red: Band 3, Green: Band 2 and Blue: Band 1.

Try different band combinations and contrast enhancement
Try different band combinations and contrast enhancement

I hope this helps make Landsat imagery more accessible to everyone.

Publishing GIS data on the web: Lizmap

Lizmap is Open Source software that has the capacity to display GIS data on the web with minimal effort once it has been configured correctly. This software sets itself apart by making use of the QGIS Server in a unique and user friendly way. In essence this allows the user to configure a map in QGIS desktop software and simply run a basic plugin to publish their data on the web with a Google Maps satellite or street layer as a base map.

I have setup Lizmap software on my server to demo its capabilities. I left the original demo map of transport in Montpellier, France and added my own with data from Lindon – one of our properties. I hope to demo its mobile device capabilities soon.

Check out my Lizmap demo here

 

Lizmap Web Client displaying Lindon property map.
Lizmap Web Client displaying Lindon property map.

More information about Lizmap:
Lizmap Web Client Documentation
Lizmap QGIS Plugin Documentation

http://rcm-na.amazon-adsystem.com/e/cm?t=agma-20&o=1&p=8&l=as1&asins=0596008651&ref=tf_til&fc1=000000&IS2=1&lt1=_blank&m=amazon&lc1=0000FF&bc1=FFFFFF&bg1=FFFFFF&f=ifr

Combining Yield Data

Yield data is collected from most modern harvesters. This data can be used to analyse a crops individual performance, compare with previous crops and combined with several years of data to generate statistical maps. In this post I show you three individual years yield maps for one paddock: wheat, sorghum and chickpeas. Using this three seasons of data I then ‘stack’ them on top of each other to produce maximum, minimum and mean. These statistics can then be used as a starting point to develop paddock zones to apply variable rate fertiliser or seed.

MODIS through Google Earth Engine to monitor crops

What Google is doing with the earth’s satellite imagery is really quite amazing. Check out Google’s Earth Engine. On the landing page you get some great videos. I recommend Amazon Deforestation: Timelapse and Drying of the Aral Seas: Timelapse.

Once you have had enough of these visit their data catalog where you can load any of the common earth observing satellites into a Google Workspace which is the Google Map environment with additional features allowing you to display the satellite images.

Google say they have developed their Earth Engine mainly for monitoring deforestation which is valuable but there are several ways this technology is useful in agriculture. For example load the MODIS Daily EVI into the Workspace. Once in the workspace go to a farming area you know well and you can monitor how much (healthy) crop is established. Although MODIS imagery is free to download, Google take out all the processing time and make it easy to keep an eye on the district.

MODIS EVI captured on 2012-12-05. Showing area north of Moree NSW Australia. Google Earth Engine used to display image.
MODIS EVI captured on 2012-12-05. Showing area north of Moree NSW Australia. Google Earth Engine used to display image.

QGIS SEXTANTE Modeler: Automate yield data processing

This is a followup post to Processing harvest yield data to grid / raster where I outlined my yield data processing steps. This post shows how I use the Modeler in SEXTANTE to automate this process utilising SAGA algorithms.

The idea of the modeler for me is to save time and keep consistent processes. In the screen shot below you can see that the model requires the yield data point file (where the Dry Yield table column needs selecting), paddock boundary and the grid and filtered points output file location. The model takes care of all the processing steps.

Select yield data point file to be processed and choose output location. The modeler does the rest!
Select yield data point file to be processed and choose output location. The modeler does the rest!

While constructing the model the only tricky part was making sure I had a separate ‘Extent’ field to input into the ‘Shapes to Grid’ algorithm. Even though I always leave this blank, it is required for the algorithm to realise the extent is just the minimum covering extent of the input shape file.

SEXTANTE model editing screen displaying input fields and algorithms involved in the yield data processing
SEXTANTE model editing screen displaying input fields and algorithms involved in the yield data processing
Yield Map: L3 Chickpeas 2012
Yield Map: L3 Chickpeas 2012

A few years ago I used to use model builder in ArcGIS and really missed it when I no longer had access to ESRI software. The SEXTANTE Modeler fills that gap and allows a combination of algorithms from several different toolboxes.

Some more advanced applications of Sextante Modeler are available here.

Processing harvest yield data to grid / raster

The process to convert point data to raster (or grid) for can differ dramatically depending on who you ask and the purpose of the conversion. Generally data is processed for two main functions. The first being ‘stacking’ or ‘layering’ to later develop prescription maps or other further processing. The second is improved visual representation of the data. My process aims to satisfy both with minimal processing time. My method does contain some compromises; since I use a coarse resolution it is not as pleasing on the eye as some other methods. In addition, my data smoothing technique is quite broad which means you will loose some detail in the grid.

At this point in time I follow these steps:

  1. Clean up point data either manually manually or with a filter to remove any obvious errors such as where header turns in paddock. EDIT: I am currently trying Yield Editor for this step.
  2. Produce a grid over the top of the point data. Any cells that share a points average the point values.
  3. Gaps in the grid are filled in.
  4. Gaussian Filter is used to ‘smooth’ the data.

In QGIS this is the process:

Note: You will need to have SEXTANTE and additional toolboxes setup to follow these steps. Instructions are available here.

  1. Load paddock boundary (shape file polygon)
  2. Load yield data points (You will need these in shape file format – use export function in SMS or FOViewer if your yield monitor does not produce ESRI shape file)
  3. SEXTANTE > SAGA Toolbox > Shapes – Points > Points Filter
    1. Radius: 100
    2. Minimum Number of Points: 25
    3. Maximum Number of Points: 200
    4. Quadrants: No
    5. Filter Criterion: Remove Below Percentile
    6. Percentile: 15-20
  4. SEXTANTE > SAGA Toolbox > Shapes – Points > Points Filter (use output from step 3)
    1. Radius: 100
    2. Minimum Number of Points: 25
    3. Maximum Number of Points: 200
    4. Quadrants: No
    5. Filter Criterion: Remove Above Percentile
    6. Percentile: 90
  5. SEXTANTE > SAGA Toolbox > Grid – Gridding > Shapes to Grid
    1. Preferred Target Grid Type: Floating Point
    2. Cell Size: 15
    3. Method for Multiple Values: Mean
  6. SEXTANTE > SAGA Toolbox > Grid – Tools > Close Gaps
  7. SEXTANTE > SAGA Toolbox > Shapes – Grid > Clip Grid with Polygon
  8. SEXTANTE > SAGA Toolbox > Gaussian Filter
    1. Standard Deviation: 3
    2. Search Mode: Circle
    3. Search Radius: 50

Make sure you do a visual comparison with original point data to ensure the final grid is a good representation of the original data.

Visually check to see that processed grid file is a fair representation of the original point data.
Visually check to see that processed grid file is a fair representation of the original point data. Note I was lazy in this photo and used un-filtered yield data.

To generate useful color maps to represent the grid data try http://colorbrewer2.org/. One day I will explain this process, but it is not too hard to try it yourself.

Now when we get clever we can use the SEXTANTE Modeler to automate this process. More on this later! –Edit: Here is the modeler info.