[geotk] Reading netcdf GMT GRD as coverage

Martin Desruisseaux martin.desruisseaux at geomatys.fr
Fri Feb 18 12:16:07 EST 2011


Hello Brian

I just commited a missing code logic in Geotk. In about 5 hours from now, you 
should have an updated geotk-SNAPSHOT with it. Please give a new try to your 
code. You should now see the following header in the ASCII file:

NCOLS      5083
NROWS      4067
XLLCORNER  422960.2775983825
YLLCORNER  4118889.9067777675
CELLSIZE   50.0

The "NODATA_VALUE" attribute is missing because it was not declared in the first 
place (the NetCDF file contains NaN values instead, which is fine). Geotk used 
-9999 by default as documented there:

    * http://www.geotoolkit.org/apidocs/org/geotoolkit/image/io/plugin/AsciiGridWriter.html

However you can have more control on the metadata to be written by overloading 
the following method:

     @Override
     protected void completeImageMetadata(IIOMetadata metadata, GridCoverage coverage) throws IOException {
         System.out.println(metadata);
     }

You should see a tree like below:

geotk-coverageio_3.07
├───RectifiedGridDomain
│   ├───origin="422960.2775983825 4118889.9067777675"
│   ├───CoordinateReferenceSystem
│   │   ├───name="EPSG:WGS 84 / UTM zone 10N"
│   │   ├───type="projected"
│   │   ├───Datum
│   │   │   ├───name="EPSG:World Geodetic System 1984"
│   │   │   ├───type="geodetic"
│   │   │   ├───Ellipsoid
│   │   │   │   ├───name="EPSG:WGS 84"
│   │   │   │   ├───axisUnit="m"
│   │   │   │   ├───semiMajorAxis="6378137.0"
│   │   │   │   └───inverseFlattening="298.257223563"
│   │   │   └───PrimeMeridian
│   │   │       ├───name="EPSG:Greenwich"
│   │   │       ├───greenwichLongitude="0.0"
│   │   │       └───angularUnit="deg"
│   │   ├───CoordinateSystem
│   │   │   ├───name="EPSG:Cartesian 2D CS. Axes: easting, northing (E,N). Orientations: east, north. UoM: m."
│   │   │   ├───type="cartesian"
│   │   │   ├───dimension="2"
│   │   │   └───Axes
│   │   │       ├───CoordinateSystemAxis
│   │   │       │   ├───name="EPSG:Easting"
│   │   │       │   ├───axisAbbrev="E"
│   │   │       │   ├───direction="east"
│   │   │       │   └───unit="m"
│   │   │       └───CoordinateSystemAxis
│   │   │           ├───name="EPSG:Northing"
│   │   │           ├───axisAbbrev="N"
│   │   │           ├───direction="north"
│   │   │           └───unit="m"
│   │   └───Conversion
│   │       ├───name="EPSG:UTM zone 10N"
│   │       ├───method="Transverse Mercator"
│   │       └───Parameters
│   │           ├───ParameterValue
│   │           │   ├───name="central_meridian"
│   │           │   └───value="-123.0"
│   │           ├───ParameterValue
│   │           │   ├───name="scale_factor"
│   │           │   └───value="0.9996"
│   │           └───ParameterValue
│   │               ├───name="false_easting"
│   │               └───value="500000.0"
│   ├───OffsetVectors
│   │   ├───OffsetVector
│   │   │   └───values="50.0 0.0"
│   │   └───OffsetVector
│   │       └───values="0.0 -50.0"
│   └───Limits
│       ├───low="0 0"
│       └───high="5082 4066"
├───SpatialRepresentation
│   ├───numberOfDimensions="2"
│   ├───centerPoint="550035.2775983824 4017214.9067777675"
│   └───pointInPixel="upperLeft"
└───ImageDescription
     └───Dimensions
         └───Dimension
             ├───descriptor="z"
             ├───scaleFactor="1.0"
             └───offset="0.0"


To add an explicit fill value, use the following code:

     @Override
     protected void completeImageMetadata(IIOMetadata metadata, GridCoverage coverage) throws IOException {
         MetadataAccessor accessor = new MetadataAccessor(metadata, SpatialMetadataFormat.FORMAT_NAME, "ImageDescription/Dimensions", "Dimension");
         accessor.selectChild(0); // 0 is the band number (only one in this case)
         accessor.setAttribute("fillSampleValues", -9999);
     }


     Martin


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/geotoolkit/attachments/20110218/f408a497/attachment.html


More information about the Geotoolkit mailing list