Parker Abercrombie

Serving WCS elevation data to World Wind using GeoServer

Posted on July 20, 2014 by Parker Abercrombie

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

Posted on July 15, 2014 by Parker Abercrombie

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 scidb | Leave a comment

SciDB tips

Posted on May 17, 2014 by Parker Abercrombie

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 scidb | Leave a comment

How to call an R function from C

Posted on May 2, 2014 by Parker Abercrombie

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 R | Leave a comment

Serving BIL Elevation Data with Geoserver

Posted on April 25, 2014 by Parker Abercrombie

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

Install the GeoServer BIL plugin.
Configure elevation raster data as a layer in GeoServer.
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.

Download GeoServer and unzip the archive. I’ll refer to the directory of the unzipped GeoServer install as GEOSERVER_ROOT
Download the DDS/BIL plugin and extract the archive into GEOSERVER/webapps/geoserver/WEB-INF/lib.
Start the server by running GEOSERVER/bin/startup.sh.
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 Geoserver, World Wind | Leave a comment

First day in Morocco

Posted on May 17, 2007 by Parker Abercrombie

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