Tag Archives: image processing

Notes on processing large files with GDAL

Context

I have recently had a project where I was provided with a reasonable number of large aerial image files. These were 600-800MB each and there were 10-15 of them. I needed to cut out the area of interest (which covered some part of all of the images), mosaic the images into a single representation of the site and find a way to share the data with non-specialists.

The workflow in my head was Clip – Mosaic – Webmap

Things to note

The images were too slow to render for me to mess about trying to do this using a GUI system, and the GDAL functions in QGIS3 (for me at least) don’t seem to be up and running properly (milage may vary).

Certain things that I have found when doing this in GDAL are that it’s best to define a data type for the bands, and to set the data structure as cloud optimised. The easiest way to create cloud optimised geotiffs is to use the following:

gdal_translate in.tif out.tif -co TILED=YES -co COPY_SRC_OVERVIEWS=YES -co COMPRESS=DEFLATE

Initially I was going to use gdal_merge.py to mosaic the clipped images as that is what is used in QGIS. However, gdal_merge.py reads everything into memory and I quickly found out that 16GB of RAM wasn’t enough as my swap partition started to be used and the whole thing ground to a halt.

The trick is to use the CPU if you can. First make a virtual raster using gdalbuildvrt, and then use gdal_translate (which utilises the CPU) to change the vrt file to whatever format you want (e.g. cloud optimised geotiff). This was actually really fast, and used the structure found in the following commands:

