Category Archives: Uncategorized

Converting NEODAAS Mercator Projection netCDF files to GeoTiffs for use in QGIS / ArcMap

NetCDF files are a common format for distributing Earth Observation data and allow the ability to store a number of variables alongside metadata. However, using netCDF files in a GIS is not always as easy as it could be.

The NERC Earth Observation Data Acquisition and Analysis Service (NEODAAS) routinely produce products such as Chlorophyll from EO data and store as netCDF files. For the UK they use a Mercator projection within a netCDF file storing the latitude and longitude of each pixel within separate arrays. Unfortunately QGIS and ArcMap are often unable to read this information so don’t read data into the correct location making it difficult to use with other datasets.

To read data into the correct location I wrote a script which converts the latitude and longitude values in the netCDF file into tie points and then uses these to warp a GeoTiff into the correct location. It exports a single variable at a time.

To use for creating a GeoTiff from Chlorophyll data:

python reproject_neodaas_netcdf.py \
   sentinel3a_olci_all_products_L3 median_uk_7d_20180720_20180726.nc \
   sentinel3a_olci_all_products_L3-median_uk_7d_20180720_20180726_chl.tif \
   --variable CHL_OC4ME

The script is below. It requires the GDAL and netCDF Python libraries.

Note although the files should line up with other datasets you are warping the data which may introduce some errors. For best results, if you have a research project funded by NERC or are eligible for a NERC research grant or training award you can contact NEODAAS to discuss specific processing requirements.

Introduction to RSGISLib Training Course

A course introducing RSGISLib was recently given in Japan by Pete Bunting. Material from the course is available to download from the following links:

The course covers pre-processing of PALSAR SAR data followed by a demonstration of pixel-based and object-based classification. All the required data and scripts are included in the package.

If you have any questions on the course please use the RSGISLib support Google group: https://groups.google.com/forum/#!forum/rsgislib-support

Convert LaTeX to Word

A big problem with writing in LaTeX is collaborating with colleagues who don’t use it. One option is to generate a Word .docx version and use the comments and track changes features in Word / LibreOffice. This does require manually copying the changes back to LaTeX so isn’t quite as nice as using latexdiff (see earlier post) but is slightly easier than adding comments to a PDF.

