Serving WCS elevation data to World Wind using GeoServer

I wrote a post earlier about how to configure GeoServer to serve elevation data to World Wind as a Web Map Service (WMS). Serving elevations over WMS is possible, but requires configuration on both the server and client. World Wind recently introduced a feature to consume elevation data from a Web Coverage Service (WCS). This is a much easier to set up, and requires only server side configuration. WCS support was adding in June 2014, is not yet in the stable World Wind release, so you’ll need a recent snapshot of World Wind.

There is one caveat: World Wind requires that the WCS server return data in a lat-lon coordinate reference system (specifically EPSG:4326). In theory, GeoServer can reproject nearly any source data into this coordinate system, but in practice I’ve had trouble getting the reprojected layer to work with World Wind. You can work around this problem by using GDAL to reproject your source data to the right projection. This will also improve performance, since the server won’t have to reproject data on-the-fly.

Use the gdalwarp command line tool to reproject a file. For example, to reproject the Spearfish elevation GeoTIFF that ships with GeoServer, run this command:

gdalwarp -t_srs EPSG:4326 sfdem.tif sfdem_geographic.tif

Now you can add the warped GeoTIFF as a GeoServer data store, publish the layer, and connect using World Wind. You shouldn’t need any further configuration on the server. See the WCSElevations example in World Wind for how to connect to the WCS server.

Posted in Uncategorized | Leave a comment

Typeset SciDB queries in LaTeX

I often use the listings package in LaTeX for formatting source code in documents. I recently needed to include some SciDB queries in a report that I wrote in LaTeX, so a put together a listing language definition file for the SciDB query languages. The style supports both the Array Query Language (AQL) and the Array Functional Language (AFL). The code is available on github.

Posted in Databases | Tagged | Leave a comment

SciDB tips

I’m working on a project that will use SciDB for an analysis of MODIS data. This post is a dumping ground for small things that I wish I had known at the start of this process. Most of these are obvious in retrospect, but I was confused and maybe others are as well.

Boolean logic operators

I couldn’t find any reference to Boolean logic operators in the SciDB documentation. It turns out that these are the same as in SQL: and, or, not. For example:

SELECT FROM my_array WHERE a > 0 AND (b = 1 OR c = 2);

List of all SciDB functions

The SciDB user manual breaks functions into categories, but sometimes it’s more useful to see a list of all the available functions. This is easy to get using the list operator:

AFL% list('functions');

Use list('operators'); to see available operators.

Use null with the between operator

The between operator selects a chunk of an array by restricting one or more dimensions to a range. But you have to specify a high and low value for each dimension. I started out by copying the min and max possible values for dimensions that I did not want to filter. I later learned that you can pass null for dimensions that you don’t want to filter (the manual doesn’t seem to mention this). For example, I have an array of vegetation index data from MODIS MOD13Q1 product. The schema of my array is:

mod13q1 <evi:int16> [year  = 2001:2001,  1,      0,
                     month = 1:12,       12,     0,
                     day   = 1:31,       10,     0,
                     i     = 0:23039999, 170000, 0]

This array have four dimensions: year, month, day, and i (the pixel index). If I wanted to extract measurements for January to March in all years I could use this query:

between(mod13q1, null, 1, null, null, null, 3, null, null)

Set the missing reason code

SciDB supports multiple types of missing data. Each null cell can have a missing reason code that indicates why the data is missing (for example, instrument failure, cloud cover, etc.) The SciDB manual explains how to represent the missing data code when loading data into SciDB, but it wasn’t clear how to set the missing data code once data was in a SciDB array. I had large array already loaded into SciDB, and I wanted to replace a certain value with a null value and a missing reason code. Turns out that the function missing can do this. In my case, I wanted to replace the value -1000 with the missing value code 1, which I accomplished with the following query:

UPDATE gimms SET ndvi = missing(1) WHERE ndvi = -1000;
Posted in Data analysis, Databases | Tagged | Leave a comment

How to call an R function from C

Last year I worked on an analysis project in which we needed to call a statistical function written in R from a C program. There are some great tools to integrate R and C++ (for example Rccp and Rinside), but for this particular project we wanted a plain C solution with minimal dependencies. Here I’ll describe the basics of how to call R from plain C and provide a minimal working example (full code available on github). For more information Writing R Extensions is a wealth of information, but it’s not easy reading for a newcomer to R internals. There are also some great stackoverflow threads on this topic: Calling R Function from C++, and R from C — Simplest Possible Helloworld.

Suppose that we want to create an array in C, pass that array to a function in R, and receive the result of R’s computation. For the sake of illustration we will call a function that simply adds one to each element in the array:

add1 <- function(a) {
  cat("R received: ", a, "\n");

  return(a + 1)
}

To call this function from C we need to perform a few steps:

