[geotk] convert/reproject netcfd to png

[TazMainiac]

I have a netcdf file with GTG data in it that I need to warp and convert to
a PNG for EPSG:4326 and EPGS:3857.

Information from the file:

Variables:
    var Clear_air_turbulence_CAT_altitude_above_msl
        attrs:
            long_name = Clear air turbulence (CAT) @ Specific altitude
above mean sea level
            units = %
            abbreviation = CAT
            grid_mapping = LambertConformal_Projection
            Grib_Variable_Id = VAR_0-19-22_L102
            Grib2_Parameter = null
            Grib2_Parameter_Discipline = Meteorological products
            Grib2_Parameter_Category = Physical atmospheric Properties
            Grib2_Parameter_Name = Clear air turbulence (CAT)
            Grib2_Level_Type = null
            Grib2_Generating_Process_Type = Forecast
            _ChunkSize = null
            missing_value = null
            _FillValue = null
        Dimensions:
            time - length: 6
            altitude_above_msl - length: 36
            y - length: 337
            x - length: 451

I'm able to use NetcdfImageReader to extract a Raster from the file for
each elevation and time using SpatialImageReadParam:

        NetcdfImageReader aReader = new NetcdfImageReader(null);
        aReader.setInput(aNetcdfFileName);

        SpatialImageReadParam aSpatialImageReadParams= aReader
.getDefaultReadParam();

        DimensionSlice timeSlice =
aSpatialImageReadParams.newDimensionSlice();
        timeSlice.addDimensionId("time");

        DimensionSlice depthSlice =
aSpatialImageReadParams.newDimensionSlice();
        depthSlice.addDimensionId( "altitude_above_msl" );

        GridEnvelope aGridEnv = aReader.getGridEnvelope(anImageIndex);
        int numberOfElevationSlices = aGridEnv.getHigh(2) -
aGridEnv.getLow(2) + 1;
        int numberOfTimeSlices      = aGridEnv.getHigh(3) -
aGridEnv.getLow(3) + 1;

        for (int time = 0; time < numberOfTimeSlices; time++)
        {
            timeSlice.setSliceIndex(time);
            for (int elevation = 0; elevation < numberOfElevationSlices;
elevation++)
            {
                depthSlice.setSliceIndex(elevation);
                Raster raster = aReader.readRaster(anImageIndex,
aSpatialImageReadParams);
                // do stuff with raster
            }
        }

I use ImageCoverageReader to get the original GridCoverage2D, which gives
me the src CRS.

My target CRS(s) were generated like:

EGPS_4326_CRS = org.geotoolkit.referencing.CRS.decode("EPSG:4326", true);
EGPS_3857_CRS = org.geotoolkit.referencing.CRS.decode("EPSG:3857", true);

Then I use CRS.findMathTransform to get:

Concat_MT[Param_MT["Affine",
    Parameter["num_row", 3],
    Parameter["num_col", 5],
    Parameter["elt_2_2", 0.0],
    Parameter["elt_2_4", 1.0]],
  Param_MT["LambertConformal",
    Parameter["grid_mapping_name", "lambert_conformal_conic"],
    Parameter["latitude_of_projection_origin", 25.0],
    Parameter["longitude_of_central_meridian", 265.0],
    Parameter["standard_parallel", 25.0],
    Parameter["earth_radius", 6371229.0]],
  Param_MT["Geographic/geocentric conversions",
    Parameter["dim", 2],
    Parameter["semi_major", 6371007.0, Unit["metre", 1]],
    Parameter["semi_minor", 6371007.0, Unit["metre", 1]]],
  PARAM_MT["Geographic/geocentric conversions",
    Parameter["dim", 2],
    Parameter["semi_major", 6378137.0, Unit["metre", 1]],
    Parameter["semi_minor", 6356752.314245179, Unit["metre", 1]]]]

I then use GridCoverageBuilder to take the wrap each of original
time/elevation Rasters and get a new GridCoverage2D:

          GridCoverageBuilder builder = new GridCoverageBuilder();
          builder.setExtent(originalRaster.getWidth(),
originalRaster.getHeight());
          builder.setPixelAnchor(PixelInCell.CELL_CENTER);
          builder.setCoordinateReferenceSystem(srcCrs);
          builder.setRenderedImage((WritableRaster)originalRaster);
          builder.setGridToCRS(mathTransfGridToCrs);
          GridCoverage2D translatedGridCoverage =
builder.getGridCoverage2D();

Then I use Operations.DEFAULT.resample to get a resampled GridCoverage2D
which I extract the Raster from again:

                    GridCoverage2D targetGridCov = (GridCoverage2D)
Operations.DEFAULT.resample(translatedGridCoverage, targetCrs);
          Raster targetRaster = targetGridCov.getRenderedImage().getData();

Finally I convert that Raster using a color map and generate the PNG and
envelope information.

All of this seems extremely convoluted and I'm pretty sure it's not working
the way I expect.  The envelope data for the 4326 conversion seems OK:

    Original Envelope: POINT(-3338.927816840278
-595.6653045472644)/POINT(2769.906332465278 3969.028890973045)
        4326 Envelope: POINT(-126.18477714798462
16.312191763953855)/POINT(-57.26792364467286 55.696065525552626)

But the 3857 envelope is bogus:

        3857 Envelope: POINT(-1.404682513797635E7
1840904.7739340013)/POINT(-6375036.098913054
7498147.296812098)

As I said, this seems like it's overly complex and I'm probably doing it
wrong.

Additionally - is there a way I access the elevation and time index values
using geotoolkit?  I'm currently using ucar GridDataSet to get the
dateTimeAxis and vertAxis values.  I don't see a way to get them with
geotoolkit?

Thanks
TazMainiac
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/geotoolkit/attachments/20160621/a849a4a1/attachment.html>