The best program I’ve found for converting LaTeX to Word is the open source (GPL) command line tool, pandoc (http://pandoc.org/).

Basic usage is quite straightforward:

pandoc latex_document.tex -o latex_document_word_version.docx

The conversion isn’t perfect, figures and tables can get a bit mangled, but it does a good job with the text.

Pandoc can convert between many different formats, including from markdown and reStructuredText (commonly used for software documentation) so it is worth having installed.

Create a CSV with the coordinates from geotagged photos

Recently I took a lot of photos on my phone during fieldwork (survey for the NERC-ARF LiDAR and hyperspectral calibration flights) and wanted to extract the GPS coordinates from each photo so I could load them into QGIS using the add delimited text dialogue. The aim wasn’t to have precise locations for each photo (GPS positions for each of the points surveyed were recorded separately) but to give a quick idea of the locations we’d visited before the main GPS data were processed.

I while ago (2012!) I wrote a script for this task. The script (CreateJPEGKMZ) requires pillow and the imagemagick command line tools.

It code isn’t particularly tidy (despite my recent updates) but it basically pulls out the GPS location from the EXIF tags and writes this to a CSV file. To create a KMZ a thumbnail image is created and a corresponding KML file. Both the thumbnail and KML are then zipped together to generate the KMZ which can be opened in GoogleEarth.

To use the script to write a CSV file:

CreateJPEGKMZ.py --outcsv gps_points.csv input_jpeg_files

To also create KMZ files for each photo:

CreateJPEGKMZ.py --outcsv gps_points.csv \
                 --outkmz output_kmz_files 
                 input_jpeg_files

This was also a good lesson in the benefits of having scripts in version control and publicly available.

Convert EASE-2 grid cell to latitude and longitude using Python

The EASE-2 grid is used for a number of NASA datasets including SMAP. It is described in the following paper:

Brodzik, M. J., Billingsley, B., Haran, T., Raup, B., & Savoie, M. H. (2012). EASE-Grid 2.0: Incremental but Significant Improvements for Earth-Gridded Data Sets. ISPRS International Journal of Geo-Information, 1(3), 32–45. http://doi.org/10.3390/ijgi1010032

Files with the centre coordinate of each cell for EASE-2 grids at different resolutions are available from https://nsidc.org/data/ease/tools.html as well as tools for conversion.

To read these files into Python the following steps can be used:

  1. Download the relevant gridloc file.

    The FTP link for the grid location files is: ftp://sidads.colorado.edu/pub/tools/easegrid2/

    For this example I’ve chosen the 36km cylindrical EASE-2 grid (gridloc.EASE2_M36km.tgz)

  2. Un-tar using:
    tar -xvf gridloc.EASE2_M36km.tgz
    
  3. Read the files into Python:
    import numpy
    
    # Read binary files and reshape to correct size
    # The number of rows and columns are in the file name
    lats = numpy.fromfile('EASE2_M36km.lats.964x406x1.double', 
                          dtype=numpy.float64).reshape((406,964))
    lons = numpy.fromfile('EASE2_M36km.lons.964x406x1.double', 
                          dtype=numpy.float64).reshape((406,964))
    
    # Extract latitude and longitude
    # for a given row and column 
    grid_row = 46
    grid_column = 470
    
    lat_val = lats[grid_row, grid_column]
    lon_val = lons[grid_row, grid_column]
    

Add an ENVI header to JAXA Global Mangrove Watch PALSAR tiles

This is an update to a script I wrote a while ago for working with PALSAR tiles from JAXA’s Global Mangrove Watch.

Untar files
Using GNU parallel [1] you can untar a directory of files at once:

ls *tar.gz | parallel -j 2 tar -xf 

Which will untar two files at a time (depending on how fast your disk is you may wish to set more running at once).

Add ENVI format header
Adding an ENVI header to each of the files (using information from the JAXA header file) allows the files to be easily read using GDAL. Following a similar method outlined in previous posts the script to add ENVI headers can be downloaded from add_gmw_palsar_header.py

It is run using:

python add_gmw_palsar_header.py KC_*

This will add a header for all files in each directory matching ‘KC*’

Create mosaic using RSGISLib

Using the ‘rsgislibmosaic.py’ utilities from RSGISLib, a mosaic can be created for each dataset using:

rsgislibmosaic.py -i . -s '*HH' -o palsar_hh.kea -ot UInt16
rsgislibmosaic.py -i . -s '*HV' -o palsar_hv.kea -ot UInt16
rsgislibmosaic.py -i . -s '*_linci' -o palsar_linci.kea -ot Byte
rsgislibmosaic.py -i . -s '*_mask' -o palsar_mask.kea -ot Byte
rsgislibmosaic.py -i . -s '*_date' -o palsar_date.kea -ot UInt16

A similar procedure should work for other datasets by modifying the Python script used to produce the ENVI header files.

[1] O. Tange (2011): GNU Parallel – The Command-Line Power Tool, ;login: The USENIX Magazine, February 2011:42-47.

Tracking changes in a LaTeX document

One of the problems people often have with using LaTeX for collaborative writing is that it is difficult to track changes in a document, like in Word.

As .tex files are text documents version control such as Git or Mercurial can be used to keep track of changes in the LaTeX source. However, looking at differences in .tex files is not as easy as having a formatted copy of the document with the changes marked, in particular for people not used to LaTeX.

The latexdiff command can be used to create a LaTeX document with the changes marked, from which a PDF can be created showing where the text has changed. For example:

latexdiff_output

Basic usage is:

latexdiff origionaldoc.tex changeddoc.tex > changes.tex

By default the command will print everything to the terminal so the output needs to be redirected (using >) to a file.

When writing papers we commonly have a shell script to generate a change tex file, create a PDF and remove temp files.

# Create diff file
latexdiff --exclude-textcmd "section,subsection,sub subsection" \
   201403_GOBIA_RSGISLib_RIOS_origional.tex \
   201403_GOBIA_RSGISLib_RIOS.tex > \
   201403_GOBIA_RSGISLib_RIOS_changes_temp.tex

# Create PDF
pdflatex 201403_GOBIA_RSGISLib_RIOS_changes_temp.tex
bibtex 201403_GOBIA_RSGISLib_RIOS_changes_temp
pdflatex 201403_GOBIA_RSGISLib_RIOS_changes_temp.tex
pdflatex 201403_GOBIA_RSGISLib_RIOS_changes_temp.tex

# Rename PDF and move extra files.
mv 201403_GOBIA_RSGISLib_RIOS_changes_temp.pdf 201403_GOBIA_RSGISLib_RIOS_changes.pdf
rm 201403_GOBIA_RSGISLib_RIOS_changes_temp.*

Note sometimes changes for particular tex commands can cause problems when creating a PDF, producing an error along these lines:

! Argument of \UL@word has an extra }.
    <inserted text>
                \par 
l.292 ...\DIFaddbegin \DIFadd{Methods}\DIFaddend }
                                                  \label{sec:method}

To stop this happening you can exclude particular commands (e.g., section) using the ‘–exclude-textcmd’ flag. (Thanks to this post on stack exchange for the tip).

Reading KEA files in MATLAB

The KEA file format is built on HDF5 so anything which can read HDF5 can read a KEA file, including MATLAB. I’m not quite sure why you’d want to use MATLAB instead of using Python, in particular as you can read KEA files and many other formats through GDAL. However, I’ve been claiming it’s possible to use KEA with MATLAB for a while so I though it would be a good idea to post an example!

This example uses CASI data over Injune, which is part of the test dataset for RSGISLib, you can download from here.

1. Get file information

info = h5info('injune_p142_casi_sub_ll.kea');
% Get bandnames
info.Groups.Name

2. Read in spatial information
If you don’t care about the geospatial information you can skip this step.

% Read in resolution
res = h5read('injune_p142_casi_sub_ll.kea','/HEADER/RES');
xRes = res(1);
yRes = res(2); 

