[geotk] convert/reproject netcfd to png

Martin Desruisseaux martin.desruisseaux at geomatys.com
Tue Jun 21 16:35:46 PDT 2016


Hello Taz

You have investigated the steps a lot :-). I will take a closer look at
them tomorrow. A quick note in the meantime: the Geotk NetCDF reader is
indeed more complex and convolved than it should be. One bad design
choice that I did was to try to fit the NetCDF reader into the standard
javax.imageio API. This was an appropriate goal for two-dimensional
raster formats. But trying to pass N-dimensional data through the
two-dimensional javax.imageio API has proven to be much more convolved
that I initially though.

UCAR has not done this mistake, since their API is independent from
javax.imageio. But with the UCAR's API we face another difficulty, which
is that their data model is very different than the OGC/ISO models.
Consequently a lot of complexity happen in the adapters that wrap a
UCAR's CoordinateSystem object into an OGC/ISO CoordinateReferenceSystem
object.

Our NetCDF reader will be massively refactored this summer in the Apache
SIS project. Its API will be different, hopefully much less convolved.

I will come back tomorrow more precisely on the issues that you raised.

    Regards,

        Martin


Le 21/06/16 à 23:22, TazMainiac a écrit :
> 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
> <tel: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/20160622/cf06b693/attachment-0001.html>


More information about the Geotoolkit mailing list