Hi,
with this mail I'm trying to summarize a few findings of mine
on how to deal with 3D vector data, and in particular with
the 2D+1 subset of it, in Geotools.
With 2D+1 data I mean the usual JTS geometry types, whose
coordinates do have a not null, nor NaN, Z component.
Expanding a bit that might also include the usual measured
geometries, so include the xyz, xym, xyzm cases that both
shapefiles and postgis do handle.
The main issues trying to deal with more dimension are
dealing the the CRS representation, the transformations,
and the envelopes.
The cleanest path seemed to have a real 3D CRS and thus
stick the information about the dimensionality in
the CRS itself.
CRS wise there are two kinds of 3D systems, 3d
geographic, like 4327, and compounds, that do take
a 2d crs and stick a vertical component onto it.
Here are two examples of them:
EPSG:4327, the 3d brother of 4326 (but unusable due to
its DMS unit of measure):
GEOGCS["WGS 84 (geographic 3D)",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84", 6378137.0, 298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
UNIT["DMS", 0.00000484813681109536],
AXIS["Geodetic longitude", EAST],
AXIS["Geodetic latitude", NORTH],
AXIS["Ellipsoidal height", UP],
AUTHORITY["EPSG","4327"]]
EPSG:7401, a compound one:
COMPD_CS["NTF (Paris) / France II + NGF Lallemand",
PROJCS["NTF (Paris) / France II",
GEOGCS["NTF (Paris)",
DATUM["Nouvelle Triangulation Francaise (Paris)",
SPHEROID["Clarke 1880 (IGN)", 6378249.2, 293.4660212936269,
AUTHORITY["EPSG","7011"]],
AUTHORITY["EPSG","6807"]],
PRIMEM["Paris", 2.5969213, AUTHORITY["EPSG","8903"]],
UNIT["grade", 0.015707963267948967],
AXIS["Geodetic latitude", NORTH],
AXIS["Geodetic longitude", EAST],
AUTHORITY["EPSG","4807"]],
PROJECTION["Lambert Conic Conformal (1SP)", AUTHORITY["EPSG","9801"]],
PARAMETER["central_meridian", 0.0],
PARAMETER["latitude_of_origin", 52.0],
PARAMETER["scale_factor", 0.99987742],
PARAMETER["false_easting", 600000.0],
PARAMETER["false_northing", 2200000.0],
UNIT["m", 1.0],
AXIS["Easting", EAST],
AXIS["Northing", NORTH],
AUTHORITY["EPSG","27582"]],
VERT_CS["NGF Lallemand",
VERT_DATUM["Nivellement General de la France - Lallemand", 2005,
AUTHORITY["EPSG","5118"]],
UNIT["m", 1.0],
AXIS["Gravity-related height", UP],
AUTHORITY["EPSG","5719"]],
AUTHORITY["EPSG","7401"]]
So, the basic idea was to use a real 3D CRS and
have everything else learn how to handle it.
Unfortunately I've stumbled into a number of issues:
- all of our code uses ReferencedEnvelope, which is
limited to 2D data. Basically, every time there is
code trying to build one we should take care of
using CRS.getHorizontalCRS to extract only the
2D part of it before building or transforming
a ReferencedEnvelope (this means, many, many places,
not only datastores, but also wrappers, rendering, and
in GeoServer, UI and configuration).
Either that, or stop using
ReferencedEnvelope altogheter and use a n-dimensional
ready alternative instead
- a few parts of our code try to invert MathTransform
objects. This is usually related to computing the
source bounds to be queries given the destination
bounds. Too bad a 3d -> 2d transformation (the
common case in rendering) is not invertible (due
to information loss).
- the referencing subsystem has issues of its own
trying to build transformations that go from
a compound CRS to another CRS that does not share
the same horizontal component. For example, using the above
compound srs, the following code will fail on
the last statemetn asserting there is no
transformation between f3d and wgs84:
CoordinateReferenceSystem wgs84 = CRS.decode("EPSG:4326");
CoordinateReferenceSystem f2d = CRS.decode("EPSG:27582");
CoordinateReferenceSystem f3d = CRS.decode("EPSG:7401");
System.out.println(CRS.findMathTransform(f3d, f2d));
System.out.println(CRS.findMathTransform(f2d, wgs84));
System.out.println(CRS.findMathTransform(f3d, wgs84));
- the same goes if we try to compute the
transformation between two system with a different
vertical component (that consistently fails), that
always consistently fails
So, in the end I was left wondering, what are people
already dealing with this use case doing?
What are postgis, shapefiles, OGR, GML and CityGML
doing? A cursory investigation shows that a number
of them always reports a 2D srs and provides
3d coordinates, and has a way to tell the user
what the number of dimensions are independent of
the CRS itself.
Examples:
- I have a couple shapefiles provided by people
interested in 2d+1 data support, and generated by
ESRI tools. The .prj contains a 2D srs, the
data in the shapefile is arcz or polyz
- postgis geometry_columns has an explicit
coord_dimension column, and the spatial_ref_sys
table has no 3d reference systems
- gml3 posList has both an explicit srsName and
a srsDimension attribute, some CityGML 1.0
examples available online miss the srsName completely
but do specify the srsDimension one (others do miss
both of them). Mind that the feature collection
envelope is specified in 3d as well.
As an example, see:
http://demo.snowflakesoftware.com:8080/hangar_terrain/GOPublisherWFS?service=wfs&version=1.1.0&request=GetFeature&maxfeatures=10&typename=Building
(32617 is UTM 17N, a 2d CRS)
- ogr2ogr is able to generate 3d kml or 3d gml
out of a xyz shapefile. Given the following
code:
ogr2ogr -f kml -s_srs "EPSG:2907" -t_srs "EPSG:4326"
3dfloor.kml 3Dfloor.shp
the output is a kml file with the z untouched, and
the planar coordinates transformed. Same goes for
GML2 output, that declared an srs of 4326 but uses
coordinates with 3 dimensions
- Oracle spatial started to have real 3d coordinate
systems in 11G, before it did 2.5d with flat CRS
too (and still does in 11G unless you explicitly ask
for full 3D systems).
So, it seems the rest of the world is doing 3D by
declaring 2D srs and just using 3d coordinates...
can't we do the same? ;-)
How do we do the same, thought?
A few ideas to store the dimensions:
- JTS CoordinateSequence does have a getDimension
method. We can use custom coordinate sequences
that allow us to specify the dimensionality
there
- GeometryDescriptor could contain a "well known"
COORDINATE_DIMENSION hint inside its userData
section
For what coordinate handling is concerned,
we can finally make use of the FEATURE_2D hint
we already have to make clients ask for flattened
data when needed, and use 3d otherwise.
As for geometry transformations are concerned, we
should rewrite the transformation code so that it
transforms the flat part but preserves the z and
m if available (the current code will wipe them
out instead).
Opinions?
Cheers
Andrea
--
Andrea Aime
OpenGeo - http://opengeo.org
Expert service straight from the developers.
------------------------------------------------------------------------------
Register Now & Save for Velocity, the Web Performance & Operations
Conference from O'Reilly Media. Velocity features a full day of
expert-led, hands-on workshops and two days of sessions from industry
leaders in dedicated Performance & Operations tracks. Use code vel09scf
and Save an extra 15% before 5/3. http://p.sf.net/sfu/velocityconf
_______________________________________________
Geotools-devel mailing list
Geotools-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/geotools-devel