Set up an embedded R environment

#include <Rinternals.h>
#include <Rembedded.h>

// Intialize the embedded R environment.
int r_argc = 2;
char *r_argv[] = { "R", "--silent" };
Rf_initEmbeddedR(r_argc, r_argv);

Load the R function definition.

We can do this by writing a C version of the R source function:

/**
* Invokes the command source("foo.R").
*/
void source(const char *name)
{
    SEXP e;

    PROTECT(e = lang2(install("source"), mkString(name)));
    R_tryEval(e, R_GlobalEnv, NULL);
    UNPROTECT(1);
}

This snippet introduces several important concepts. R internals use SEXP (S-expressions) to represent function calls and arguments. In this example, the S-expression is basically a pair in which the first element is the function to invoke and the second element is the arguments to the function. The above code builds an SEXP to call the function source with the argument name, and then evaluates that expression in the R environment.

Also notice the use of the PROTECT macro. This macro tells R than an object will be used by C code, so it must not be garbage collected. PROTECT adds the object to a protection stack, so we are done with the object we must call UNPROTECT(1), which pops one object off of the stack. (See Writing R Extensions 5.9.1).

Call the function

Suppose that we want to invoke our R function on the following C array:

int a[] = { 1, 2, 3, 4, 5 };
int alen = 5;

First we need to create an R vector and copy the data from the C array to the vector:

// Load the R function (source function defined above)
source("func.R")

// Allocate an R vector and copy the C array into it.
SEXP arg;
PROTECT(arg = allocVector(INTSXP, alen));
memcpy(INTEGER(arg), a, alen * sizeof(int));

Note the use of the INTEGER macro, which gives us an int* view of the SEXPR. There are corresponding macros for other data types (e.g., REAL).

Next we need to set up and evaluate an SEXP for the add1 function call:

// Setup a call to the R function
SEXP add1_call;
PROTECT(add1_call = lang2(install("add1"), arg));

// Execute the function
int errorOccurred;
SEXP ret = R_tryEval(add1_call, R_GlobalEnv, &errorOccurred);

And finally we can get a C pointer to the function’s result and do something with the results in our C code:

if (!errorOccurred)
{
    double *val = REAL(ret);

    printf("R returned: ");
    for (int i = 0; i < LENGTH(ret); i++)
        printf("%0.1f, ", val[i]);
    printf("\n");
}

Don’t forget to clean up protected memory:

// Unprotect add1_call and arg 
UNPROTECT(2);

// Release R environment
Rf_endEmbeddedR(0);

Compile and run the example

Note that the R_HOME environment variable must be defined when the program is run.

export R_HOME=/Library/Frameworks/R.framework/Resources
export LD_LIBRARY_PATH=$R_HOME/lib

cc -o r_test -g -I$R_HOME/include -L$R_HOME/lib -lR -lRblas r_test.c

./r_test

And the output is:

R received:  1 2 3 4 5
R returned: 2.0, 3.0, 4.0, 5.0, 6.0,

The full example is available on github.

Posted in Programming | Tagged | Leave a comment

Serving BIL Elevation Data with Geoserver

I’ve recently been looking into how to use the GeoServer map server to serve elevation data to the NASA World Wind virtual globe. It is possible to do this, but there are some complications. In this article I’m going to explain some of the issues.

Update: World Wind recently added support for consuming elevation data from a Web Coverage Service (WCS), which gets around many of the challenges of working with BIL. For many applications this may be a better solution than using WMS and BIL. See my post on working with WCS elevation data in GeoServer for more details and some caveats.

Update: I recently submitted a GeoServer patch that adds a configuration panel to the BIL plugin. This allows the server administrator to configure GeoServer to work with World Wind without any client side configuration. See DDS/BIL Plugin Documentation for details.

World Wind retrieves elevation data on demand from a Web Map Service (WMS), and expects this data to be available in the Band Interleaved Line (BIL) format. GeoServer, out the box, does not support output in BIL, but there is a plugin that adds support for this format. So the steps to serve elevation to World Wind are

  1. Install the GeoServer BIL plugin.
  2. Configure elevation raster data as a layer in GeoServer.
  3. Configure an elevation model in World Wind to connect to the server.

Let’s look at an example of how to do this, and then discuss some of the pitfalls. GeoServer comes with a sample data set of elevation data for Spearfish, South Dakota. This data is provided as a GeoTIFF, and we will use the BIL plugin to convert it on-the-fly to BIL. First we need to configure the server to serve this data as BIL, and then connect to the WMS layer using World Wind.

  1. Download GeoServer and unzip the archive. I’ll refer to the directory of the unzipped GeoServer install as GEOSERVER_ROOT
  2. Download the DDS/BIL plugin and extract the archive into GEOSERVER/webapps/geoserver/WEB-INF/lib.
  3. Start the server by running GEOSERVER/bin/startup.sh.
  4. Configure an elevation model in World Wind to use the layer served by GeoServer. Here’s one way to do this:
