source:
trunk/workshop-foss4g/geography.rst
@
48
Revision 48, 13.7 KB checked in by nbozon, 13 years ago (diff) |
---|
Section 17: Coordonnées géographiques
Il est trÚs fréquent de manipuler des données à coordonnées "géographiques" ou de "longitude/latitude".
Au contraire des coordonnées de type Mercator, UTM ou Stateplane, les coordonnées géographiques ne représentent pas une distance linéaire depuis une origine, tel que dans un plan. Elles décrivent la distance angulaire entre l'équateur et les pÎles. Dans les sytÚmes de coordonnées sphériques, un point est spécifié par son rayon (distance à l'origine), son angle de rotation par rapport au méridien plan, et son angle par rapport à l'axe pÎlaire.
Vous pouvez continuer à utiliser des coordonnées géographiques comme des coordonnées cartésiennes approximatives pour vos analyses spatiales. Par contre les mesures de distances, d'aires et de longueur seront éronées. Etant donné que les coordonnées spériques mesurent des angles, l'unité est le dégré. Par exemple, les résultats cartésien approximatifs de tests tels que 'intersects' et 'contains' peuvent s'avérer terriblement faux. Par ailleurs, plus une zone est située prÚs du pÎle ou de la ligne de date internationale, plus la distance entre les points est agrandie.
Voici par exemple les coordonnées des villes de Los Angeles et Paris.
- Los Angeles: POINT(-118.4079 33.9434)
- Paris: POINT(2.3490 48.8533)
La requête suivante calcule la distance entre Los Angeles et Paris en utilisant le systÚme cartésien standard de PostGIS :command:`ST_Distance(geometry, geometry)`. Notez que le SRID 4326 déclare un systÚme de références spatiales géographiques.
SELECT ST_Distance( ST_GeometryFromText('POINT(-118.4079 33.9434)', 4326), -- Los Angeles (LAX) ST_GeometryFromText('POINT(2.5559 49.0083)', 4326) -- Paris (CDG) );
121.898285970107
Aha! 121! Mais, que veut dire cela ?
L'unité pour SRID 4326 est le degré. Donc la réponse signifie 121 degrés. Sur une sphÚre, la taille d'un degré "au carré" est assez variable. Elle devient plsu petite au fur et à mesure que l'on s'éloigne de l'équateur. Pensez par exemple aux méridiens sur le globe qui se ressÚrent entre eux au niveau des pÎles. Donc une distance de 121 degrés ne veut rien dire !
Pour calculer une distance ayant du sens, nous devons traiter les coordonnées géographiques non pas come des coordonnées cartésiennes approximatives, mais plutÎt comme de réelles coordonnées sphériques. Nous devons mesurer les distances entre les points comme de vrais chemins par dessus uen sphÚre, comme une portion d'un grand cercle.
Depuis sa version 1.5, PostGIS fournit cette fonctionnalité avec le type geography.
Note
Différentes bases de données spatiales développent différentes approches pour manipuler les coordonnées géographiques.
- Oracle essaye de mettre à jour la différence de maniÚre transparente en lanacant des calculs lorsuqe le SRID est géographique.
- SQL Server utilise deux types spatiaux, "STGeometry" pour les coordonnées cartésiens et STGeography" pour les coordonnées géographqiues.
- Informix Spatial est une pure extension cartésienne d'Informix, alors qu'Informix Geodetic est une pure extension géographique.
- Comme SQL Server, PostGIS utilise deux types: "geometry" et "geography".
En utilisant le type geography plutot que geometry, essayon sà nouveau de mesurer la distance entre Los Angeles et Paris. Au lieu de la commande :command:`ST_GeometryFromText(text)`, nous utiliserons cette fois :command:`ST_GeographyFromText(text)`.
SELECT ST_Distance( ST_GeographyFromText('POINT(-118.4079 33.9434)'), -- Los Angeles (LAX) ST_GeographyFromText('POINT(2.5559 49.0083)') -- Paris (CDG) );
9124665.26917268
Toutes les valeurs retournées étant en mÚtres, notre réponse est donc 9124 kilomÚtres.
Les versions plus anciennes de PostGIS supportaient uniquement des calculs sur sphÚre trÚs basiques comme la fonction :command:`ST_Distance_Spheroid(point, point, measurement)`. Celle-ci est trÚs limitée et ne fonctionne uniquement sur des points. Elle ne supporte pas non plus l'indexation au niveau des pÎles ou de la ligne de date internationale.
Le besoin du support des autres types de géométries se fit ressentir lorsqu'il s'agissait de répondre à des questions du type "A quelle distance la ligne de vol d'un avion Los Angeles/Paris passe-t-elle de l'Islande?"
Working with geographic coordinates on a cartesian plane (the purple line) yields a very wrong answer indeed! Using great circle routes (the red lines) gives the right answer. If we convert our LAX-CDG flight into a line string and calculate the distance to a point in Iceland using geography we'll get the right answer (recall) in meters.
SELECT ST_Distance( ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)'), -- LAX-CDG ST_GeographyFromText('POINT(-21.8628 64.1286)') -- Iceland );
531773.757079116
So the closest approach to Iceland on the LAX-CDG route is a relatively small 532km.
The cartesian approach to handling geographic coordinates breaks down entirely for features that cross the international dateline. The shortest great-circle route from Los Angeles to Tokyo crosses the Pacific Ocean. The shortest cartesian route crosses the Atlantic and Indian Oceans.
SELECT ST_Distance( ST_GeometryFromText('Point(-118.4079 33.9434)'), -- LAX ST_GeometryFromText('Point(139.733 35.567)')) -- NRT (Tokyo/Narita) AS geometry_distance, ST_Distance( ST_GeographyFromText('Point(-118.4079 33.9434)'), -- LAX ST_GeographyFromText('Point(139.733 35.567)')) -- NRT (Tokyo/Narita) AS geography_distance;
geometry_distance | geography_distance -------------------+-------------------- 258.146005837336 | 8833954.76996256
Using Geography
In order to load geometry data into a geography table, the geometry first needs to be projected into EPSG:4326 (longitude/latitude), then it needs to be changed into geography. The :command:`ST_Transform(geometry,srid)` function converts coordinates to geographics and the :command:`Geography(geometry)` function "casts" them from geometry to geography.
CREATE TABLE nyc_subway_stations_geog AS SELECT Geography(ST_Transform(the_geom,4326)) AS geog, name, routes FROM nyc_subway_stations;
Building a spatial index on a geography table is exactly the same as for geometry:
CREATE INDEX nyc_subway_stations_geog_gix ON nyc_subway_stations_geog USING GIST (geog);
The difference is under the covers: the geography index will correctly handle queries that cover the poles or the international date-line, while the geometry one will not.
There are only a small number of native functions for the geography type:
:command:`ST_AsText(geography)` returns text
:command:`ST_GeographyFromText(text)` returns geography
:command:`ST_AsBinary(geography)` returns bytea
:command:`ST_GeogFromWKB(bytea)` returns geography
:command:`ST_AsSVG(geography)` returns text
:command:`ST_AsGML(geography)` returns text
:command:`ST_AsKML(geography)` returns text
:command:`ST_AsGeoJson(geography)` returns text
:command:`ST_Distance(geography, geography)` returns double
:command:`ST_DWithin(geography, geography, float8)` returns boolean
:command:`ST_Area(geography)` returns double
:command:`ST_Length(geography)` returns double
:command:`ST_Covers(geography, geography)` returns boolean
:command:`ST_CoveredBy(geography, geography)` returns boolean
:command:`ST_Intersects(geography, geography)` returns boolean
:command:`ST_Buffer(geography, float8)` returns geography [1]
:command:`ST_Intersection(geography, geography)` returns geography [1]
Creating a Geography Table
The SQL for creating a new table with a geography column is much like that for creating a geometry table. However, geography includes the ability to specify the object type directly at the time of table creation. For example:
CREATE TABLE airports ( code VARCHAR(3), geog GEOGRAPHY(Point) ); INSERT INTO airports VALUES ('LAX', 'POINT(-118.4079 33.9434)'); INSERT INTO airports VALUES ('CDG', 'POINT(2.5559 49.0083)'); INSERT INTO airports VALUES ('REK', 'POINT(-21.8628 64.1286)');
In the table definition, the GEOGRAPHY(Point) specifies our airport data type as points. The new geography fields don't get registered in the geometry_columns. Instead, they are registered in a new view called geography_columns that is automatically kept up to date without need for an :command:`AddGeom...` like functions.
SELECT * FROM geography_columns;
f_table_name | f_geography_column | srid | type -------------------------------+--------------------+------+---------- nyc_subway_stations_geography | geog | 0 | Geometry airports | geog | 4326 | Point
Note
The ability to define geometry types and SRIDs inside the table CREATE statement, and the automatic update of the geometry_columns metadata are features that have been prototyped with geography and will be added to the geometry type for PostGIS 2.0.
Casting to Geometry
While the basic functions for geography types can handle many use cases, there are times when you might need access to other functions only supported by the geometry type. Fortunately, you can convert objects back and forth from geography to geometry.
The PostgreSQL syntax convention for casting is to append ::typename to the end of the value you wish to cast. So, 2::text with convert a numeric two to a text string '2'. And 'POINT(0 0)'::geometry will convert the text representation of point into a geometry point.
The :command:`ST_X(point)` function only supports the geometry type. How can we read the X coordinate from our geographies?
SELECT code, ST_X(geog::geometry) AS longitude FROM airports;
code | longitude ------+----------- LAX | -118.4079 CDG | 2.5559 REK | -21.8628
By appending ::geometry to our geography value, we convert the object to a geometry with an SRID of 4326. From there we can use as many geometry functions as strike our fancy. But, remember -- now that our object is a geometry, the coordinates will be interpretted as cartesian coordinates, not spherical ones.
Why (Not) Use Geography
Geographics are universally accepted coordinates -- everyone understands what latitude/longitude mean, but very few people understand what UTM coordinates mean. Why not use geography all the time?
- First, as noted earlier, there are far fewer functions available (right now) that directly support the geography type. You may spend a lot of time working around geography type limitations.
- Second, the calculations on a sphere are computationally far more expensive than cartesian calculations. For example, the cartesian formula for distance (Pythagoras) involves one call to sqrt(). The spherical formula for distance (Haversine) involves two sqrt() calls, an arctan() call, four sin() calls and two cos() calls. Trigonometric functions are very costly, and spherical calculations involve a lot of them.
The conclusion?
If your data is geographically compact (contained within a state, county or city), use the geometry type with a cartesian projection that makes sense with your data. See the http://spatialreference.org site and type in the name of your region for a selection of possible reference systems.
If, on the other hand, you need to measure distance with a dataset that is geographically dispersed (covering much of the world), use the geography type. The application complexity you save by working in geography will offset any performance issues. And, casting to geometry can offset most functionality limitations.
Function List
ST_Distance(geometry, geometry): For geometry type Returns the 2-dimensional cartesian minimum distance (based on spatial ref) between two geometries in projected units. For geography type defaults to return spheroidal minimum distance between two geographies in meters.
ST_GeographyFromText(text): Returns a specified geography value from Well-Known Text representation or extended (WKT).
ST_Transform(geometry, srid): Returns a new geometry with its coordinates transformed to the SRID referenced by the integer parameter.
ST_X(point): Returns the X coordinate of the point, or NULL if not available. Input must be a point.
Footnotes
[1] | (1, 2) The buffer and intersection functions are actually wrappers on top of a cast to geometry, and are not carried out natively in spherical coordinates. As a result, they may fail to return correct results for objects with very large extents that cannot be cleanly converted to a planar representation. For example, the :command:`ST_Buffer(geography,distance)` function transforms the geography object into a "best" projection, buffers it, and then transforms it back to geographics. If there is no "best" projection (the object is too large), the operation can fail or return a malformed buffer. |