<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>