// Parse the capabilities document from the WMS
final String GET_CAPABILITIES_URL = "http://localhost:8080/geoserver/wms?request=getCapabilities";
WMSCapabilities caps = WMSCapabilities.retrieve(new URI(GET_CAPABILITIES_URL));
caps.parse();

// Configure parameters for the Spearfish elevation model.
AVList params = new AVListImpl();
params.setValue(AVKey.LAYER_NAMES, "sf:sfdem");
params.setValue(AVKey.IMAGE_FORMAT, "application/bil32");
params.setValue(AVKey.BYTE_ORDER, AVKey.BIG_ENDIAN);
params.setValue(AVKey.MISSING_DATA_SIGNAL, -9.999999933815813E36);

Factory factory = (Factory) WorldWind.createConfigurationComponent(AVKey.ELEVATION_MODEL_FACTORY);
final ElevationModel spearfish = (ElevationModel) factory.createFromConfigSource(caps, params);

SwingUtilities.invokeLater(new Runnable()
{
    public void run()
    {
        // Get the WorldWindow's current elevation model.
        Globe globe = AppFrame.this.getWwd().getModel().getGlobe();

        // Replace elevation model with imported elevations. This makes the elevation 0 everywhere
        // except in the region imported, so it is easy to tell that the elevations are being pulled
        // from geoserver. For production use create a CompoundElevationModel and add the new elevations
        // to the compound model.
        globe.setElevationModel(spearfish);

        // Set the view to look at the imported elevations.
        Position spearfishSouthDakota = Position.fromDegrees(44.4709, -103.6812, 10000);
        getWwd().getView().setEyePosition(spearfishSouthDakota);
   }
});

(Full example code available on github.)

The example above is sufficient to load elevations from GeoServer. However, you will notice that the example required several parameters to be configured on the World Wind side. Unfortunately, the client needs three pieces of information to interpret the elevation data, and GeoServer does not provide this info. Specifically, we need to know the data type, the byte order, and the missing data value of the BIL files. To understand these settings it is helpful to understand the details of the BIL file format.

What is BIL anyway?

BIL stands for Band Interleaved Line, and is the file format that World Wind uses to store elevation data. More generally, BIL is a simple data format for multi-band raster data. There is a good description of the file format in the ArcGIS User Manual.

BIL files are often accompanied by a text header file that specifies the data type, missing data value, geographic extent, band names, etc. There are (at least) two different formats for these files: the ESRI header format and the ENVI header format. Without a header file a BIL file is just an opaque bunch of bytes and the client must know the data type, array dimensions and so on in order to interpret the file. When GeoServer returns BIL data to World Wind it provides only the raw bytes, not the header file. So the programmer needs to provide the missing information.

Data type

When World Wind asks a WMS server for a map it specifies what type of file the server should return. There are several MIME types used for BIL files:

  • image/bil
  • application/bil
  • application/bil16
  • application/bil32

If you use the application/bil16 or application/bil32 then World Wind can infer the data type from the MIME type (either 16 bit integer or 32 bit floating point). However, World Wind assumes that image/bil and application/bil are 16 bit integer, which is not always the case for data served by GeoServer. To avoid data type mismatches it’s best to explicitly request either application/bil16 or application/bil32.

params.setValue(AVKey.IMAGE_FORMAT, "application/bil32");

Byte order

As well as the data type, World Wind needs to know the byte order (big or little endian) of the BIL files. If not specified, World Wind will assume that the data is little endian. However, the GeoServer BIL plugin produces files in big endian order (network byte order). So the client code must specify this byte order:

params.setValue(AVKey.BYTE_ORDER, AVKey.BIG_ENDIAN);

Missing data signal

Missing data value in GeoServer interface

Refer to the “Coverage Band Details” section in the GeoServer Layer editor to find the missing data value for your elevation layer.

Elevation data sets often include a special value that marks pixels for which no data is available. Data might be unavailable because the data set covers only a certain area, or because the elevation could not be measured for some pixels. World Wind needs to know that it should ignore these values.

The BIL plugin does not set the missing data value, so whatever missing data value is used in the layer source data will be passed to the generated BIL files. In the case of the Spearfish DEM, this value is -9.999999933815813E36. You can see the missing data value in the GeoServer admin interface under the layer editor page. You can also find the missing data value using the gdalinfo command line tool.

By default World Wind uses -9999 as the missing data value. If this is not the value set in your elevation source data you will need to either set the missing data value in your layer parameters or change the missing data value in the source data.

params.setValue(AVKey.MISSING_DATA_SIGNAL, -9.999999933815813E36);