gdalbuildvrt output.vrt /path/to/folder/of/*.tif
gdal_translate -of GTiff output.vrt mosaic.tif

This ends up with an image that will be relatively large (1.2 GB in my case) which is still tricky to share. But fear not! We can use gdal2tiles.py to create a whole load of tiles and automatically generate some leaflet code to allow you to see the data. Then it’s just a question of moving that folder onto a web server and sending your contacts the link.

A quick thing to note/remember is that the webmap zoom factors relate to the following:

  • 0 represents the whole world (1:500,000,000)
  • 19 is very close up (1:1000)

Workflow

The workflow that I used in the end was as follows:

First crop each aerial image to the site shapefile

maindir="/path/to/folder/of/incoming/aerial/files"
cd ${maindir}

for f in *.tif

do
gdalwarp -ot Byte -of GTiff -cutline siteBoundary.shp 
-crop_to_cutline -dstnodata 0.0 -overwrite ${f} 
outputfolder/${f}_cropped.tif-co TILED=YES -co COPY_SRC_OVERVIEWS=YES
 -co COMPRESS=DEFLATE
done

Then build the virtual image

gdalbuildvrt output.vrt /path/to/folder/of/*.tif

Output the virtual image to a cloud optimised geotiff

gdal_translate -ot Byte -of GTiff -co TILED=YES 
-co COPY_SRC_OVERVIEWS=YES -co COMPRESS=DEFLATE output.vrt mosaic.tif

Create the web tiles and map file using the following command

gdal2tiles.py -s EPSG:27700 -z 11-19 mosaic.tif

Sources

https://lostingeospace.blogspot.com/2011/04/rapid-mosaicking-with-gdal.html

http://www.cogeo.org/

Advertisements
Tagged , , ,

Convert 16bit to 8bit GeoTiff

Note to self: use the following command to scale 16 bit spatial raster data down to 8 bit.

gdal_translate -a_nodata 0 -of GTiff -ot Byte -scale 0 65535 0 255 input.tif 8bit_output.tif

Tagged , , , ,

ARCSI

At the end of 2014 I used ARCSI (Atmospheric and Radiometric Correction of Satellite Imagery) for the first time. ARCSI is an open software project that provides a command line tool for the atmospheric correction of Earth Observation imagery. It provides a pretty much automatic way of running 6S.

Details of the software can be found here: https://bitbucket.org/petebunting/arcsi

At tutorial can be found here: https://spectraldifferences.wordpress.com/2014/05/27/arcsi/ 

This post is based largely on the tutorial linked to above but I also try to pull together some of the tips I read about in the help forums here: https://groups.google.com/forum/#!forum/rsgislib-support

Set up

All of the following instructions are run on an installation of the Ubuntu 14.04 operating system.

First, download the Anaconda 3.4 python distribution: http://docs.continuum.io/anaconda/

Install the Anaconda version of Python and the conda package manager using the command line. Once the bash script has run then install ARCSI and TuiView (a fast image viewer) using conda. If a warning comes up regarding gdal then install the gdal-data package.

bash Anaconda3-2.1.0-Linux-x86_64.sh

conda install -c https://conda.binstar.org/osgeo arcsi tuiview
conda update rsgislib arcsi

conda install -c jjhelmus gdal-data

For a successful installation I needed to change the default installation path for Anaconda from ~/anaconda3 to ~/anaconda

You should now be able to check your installation using the following command:


arcsi.py -h | less

To finalise the set up you need t point the GDAL driver path to the KEA drivers. Run the following to set up the path:


export GDAL_DRIVER_PATH=~/anaconda/gdalplugins:$GDAL_DRIVER_PATH

GDAL_DATA="/home/username/anaconda/share/gdal"

export GDAL_DATA

where username is the user account on the Ubuntu installation.

Run the code

When you run ARCSI, you’ll enter a command similar to the following:


arcsi.py -s ls8 -f KEA --stats -p RAD TOA SREF --aeropro NoAerosols --atmospro Tropical --aot 0.25 -o dir/to/outputs -i metadatafile_MTL.txt

To break this down a bit it first calls the arcsi.py script, passing to it the following parameters:

  • Sensor (-s) – landsat 8
  • Output image format (-f) – KEA
  • Parameters to compute (-p) – RAD Radiance conversion, TAO Top of atmosphere, SREF surface reflectance using 6S (for the full range please consult the official documentation)
  • Output directory (-o) – dir/to/outputs, into which the three computed outputs will be saved
  • Information (-i) – Landsat metadata file, using the format provided by the images available through Earth Explorer

 

Further information of all the parameters can be found in the help forums, the arcsi help command and in the code.

If an error is reported when you first run the ARCSI command, it might be that LIBGFORTRAN.SO.3 cannot be found. This will lead to the 6S model failing. You will need to install the required files using the following command:


sudo apt-get install libgfortran3

Dealing with outputs

The best way to view the output is to start the supplied viewer using the  ‘tuiview‘  command. This is an intuitive to use and very responsive viewer that handles the KEA format natively. To transform the KEA format outputs into a format readable by a wide range of GIS software, use the GDAL Translate command.


gdal_translate -of GTiff keafile.kea outputfile.tif

Tagged , , , ,

Cool image

Whilst whiling away my life on the web I came across this stunning image: http://eijournal.com/showcase-articles/atmospheric-lows-and-highs-create-catalina-eddy

It reminds me of why I love remote sensing, and in particular why satellite earth observation will never be beaten. Physical geography, I bow to your awesomness.

Tagged ,

Numpy

I have spent the last day or so looking at the interaction between PIL and numpy. PIL is useful for quickly loading and manipulating images (size, subsets etc) but if you really want to get to grips with image processing you need to understand numpy and the array structure. If you have come from Matlab then it can be a bit confusing as the syntax is “similar but different”. Anyway, here are some basic snippets of code:

import Image as I # Import PIL
import numpy as N # Import Numpy
import pylab as P # Import Matplotlib functions
a = I.open("/home/al/Documents/temp_image_proc/bird/image4.jpg", "r") # Open the image
b = N.asarray(a) # Convert to Numpy array
red = b[:,:,0] # Extract red band
P.imshow(red) # Display the image
P.show()

So we import the Image functions from PIL, numpy and the matplotlib plotting functions in pylab. Then we open up an image (.jpg format) convert it to an array, subset out the red band and display it. The image is shown here:

Continue reading

Tagged , , , , ,
Advertisements