16. データの投影

地球は平らでなく、平面の紙地図 (またはコンピュータ画面) に置くための簡単な方法はありません。長所と短所を持つあらゆる種類の独創的な開放を考え出しました。全てのオブジェクトが相互に相対的なサイズを持つように、面積を保存する投影があります。また、メルカトル図法のように角度を保持する (正角) ものがあります。他にも、いくつかのパラメータで歪みがほとんどない、いいとこ取りのものがあります。全ての投影法に共通するのが、球面の世界を平面のデカルト座標系に変換していることと、投影法の選択はデータをどのように使うかに依存することです。

ここまでで ニューヨーク市のデータのロード の時に投影法を目にしています (厄介な SRID 26918 を思い出して下さい)。しかしながら、時々、空間参照系間の変換を行う必要があります。PostGIS は ST_Transform(geometry, srid) 関数を使って、データの投影を変更に対応する機能が組み込まれています。ジオメトリの空間参照系識別子の管理のために PostGIS では ST_SRID(geometry) と ST_SetSRID(geometry, srid) 関数が提供されています。

ST_SRID 関数でデータの SRID を確認することができます:

SELECT ST_SRID(geom) FROM nyc_streets LIMIT 1;

26918

そして "26918" の定義は何でしょう? "データロードのセクション" で見た通り、定義は spatial_ref_sys にあります。事実 二つの 定義があります。"Well-Known Text" (WKT) 定義は srtext カラムにあり、二つ目の "proj.4" フォーマットの定義は proj4text カラムにあります。

SELECT * FROM spatial_ref_sys WHERE srid = 26918;

PostGIS の投影変換エンジンは最適な投影を spatial_ref_sys テーブルから発見しようと試みます:

  • auth_name / auth_srid If proj can find a valid "authority name" and "authority srid" in its internal catalogue, it will use that to generate a projection definition.

  • srtext If proj can parse and form a definition object from the srtext it will use that.

  • proj4text Finally, proj will attempt to process the proj4text.

All this redundancy means that all you need to create a new projection in PostGIS is either a valid srtext string or proj4text string. All the common authority name/code pairs are already loaded in the table by default.

If you have a choice when creating a custom projection, fill out the srtext column, since that column is also used by external programs like GeoServer, QGIS, and FME and others.

16.1. Comparing Data

Taken together, a coordinate and an SRID define a location on the globe. Without an SRID, a coordinate is just an abstract notion. A "Cartesian" coordinate plane is defined as a "flat" coordinate system placed on the surface of Earth. Because PostGIS functions work on such a plane, comparison operations require that both geometries be represented in the same SRID.

If you feed in geometries with differing SRIDs you will just get an error:

SELECT ST_Equals(
         ST_GeomFromText('POINT(0 0)', 4326),
         ST_GeomFromText('POINT(0 0)', 26918)
         );

ERROR:  ST_Equals: Operation on mixed SRID geometries (Point, 4326) != (Point, 26918)

注釈

Be careful of getting too happy with using ST_Transform for on-the-fly conversion. Spatial indexes are built using SRID of the stored geometries. If comparison are done in a different SRID, spatial indexes are (often) not used. It is best practice to choose one SRID for all the tables in your database. Only use the transformation function when you are reading or writing data to external applications.

16.2. Transforming Data

If we return to our proj4 definition for SRID 26918, we can see that our working projection is UTM (Universal Transverse Mercator) of zone 18, with meters as the unit of measurement.

SELECT srtext FROM spatial_ref_sys WHERE srid = 26918;

PROJCS["NAD83 / UTM zone 18N",
  GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
      SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],
      TOWGS84[0,0,0,0,0,0,0],
      AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4269"]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["latitude_of_origin",0],
  PARAMETER["central_meridian",-75],
  PARAMETER["scale_factor",0.9996],
  PARAMETER["false_easting",500000],
  PARAMETER["false_northing",0],
  UNIT["metre",1,AUTHORITY["EPSG","9001"]],
  AXIS["Easting",EAST],AXIS["Northing",NORTH],
  AUTHORITY["EPSG","26918"]]

Let's convert some data from our working projection to geographic coordinates -- also known as "longitude/latitude".

To convert data from one SRID to another, you must first verify that your geometry has a valid SRID. Since we have already confirmed a valid SRID, we next need the SRID of the projection to transform into. In other words, what is the SRID of geographic coordinates?

The most common SRID for geographic coordinates is 4326, which corresponds to "longitude/latitude on the WGS84 spheroid". You can see the definition here:

You can also pull the definitions from the spatial_ref_sys table:

SELECT srtext FROM spatial_ref_sys WHERE srid = 4326;

GEOGCS["WGS 84",
  DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
  UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
  AUTHORITY["EPSG","4326"]]

Let's convert the coordinates of the 'Broad St' subway station into geographics:

SELECT ST_AsText(ST_Transform(geom,4326))
FROM nyc_subway_stations
WHERE name = 'Broad St';

POINT(-74.01067146887341 40.70710481558761)

If you load data or create a new geometry without specifying an SRID, the SRID value will be 0. Recall in ジオメトリ, that when we created our geometries table we didn't specify an SRID. If we query our database, we should expect all the nyc_ tables to have an SRID of 26918, while the geometries table defaulted to an SRID of 0.

To view a table's SRID assignment, query the database's geometry_columns table.

SELECT f_table_name AS name, srid
FROM geometry_columns;

        name         | srid
---------------------+-------
 nyc_census_blocks   | 26918
 nyc_homicides       | 26918
 nyc_neighborhoods   | 26918
 nyc_streets         | 26918
 nyc_subway_stations | 26918
 geometries          |     0

However, if you know what the SRID of the coordinates is supposed to be, you can set it post-facto, using ST_SetSRID on the geometry. Then you will be able to transform the geometry into other systems.

SELECT ST_AsText(
    ST_Transform(
      ST_SetSRID(geom,26918),
      4326)
    )
  FROM geometries;

16.3. 関数リスト

ST_AsText: ジオメトリ/ジオグラフィの SRID メタデータを持たない Well-Known Text (WKT) 表現を返します。

ST_SetSRID(geometry, srid): Sets the SRID on a geometry to a particular integer value.

ST_SRID: ST_Geometry の、spatial_ref_sys テーブルで定義されている空間参照系識別子を返します。

ST_Transform(geometry, srid): 整数パラメータで指定される SRID に座標変換した新しいジオメトリを返します。