NetCDF: projected x/y dimensions in kilometers results in incorrect output


This issue was first reported as against version 16.4. The issue does exist against 18.2. Version 18.2 did a slight refactor where one of the impacted classes was moved to main out of render. However, the code remain largely unchanged.

I have an ImageMosaic store backed by a Netcdf coverage. The X/Y dimensions are defined in kilometers.

The native CRS is based on a rotated polar_stereographic projection, with units defined as 'km':

900914=PROJCS["polar_stereographic", GEOGCS["unknown", DATUM["unknown", SPHEROID["unknown", 6371229.0, 0.0]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH]], PROJECTION["Polar_Stereographic"], PARAMETER["central_meridian", -135.0], PARAMETER["latitude_of_origin", 90.0], PARAMETER["scale_factor", 0.9330127239227295], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], UNIT["km", 1000.0], AXIS["Easting", EAST], AXIS["Northing", NORTH],AUTHORITY["EPSG","900914"]]

First issue, on gt-main:
The following file:
The validArea for the projection defines x bounds as [-Double.MAX_VALUE,Double.MAX_VALUE]
When this envelope is reprojected to the target CRS in the user's request, it results in a much smaller area than what is expected. I am not sure if this is caused by an overflow due to the MAX_VALUE, or by a bug in the AffineTransform2D method that is executed when the reprojection is performed, since we specify UNIT['km', 1000.0] in the target CRS.

This issue causes the eventual image rendered, or WCS output, to be cut off.

I was able to get around the issue replacing the X bounds in the validArea to [-180, 180] (since WGS84 is specified in this envelope). This seemed to fix the issue I was seeing, and the entire domain of the dataset is returned by Geoserver.

The second issue involves the following if block:

Since a UNIT['km', 1000.0] is specified, the instance of the MathTransform is actually a ConcatenatedTransform, which is not handled in this method. This causes the parseCRS method to return null, and eventually throw a NullPointerException.

My fix was to check for a ConcatenatedTransform right before the if block, and extract the transform1, which represents the reprojection, before continuing:

MathTransform transform = projection.getMathTransform();


  • FIX

  • In the case the projection defines a unit different than 'm', for example:

  • UNIT['km', 1000.0]

  • A 'ConcatenatedTransform' object is built. It wraps the two different Transforms necessary:

  • transform1 = the reprojection (i.e. PolarStereographic)

  • transform2 = AffineTransform (to handle the unit conversion from 'm' to 'km')

  • We must extract the ProjectedCRS from the ConcatenatedTransform object.
    MathTransform transform2 = null;
    if (transform instanceof ConcatenatedTransform)
    // TODO Consider verifying that transform2 is not necessary in this method
    // The AffineTransform can be used to determine the units of the NetCDFCoordinate
    // objects below, instead of hardcoding to 'm'.
    transform2 = ((ConcatenatedTransform) transform).transform2;
    // Get the ProjectedCRS
    transform = ((ConcatenatedTransform) transform).transform1;
    /* END FIX */

if (transform instanceof TransverseMercator) {

A third issue exists concerning the x/y NetCDFCoordinate instances further below being hardcoded to the 'm' units. This causes the x/y variables in the Netcdf output to be labeled with the 'm' units, instead of 'km'. Note that the AffineTransform instance, as shown above, can be used to determine which units should be used, instead of hard-coding 'm'. However, I noticed that 'm' seems to be hardcoded in the baseline in multiple places, so I am not yet sure how much refactoring that would involve.


Running GeoServer 2.12.2 and GeoTools 16.4




Jerry Wilwerding



Affects versions