Changing the missing data value using GDAL

It is straight forward to change the missing data value of a raster file using the gdalwarp tool. To translate the Spearfish DEM to use -9999 as the missing data value (the World Wind default):

gdalwarp -dstnodata -9999 sfdem.tif sfdem_out.tif

Summary

It is possible to configure GeoServer to serve elevation data in a format that World Wind can consume, but it requires some configuration on both sides. If elevations are not rendered correctly you have probably misconfigured either the data type, byte order, or missing data settings in client code.

Posted in Geospatial | Tagged , | Leave a comment

First day in Morocco

camels

I arrived in Morocco today, via ferry from Spain. I’m spending the first night at a youth hostel in Tangier. First impressions are that Morocco is intense. So many people try to make a living off of tourists that you constantly have to fend off offers from ‘guides’ and hustlers. Leaving the port I declined at least half a dozen offers to carry my bags, find me a taxi/bus/hotel, etc. On the way to the hostel I had at least as many offers of hashish.

First impressions continued, the Rough Guide to Morocco has yet to prove it’s worth. It listed the wrong bus times, and the city map of Tangier is misleading and got me lost. All that said, Morocco is quite exotic, and very different from Europe. I imagine it will get even more different as I move away from Tangier.

parker-marakesh

Posted in travel | Leave a comment

Lost…er…’Off-route exploration’ in the Sierra Nevada

I just finished a four day hiking trip in the Spanish Sierra Nevada with my Dutch friend Jurgen. I don’t really get Spanish trail maps. They use the same marking for a huge well groomed trail with suspension bridges over each little stream as they do for an overgrown goat track. Or more often the marking seems to be a suggested cross country route, since there’s no trail at all! With some bushwhacking we found our way to a mountain hut just below the peak of Valeta (a big mountain). To get there we undertook a torturous climb up a snowy ski slope, and were absolutely exhausted when we finally reached the hut (even though it was only two in the afternoon). We shared the hut with two Spanish guys and two guys from Israel. The trip ended with an adventurous glissade (technical term for sliding on your ass) down the ski slope.

Posted in travel | Leave a comment

Sunseed

I’ve spent the last month doing volunteer work at Sunseed Desert Technology (www.sunseed.co.uk) in southern Spain (near Almeria). It’s a cool place, with friendly people. The work is quite varied. Cooking, gardening, laying tile, building (or more often repairing) stuff. I’ve enjoyed it. Unfortunately there are no Spanish people here, and my Spanish hasn’t improved a bit for having spent a month in Spain.

This part of Spain is quite beautiful, in a dry desert way. Lots of crumbling, ruined buildings and cactus. A couple weeks back a few of us went to the mountain for the weekend. It was amazing to hike up from the desert into the snow. We even had four inches of snow overnight, as we staying in a mountain shelter.

On Tuesday I move on to…somewhere else. Maybe back to Granada to start with.

Posted in travel | Leave a comment

First week in Spain

I arrived in Spain last Thursday, in Madrid. I hadn’t bothered to book a bed in advance, since it was mid-week during the off season. And of course I arrived in Madrid after a late night and a day of travel to find that all the hostels were booked up on account of some art show. After calling around a bit I found a bed, but only for one night. The morning I had to go through the hunt again to find a hostel with a single bed free, and again only for one night. Two days of that was enough for me and I left Madrid for Valencia, which is a nicer city anyway, and my friend Maria lives here.

Other than the stress of finding accommodation, I’m enjoying Spain. I met some nice people at the hostel in Madrid. The weather is a lot warmer than it was in Ireland too. My Spanish is coming back to me, though it needs a lot of improvement. But I guess I’m in the right country for that. Madrid I wasn’t all that impressed with. Too much of a big city. Valencia I like more. Though I’m surprised that the streets are so wide here, and that there are so many cars. Maria told me that there are so many cars in Valencia that it’s legal to double park. The person parking double has to leave the emergency brake off so that other people can roll the car out of the way.

Posted in travel | Leave a comment

Jampa Ling

I’ve spent the last week working at the Jampa Ling Tibetan Buddhist center in exchange for room and board. There’s a community of about ten people here, and occasional guests. The center is kind of in the middle of Ireland, not far from the border with Northern Ireland. It’s a very pretty place, with a couple lakes nearby and a small forest. I’ve been working in the garden, which has been nice, but also a bit chilly since the weather has been cold. It’s snowing today. The people here are very friendly. There is a Tibetan Lama in residence, and he is a really lovely man.

Right now we are all busy getting ready for a teaching weekend. Tomorrow fifty people will arrive for Rinpoche’s teaching, so all the beds need to be made, and floors cleaned and all that.

I will be here one more week, and then I go to Spain! Where it will be warm!

Posted in travel | Leave a comment