<html>
<head>
<meta content="text/html; charset=utf-8" http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
<div class="moz-cite-prefix">
<p>Hello Taz</p>
<p>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 <tt>javax.imageio</tt> API.
This was an appropriate goal for two-dimensional raster formats.
But trying to pass N-dimensional data through the
two-dimensional <tt>javax.imageio</tt> API has proven to be
much more convolved that I initially though.</p>
<p>UCAR has not done this mistake, since their API is independent
from <tt>javax.imageio</tt>. 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.<br>
</p>
<p>Our NetCDF reader will be massively refactored this summer in
the Apache SIS project. Its API will be different, hopefully
much less convolved.</p>
<p>I will come back tomorrow more precisely on the issues that you
raised.<br>
</p>
<p> Regards,</p>
<p> Martin</p>
<p><br>
</p>
Le 21/06/16 à 23:22, TazMainiac a écrit :<br>
</div>
<blockquote
cite="mid:CA+DQs550eWoyKzfB0Uioym3j58tFje9go1QSJXayRWhmJT-S-A@mail.gmail.com"
type="cite">
<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
moz-do-not-send="true" 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>
</blockquote>
<br>
</body>
</html>