[geotk] Grib file coordinate transformation

Martin Desruisseaux martin.desruisseaux at geomatys.fr
Mon Jan 30 14:51:55 EST 2012


Hello Pablo

Ideally, the NetCDF reader would have setup the ProjectedCRS automatically. In 
practice, this work is in progress right now so a temporary workaround is 
needed. Instead of:

> GridCoverage coverage = CoverageIO.read("~/myfile.grb");

we will need something more elaborated:

     NetcdfImageReader netcdfReader = new NetcdfImageReader(null);
     netcdfReader.setInput("~/myfile.grb");

     ImageCoverageReader reader = new ImageCoverageReader();
     reader.setInput(netcdfReader);

     GridCoverage2D coverage = reader.read(0, null);
     reader.dispose();
     netcdfReader.dispose();


The above code should produce the same GridCoverage2D than the one produced by 
the convenience CoverageIO method. Once we verified that it work, you can apply 
the workaround described in the "/How do I declare the actual georeferencing of 
an image file?/" FAQ on 
http://www.geotoolkit.org/modules/coverage/faq.html#specifyCRS (you can ignore 
the discussion about the widget - unless you want to see your metadata - and 
ignore all the code example below "/// If the 'grid to CRS' transform need to be 
modified.../").

At this point, you are almost done. Replace the current code:

> values = gc2D.evaluate(new Point2D.Float(600, -420), values);

by:

   DirectPosition pos = new DirectPosition2D(DefaultGeographicCRS.WGS84, longitude, latitude);
   values = gc2D.evaluate(pos, values);


However before doing so, it may be worth to check if your ProjectedCRS is the 
expected one. The magnitude of the (600, -420) coordinates make me suspect that 
they are kilometres rather than metres, in which case we need to instruct the 
CRS so. First, I suggest that you try the following:

   MathTransform tr = CRS.findMathTransform(crs, DefaultGeographicCRS.WGS84);
   System.out.println(tr.transform(pos));


where 'pos' is the above DirectPosition object. See if the result seems 
reasonable. If not because the linear units should be kilometres, we would need 
to change the construction of the CartesianCS object.

     Regards,

         Martin



More information about the Geotoolkit mailing list