[geotk] Building a GeographicCRS from scratch

[Martin Desruisseaux]

Hello Pablo

Le 01/02/12 13:15, Pablo Rozas Larraondo a écrit :
> I'm working with a Grib file that uses a Lambert Conformal Conic projection. 
> This projection doesn't correspond to any registered EPSG code but all the 
> information about the projection is in the descriptor of the grib file using 
> the standard WMO format 
> (http://dss.ucar.edu/docs/formats/grib/gribdoc/lcgrid.html).

I was not aware of this documentation, thanks for the URL.

> To solve this issue I'm thinking in using the Geotoolkit API and creating a 
> GeographicCRS and Conversion objects from scratch to project my file in 
> cartesian coordinates. There are a lot of attributes in the Grib descriptor 
> as: uvRelativeToGrid, LoVinDegrees, earthIsOblate, etc. that I somehow have to 
> translate into parameters of this objects, but I not sure about how to do it... 

To creates a GeographicCRS, we needs first the following objects in that order. 
For each case, you can use either a factory (DatumFactory, CSFactory, 
CRSFactory) or instantiate directly the object (new DefaultEllipsoid(...), new 
DefaultGeographicCRS(...), /etc/.). Direct instantiation is easier but factories 
gives you more portability; you could run your code using the Proj.4 bindings 
and compare the result for example. Following code uses direct instantiation for 
simplicity.

     import javax.measure.unit.SI;
     import org.geotoolkit.referencing.cs.*;
     import org.geotoolkit.referencing.crs.*;
     import org.geotoolkit.referencing.datum.*;
     import org.geotoolkit.referencing.operation.*;

     Ellipsoid ellipsoid = DefaultEllipsoid.createFlattenedSphere("any name of your choice", semiMajorAxisLength, inverseFlattening, SI.METRE);
     // Note: you can also use DefaultEllipsoid.createEllipsoid(...) if you known the semi-minor axis length rather than the inverse flattening.
     // If only something else than "semi-minor axis length" or "inverse flattening" is known, we can probably convert that parameter to one of
     // the above ones. For example if the parameter is "excentricity", there is some formulas for doing the conversion.

     GeodeticDatum datum = new DefaultGeodeticDatum("any name of your choice", ellipsoid);
     // Note: it is also possible to specify the prime meridian if it is something else than Greenwich.

     GeographicCRS baseCRS = new DefaultGeographicCRS("any name of your choice", datum, DefaultEllipsoidalCS.GEODETIC_2D);

     // If the units of the projected CRS is something else than metres, we need to specify that fact in a custom Cartesian CS.
     // If only the units are different and everything else (axes names and directions) are the same than the default CS, we can
     // use the 'usingUnit' convenience method.
     CartesianCS projectedCS = DefaultCartesianCS.PROJECTED.usingUnit(SI.KILOMETRE);

     // Conversion from geographic CRS to projected CRS assuming the projection parameters are known (this is an other topic).
     ParameterValueGroup parameters = ...;
     Conversion conversion = new DefiningConversion("Any name of your choice", parameters);

     // In the special case of projected CRS, we are better to use the factory even if direct instantiation is possible.
     // This is because the factory method will do more work (check for axis order, etc.).
     CRSFactory crsFactory = FactoryFinder.getCRSFactory(null);
     Map<String,?>  properties = Collections.singletonMap(ProjectedCRS.NAME_KEY, "Put here whatever name you want");
     ProjectedCRS crs = crsFactory.createProjectedCRS(properties, baseCRS, projectedCS);



Part of the hard work is to map the NetCDF projection parameters with the OGC or 
EPSG parameter names expected by Geotk. I attached to this email an initial 
mapping that I found. This is still work in progress and should be committed in 
Geotk soon, but in the main time this map may be used for making progress with 
the current Geotk version.



> Another issue I've found working with my file is the difference between the 
> NetCDF and the GeoToolkit API when they read the extreme values of the x&y 
> axis values. The difference is small but significant, is this normal?
>
> CODE NETCDF API:
> BOX2D(-704.9721645900614, -584.9804882997278, 705.0278354099386, 
> 585.0195117002722)
> (...snip...)
> CODE GEOTOOLS API:
> BOX4D(-706.2221645900614 -586.2304882997278 NaN NaN, 706.2778354099386 
> 586.2695117002722 NaN NaN)

The NetCDF API maps pixel center, will the Geotk API in this case maps the fully 
enclosing bounding box (from upper-left corner to lower-right corner). If 
everything worked well, the difference between the above BOX2D should be half a 
pixel width/height.

     Regards,

         Martin

-------------- next part --------------
#
# Mapping from NetCDF names to EPSG names.
# Projection names are unindented.
# Parameter names are indented.
#

AlbersEqualArea                     = OGC:Albers_Conic_Equal_Area,      EPSG:Albers Equal Area
  longitude_of_central_meridian     = OGC:central_meridian,             EPSG:Longitude of false origin
  latitude_of_projection_origin     = OGC:latitude_of_origin,           EPSG:Latitude of false origin
  standard_parallel                 = OGC:standard_parallel_1,          EPSG:Latitude of 1st standard parallel
  earth_radius

FlatEarth
  longitude_of_projection_origin
  latitude_of_projection_origin
  rotationAngle

LambertAzimuthalEqualArea           = OGC:Lambert_Azimuthal_Equal_Area, EPSG:Lambert Azimuthal Equal Area
  longitude_of_projection_origin    = OGC:longitude_of_center,          EPSG:Longitude of natural origin
  latitude_of_projection_origin     = OGC:latitude_of_center,           EPSG:Latitude of natural origin
  earth_radius

LambertConformal                    = OGC:Lambert_Conformal_Conic_1SP,  EPSG:Lambert Conic Conformal (1SP)
  longitude_of_central_meridian     = OGC:central_meridian,             EPSG:Longitude of natural origin
  latitude_of_projection_origin     = OGC:latitude_of_origin,           EPSG:Latitude of natural origin
  standard_parallel                 = OGC:standard_parallel_1,          EPSG:Latitude of 1st standard parallel
  earth_radius

latitude_longitude

Mercator                            = OGC:Mercator_2SP,                 EPSG:Mercator (variant B)
  longitude_of_projection_origin    = OGC:central_meridian,             EPSG:Longitude of natural origin
  latitude_of_projection_origin     = OGC:latitude_of_origin,           EPSG:Latitude of natural origin
  standard_parallel                 = OGC:standard_parallel_1,          EPSG:Latitude of 1st standard parallel

Orthographic                        = OGC:Orthographic,                 EPSG:Orthographic
  longitude_of_projection_origin    = OGC:central_meridian,             EPSG:Longitude of natural origin
  latitude_of_projection_origin     = OGC:latitude_of_origin,           EPSG:Latitude of natural origin

RotatedLatLon
  grid_south_pole_latitude
  grid_south_pole_longitude
  grid_south_pole_angle

RotatedPole
  grid_north_pole_latitude
  grid_north_pole_longitude

Stereographic                       = OGC:Stereographic,                EPSG:Stereographic
  longitude_of_projection_origin    = OGC:central_meridian,             EPSG:Longitude of natural origin
  latitude_of_projection_origin     = OGC:latitude_of_origin,           EPSG:Latitude of natural origin
  scale_factor_at_projection_origin = OGC:scale_factor,                 EPSG:Scale factor at natural origin

TransverseMercator                  = OGC:Transverse_Mercator,          EPSG:Transverse Mercator
  longitude_of_central_meridian     = OGC:central_meridian,             EPSG:Longitude of natural origin
  latitude_of_projection_origin     = OGC:latitude_of_origin,           EPSG:Latitude of natural origin
  scale_factor_at_central_meridian  = OGC:scale_factor,                 EPSG:Scale factor at natural origin

UtmProjection
  semi-major_axis
  inverse_flattening
  UTM_zone
  north_hemisphere

VerticalPerspectiveView
  longitude_of_projection_origin
  latitude_of_projection_origin
  height_above_earth
  false_easting
  false_northing
  earth_radius