Scala, Geotoolkit and UTM Zone 10N gridded data (Part 2)

by hohonuuli 2/06/2012 02:16:00 PM 0 comments

Overview

In part 1, I illustrated how to set up a project and open a GRD file as a coverage object that we can use for analysis. In this section I'll briefly show how I extract depth from the grid and use it to calculate altitude.

What we have so far

At this point we have a UTM GridCoverage. For my analysis, this grid covers most of the Monterey Canyon (wikipedia).

The Calculation

I need to calculate the altitude of my observations above the ocean bottom. Altitude is just a function of the difference between the observation depth and the depth of the bottom (i.e. bathymetry)

$\mathrm{altitude}=\mathrm{f(depth,\; bathymetry)}$
$$
Here's a quick overview of how to extract a value from the grid:
1) Create a representation of the position of interest (e.g.  36N, -124W). Like a GridCoverage, a position needs a coordinate reference system or CRS so that the software knows the frame of reference to use for your position.
2). 'evaluate' that position in the GridCoverage2D. Keep in mind that some grids have several bands. For example, an image will have 3 color bands. So evaluate doesn't return a single value, but rather an array of values with each value in the array representing a different band.

IMPORTANT: If you evaluate a position that's outside the bounds of the grid it will throw PointOutSideCoverageException.

I'll need a position object to hold my x and y coordinates as well as depth. Along with the position, I'll need a reference to the coordinate reference system that my coordinates belong to. Scala's case classes are perfect for this. Here's some very simple boiler plate to represent my observation:

import org.geotoolkit.geometry.DirectPosition2D
import org.opengis.referencing.crs.CoordinateReferenceSystem
/** 
 * @param latitude decimal degrees
 * @param longitude decimal degres
 * @param depth meters (0 is surface, + is down!! [oceanographic convention])
 * @param crs The coordinate reference system of the observation.
 */
case class Observation(latitude: Double, 
                       longitude: Double, 
                       depth: Double, 
                       crs: CoordinateReferenceSystem) {

    /** Transform to a GIS position */
    lazy val directPosition = new DirectPosition2D(crs, latitude, longitude)
}

Creating a CRS is a relatively expensive operation. I'll create one and reuse it for all my observations:

val wgs84CRS = CRS.decode("EPSG:4326")

Then I'll create my observations as needed:

val observation = Observation(36.5D, -123.001D, 2345D, wgs84CRS)

Next, I'll need a function that converts my Observation objects into altitude.

import org.geotoolkit.coverage.grid.GridCoverage2D
import org.geotoolkit.referencing.CRS
import org.opengis.geometry.DirectPosition
import org.geotoolkit.geometry.DirectPosition2D

/**
 * Function that takes an Observation and calculates altitude from the 
 * given bathymetry.
 *
 * @param gridCoverage This is the bathymetry grid.
 */
class ToAltitude(gridCoverage: GridCoverage2D) extends (Observation => Double) {

    def apply(v1: Observation): Double = {
        try {
            val position = observation.directPosition
 
            // Extracting a value from my grid
            val depths = Array.ofDim[Double](3)
            gridCoverage.evaluate(position, depths)

            // Bathymetry is - down. Depth is + down
            -depths(0) - v1.depth
        }
        catch {
            case e: Exception => {
                // You should log an exception. For brevity, I'll return
                // NaN to represent a missing value. 
                Double.NaN
            }
        }
    }

}

Now the usage is very simple. First, I'll create an instance of my ToAltitude function:

val toAltitude = new ToAltitude(utmZone10N)

And then I can apply it to my observations:

val altitude = toAltitude(observation)