<div dir="ltr"><div><div>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.<br><br>Information from the file:<br></div><span style="font-family:monospace,monospace"><br></span><div style="margin-left:40px"><span style="font-family:monospace,monospace">Variables:</span><br><span style="font-family:monospace,monospace">    var Clear_air_turbulence_CAT_altitude_above_msl</span><br><span style="font-family:monospace,monospace">        attrs:</span><br><span style="font-family:monospace,monospace">            long_name = Clear air turbulence (CAT) @ Specific altitude above mean sea level</span><br><span style="font-family:monospace,monospace">            units = %</span><br><span style="font-family:monospace,monospace">            abbreviation = CAT</span><br><span style="font-family:monospace,monospace">            grid_mapping = LambertConformal_Projection</span><br><span style="font-family:monospace,monospace">            Grib_Variable_Id = VAR_0-19-22_L102</span><br><span style="font-family:monospace,monospace">            Grib2_Parameter = null</span><br><span style="font-family:monospace,monospace">            Grib2_Parameter_Discipline = Meteorological products</span><br><span style="font-family:monospace,monospace">            Grib2_Parameter_Category = Physical atmospheric Properties</span><br><span style="font-family:monospace,monospace">            Grib2_Parameter_Name = Clear air turbulence (CAT)</span><br><span style="font-family:monospace,monospace">            Grib2_Level_Type = null</span><br><span style="font-family:monospace,monospace">            Grib2_Generating_Process_Type = Forecast</span><br><span style="font-family:monospace,monospace">            _ChunkSize = null</span><br><span style="font-family:monospace,monospace">            missing_value = null</span><br><span style="font-family:monospace,monospace">            _FillValue = null</span><br><span style="font-family:monospace,monospace">        Dimensions:</span><br><span style="font-family:monospace,monospace">            time - length: 6</span><br><span style="font-family:monospace,monospace">            altitude_above_msl - length: 36</span><br><span style="font-family:monospace,monospace">            y - length: 337</span><br><span style="font-family:monospace,monospace">            x - length: 451</span><br></div><div style="margin-left:40px"><span style="font-family:monospace,monospace"></span></div><br>I'm able to use NetcdfImageReader to extract a Raster from the file for each elevation and time using SpatialImageReadParam:<br><br><span style="font-family:monospace,monospace">        NetcdfImageReader </span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">aReader</span> = new NetcdfImageReader(null);<br>        </span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">aReader</span>.setInput(aNetcdfFileName);<br><br>        SpatialImageReadParam </span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">aSpatialImageReadParams</span>= </span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">aReader</span>.getDefaultReadParam();<br><br>        DimensionSlice timeSlice = aSpatialImageReadParams.newDimensionSlice();<br>        timeSlice.addDimensionId("time");<br><br>        DimensionSlice depthSlice = aSpatialImageReadParams.newDimensionSlice();<br>        depthSlice.addDimensionId( "altitude_above_msl" );<br><br>        GridEnvelope aGridEnv = aReader.getGridEnvelope(anImageIndex);<br>        int numberOfElevationSlices = aGridEnv.getHigh(2) - aGridEnv.getLow(2) + 1;<br>        int numberOfTimeSlices      = aGridEnv.getHigh(3) - aGridEnv.getLow(3) + 1;<br><br>        for (int time = 0; time < numberOfTimeSlices; time++)<br>        {          <br>            timeSlice.setSliceIndex(time);<br>            for (int elevation = 0; elevation < numberOfElevationSlices; elevation++)<br>            {<br>                depthSlice.setSliceIndex(elevation);<br>                Raster raster = aReader.readRaster(anImageIndex, aSpatialImageReadParams);<br></span></div><span style="font-family:monospace,monospace">                // do stuff with raster<br></span><span style="font-family:monospace,monospace">            }<br>        }<br><br></span>I use ImageCoverageReader to get the original GridCoverage2D, which gives me the src CRS.<br><br>My target CRS(s) were generated like:<br><br><div style="margin-left:40px"><span style="font-family:monospace,monospace">EGPS_4326_CRS = org.geotoolkit.referencing.CRS.decode("EPSG:4326", true);</span><br><span style="font-family:monospace,monospace">EGPS_3857_CRS = org.geotoolkit.referencing.CRS.decode("EPSG:3857", true);</span><br></div><br>Then I use CRS.findMathTransform to get:<br><br><div style="margin-left:40px"><span style="font-family:monospace,monospace">Concat_MT[Param_MT["Affine",</span><br><span style="font-family:monospace,monospace">    Parameter["num_row", 3],</span><br><span style="font-family:monospace,monospace">    Parameter["num_col", 5],</span><br><span style="font-family:monospace,monospace">    Parameter["elt_2_2", 0.0],</span><br><span style="font-family:monospace,monospace">    Parameter["elt_2_4", 1.0]],</span><br><span style="font-family:monospace,monospace">  Param_MT["LambertConformal",</span><br><span style="font-family:monospace,monospace">    Parameter["grid_mapping_name", "lambert_conformal_conic"],</span><br><span style="font-family:monospace,monospace">    Parameter["latitude_of_projection_origin", 25.0],</span><br><span style="font-family:monospace,monospace">    Parameter["longitude_of_central_meridian", 265.0],</span><br><span style="font-family:monospace,monospace">    Parameter["standard_parallel", 25.0],</span><br><span style="font-family:monospace,monospace">    Parameter["earth_radius", 6371229.0]],</span><br><span style="font-family:monospace,monospace">  Param_MT["Geographic/geocentric conversions",</span><br><span style="font-family:monospace,monospace">    Parameter["dim", 2],</span><br><span style="font-family:monospace,monospace">    Parameter["semi_major", 6371007.0, Unit["metre", 1]],</span><br><span style="font-family:monospace,monospace">    Parameter["semi_minor", 6371007.0, Unit["metre", 1]]],</span><br><span style="font-family:monospace,monospace">  PARAM_MT["Geographic/geocentric conversions",</span><br><span style="font-family:monospace,monospace">    Parameter["dim", 2],</span><br><span style="font-family:monospace,monospace">    Parameter["semi_major", 6378137.0, Unit["metre", 1]],</span><br><span style="font-family:monospace,monospace">    Parameter["semi_minor", 6356752.314245179, Unit["metre", 1]]]]</span><br></div><div><span style="font-family:monospace,monospace"></span><br>I then use GridCoverageBuilder to take the wrap each of original time/elevation Rasters and get a new GridCoverage2D:<br></div><span style="font-family:monospace,monospace"><br></span><span style="font-family:monospace,monospace">          GridCoverageBuilder builder = new GridCoverageBuilder();</span><br><span style="font-family:monospace,monospace">          builder.setExtent(originalRaster.getWidth(), originalRaster.getHeight());</span><br><span style="font-family:monospace,monospace">          builder.setPixelAnchor(PixelInCell.CELL_CENTER);</span><br><span style="font-family:monospace,monospace">          builder.setCoordinateReferenceSystem(srcCrs);</span><br><span style="font-family:monospace,monospace">          builder.setRenderedImage((WritableRaster)originalRaster);</span><br><span style="font-family:monospace,monospace">          builder.setGridToCRS(mathTransfGridToCrs);</span><br><span style="font-family:monospace,monospace">          GridCoverage2D translatedGridCoverage = builder.getGridCoverage2D();</span><br><br>Then I use Operations.DEFAULT.resample to get a resampled GridCoverage2D which I extract the Raster from again:<br><br>                    <span style="font-family:monospace,monospace">GridCoverage2D targetGridCov = (GridCoverage2D) Operations.DEFAULT.resample(</span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">translatedGridCoverage</span>, targetCrs);<br>          Raster targetRaster = </span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">targetGridCov.getRenderedImage().getData();<br><br></span></span>Finally I convert that Raster using a color map and generate the PNG and envelope information. <br><br>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:<br><br><span style="font-family:monospace,monospace">    Original Envelope: POINT(-3338.927816840278 -595.6653045472644)/POINT(2769.906332465278 3969.028890973045)<br>        4326 Envelope: POINT(-126.18477714798462 16.312191763953855)/POINT(-57.26792364467286 55.696065525552626)<br></span><br>But the 3857 envelope is bogus:<br><br><span style="font-family:monospace,monospace">        3857 Envelope: POINT(-1.404682513797635E7 1840904.<a href="tel:7739340013" value="+17739340013" target="_blank">7739340013</a>)/POINT(-6375036.098913054 7498147.296812098)</span><br><br>As I said, this seems like it's overly complex and I'm probably doing it wrong.<br><br>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?<br><br>Thanks<br>TazMainiac</div>