% Top left coordinate
topLeft = h5read('injune_p142_casi_sub_ll.kea','/HEADER/TL');
topLeftX = topLeft(1); % Longitude
topLeftY = topLeft(2); % Latitude

% Size
imageSize = h5read('injune_p142_casi_sub_ll.kea','/HEADER/SIZE');
imageSizeX = imageSize(1);
imageSizeY = imageSize(2);

% Calculate bottom right
bottomRightX = topLeftX + (double(imageSizeX) * xRes);
bottomRightY = topLeftY + (double(imageSizeY) * yRes);

lon_axis = linspace(topLeftX,bottomRightX,imageSizeX);
lat_axis = linspace(bottomRightY,topLeftY,imageSizeY);

3. Read in data
This assumes reading in bands 12, 8 and 3.

band12 = transpose(h5read('injune_p142_casi_sub_ll.kea','/BAND12/DATA'));
band8 = transpose(h5read('injune_p142_casi_sub_ll.kea','/BAND8/DATA'));
band3 = transpose(h5read('injune_p142_casi_sub_ll.kea','/BAND3/DATA'));

4. Display data
A three band composite is created, this is stretched and then displayed.

% Create composite
composite = cat(3,band12,band8,band3);

% Apply linear stretch
stretchedComposite = imadjust(composite,stretchlim(composite));

% View image
h = imagesc(lon_axis, lat_axis, stretchedComposite);

If you’ve used the image from the example the output should look something like this:

casi_data_in_matlab

More details about the KEA format (including a description of the data structure) are available in the following paper: Peter Bunting and Sam Gillingham, The KEA image file format, Computers & Geosciences, Volume 57, August 2013, Pages 54-58, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2013.03.025.

Create an animated GIF from a series of images

ImageMagick is a collection of command line tools which allows you to do lots of cool things with images. I use it a lot to create animations, in particular to add multiple images and text to a single slide. If you want to make something simpler than a movie, it is possible to generate a gif animation using ImageMagick. The simplest command is:

convert -delay 100 *files.png animation.gif

Which will create an animation from all files matching ‘*files.png’ with a 100 x 1/100 s (1 s) delay between images.

If you want to add a title for each image you can also do this with ImageMagick, using something similar to:

convert image.png -gravity south \
 -strokewidth 10 -pointsize 100 \
 -fill white -annotate 0 "Title" \
 image_with_title.png 

An example script, using RSGISLib to colour images and imagemagick to create a GIF, is available on github.

Nested proportional symbols in OpenLayers

For the visualising data from the SoilSCAPE in-situ network of soil moisture sensors I wanted to use proportional symbols, to provide a visual comparison of the soil moisture from different nodes. The data is received in real time, so I wanted the maps to automatically update and show the current soil moisture. Additionally, as there were three sensors attached to each node I wanted to use three circles, nested within each other, with a different colour for each circle, to show the current moisture for all three sensors, which made it slightly more tricky.

This is the final map for Tonzi Ranch.

To produce the maps three main components were used: OpenLayers, imagemagick, and a Python script running as a cron job.

Maps – OpenLayers

All the maps on the SoilSCAPE site use OpenLayers, which allows you to display a text file using:


    var gsat = new OpenLayers.Layer.Google(
        "Google Satellite",
        {type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22});

   // Read soil moisture sensor information (location, icon, description etc., from text file)
   var layer = new OpenLayers.Layer.Vector("Soil moisture sensors", {
        strategies: [new OpenLayers.Strategy.BBOX({resFactor: 1.1})],
        protocol: new OpenLayers.Protocol.HTTP({
    url: './all/all-points_sm.txt',
    format: new OpenLayers.Format.Text()
        })
    });

    // Add controls to map
    map = new OpenLayers.Map( 'map');

    // Add layers to map
    map.addLayers([gsat, layer]);
    map.addControl(new OpenLayers.Control.Zoom());
    map.addControl(new OpenLayers.Control.LayerSwitcher());
    map.zoomToExtent(mapBounds.transform(proj, map.getProjectionObject()));

The format for the text file is:

point	title	description	icon	iconSize	iconOffset
36.0021,-98.63108333	Node#101	Descriptions	icon.png	27	-2.5,-2.5

Real time updates – Python

The measurements from the sensors are sent to our lab via SMS and added to a MySQL database. The open layers text file is generated at regular intervals using a Python script to pull data from the database running as a cronjob. This approach is used to create maps of sensor status, with the link to the appropriate colour symbol saved in the open layers file based on the last received measurement. A description is also saved in the file with a link to a plot of measurements (also generated using a Python script running as a cronjob).

Symbols – imagemagick

To create the nested symbols imagemagick is used to create a different symbol for each node, by adjusting the size of three symbols and stacking them together. A Python script is used for this which can be downloaded from SoilSCAPE Utilities on Bitbucket. The script is run using:

python createPropSymbol.py outDIR outSymbolName.png Val1 Val2 Val3

To embed the maps within our Drupal site, an iframe is used.

More information about the SoilSCAPE project is available on the project website. Data can also be downloaded using the links at the bottom of each map. For more details on the wireless sensor aspect of the project see some of Agnelo Silva’s publications. The project is funded by NASA ESTO.