Rappresentazioni di SR

Esistono varie rappresentazioni di un sistema di riferimento spaziale. Ad es il sistema di riferimento WGS84 puo' essere rappresentato in questi modi (ed altri ancora):

EPSG 4326

WKT 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"]]

PROJ4 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

GML <gml:GeographicCRS gml:id="ogrcrs1">\n <gml:srsName>WGS 84</gml:srsName>\n <gml:srsID>\n <gml:name gml:codeSpace="urn:ogc:def:crs:EPSG::">4326</gml:name>\n </gml:srsID>\n <gml:usesEllipsoidalCS>\n <gml:EllipsoidalCS gml:id="ogrcrs2">\n <gml:csName>ellipsoidal</gml:csName>\n <gml:csID>\n <gml:name gml:codeSpace="urn:ogc:def:cs:EPSG::">6402</gml:name>\n </gml:csID>\n <gml:usesAxis>\n <gml:CoordinateSystemAxis gml:id="ogrcrs3" gml:uom="urn:ogc:def:uom:EPSG::9102">\n <gml:name>Geodetic latitude</gml:name>\n <gml:axisID>\n <gml:name gml:codeSpace="urn:ogc:def:axis:EPSG::">9901</gml:name>\n </gml:axisID>\n <gml:axisAbbrev>Lat</gml:axisAbbrev>\n <gml:axisDirection>north</gml:axisDirection>\n </gml:CoordinateSystemAxis>\n </gml:usesAxis>\n <gml:usesAxis>\n <gml:CoordinateSystemAxis gml:id="ogrcrs4" gml:uom="urn:ogc:def:uom:EPSG::9102">\n <gml:name>Geodetic longitude</gml:name>\n <gml:axisID>\n <gml:name gml:codeSpace="urn:ogc:def:axis:EPSG::">9902</gml:name>\n </gml:axisID>\n <gml:axisAbbrev>Lon</gml:axisAbbrev>\n <gml:axisDirection>east</gml:axisDirection>\n </gml:CoordinateSystemAxis>\n </gml:usesAxis>\n </gml:EllipsoidalCS>\n </gml:usesEllipsoidalCS>\n <gml:usesGeodeticDatum>\n <gml:GeodeticDatum gml:id="ogrcrs5">\n <gml:datumName>WGS_1984</gml:datumName>\n <gml:datumID>\n <gml:name gml:codeSpace="urn:ogc:def:datum:EPSG::">6326</gml:name>\n </gml:datumID>\n <gml:usesPrimeMeridian>\n <gml:PrimeMeridian gml:id="ogrcrs6">\n <gml:meridianName>Greenwich</gml:meridianName>\n <gml:meridianID>\n <gml:name gml:codeSpace="urn:ogc:def:meridian:EPSG::">8901</gml:name>\n </gml:meridianID>\n <gml:greenwichLongitude>\n <gml:angle gml:uom="urn:ogc:def:uom:EPSG::9102">0</gml:angle>\n </gml:greenwichLongitude>\n </gml:PrimeMeridian>\n </gml:usesPrimeMeridian>\n <gml:usesEllipsoid>\n <gml:Ellipsoid gml:id="ogrcrs7">\n <gml:ellipsoidName>WGS 84</gml:ellipsoidName>\n <gml:ellipsoidID>\n <gml:name gml:codeSpace="urn:ogc:def:ellipsoid:EPSG::">7030</gml:name>\n </gml:ellipsoidID>\n <gml:semiMajorAxis gml:uom="urn:ogc:def:uom:EPSG::9001">6378137</gml:semiMajorAxis>\n <gml:secondDefiningParameter>\n <gml:inverseFlattening gml:uom="urn:ogc:def:uom:EPSG::9201">298.257223563</gml:inverseFlattening>\n </gml:secondDefiningParameter>\n </gml:Ellipsoid>\n </gml:usesEllipsoid>\n </gml:GeodeticDatum>\n </gml:usesGeodeticDatum>\n</gml:GeographicCRS>\n'

Un modo molto rapido per ottenere una di queste rappresentazioni partendo da un'altra e' con l'API di GDAL/OGR. Ad es il seguente script di Python converte una rappresentazione PROJ4 a WKT:

   import sys
   from osgeo import osr
   def proj2wkt(proj4_string):
       srs = osr.SpatialReference()
       srs.ImportFromProj4(proj4_string)
       print srs.ExportToWkt()
   proj2wkt(sys.argv[1])

salvando questo semplice script ad e.s come proj2wkt.py, e lanciandolo da linea di comando:

$ python proj2wkt.py '+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +no_def' +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +no_def

output: PROJCS["unnamed",GEOGCS["International 1909 (Hayford)",DATUM["unknown",SPHEROID["intl",6378388,297],TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",2520000],PARAMETER["false_northing",0],UNIT["Meter",1]]

ovviamente e' facile modificarlo per avere anche altre rappresentazioni, oppure creare un'utility che consenta da linea di comando di ottenere una qualsiasi rappresentazione.

Notare che GRASS consente di ottenere la stessa informazione utilizzando il comando g.proj

Come ottenere le rappresentazioni del SR a partire dal file prj di uno shapefile

Lo shapefile ha la sua proiezione definita nel file prj (se presente). A volte e' utile ricondurre questa definizione ad una rappresentazione standard, ad es WKT o EPSG. Questo semplice script Python, usando l'API di GDAL/OGR, effettua proprio questa operazione:

import sys from osgeo import osr

def esriprj2standards(shapeprj_path):

  prj_file = open(shapeprj_path, 'r')
  prj_txt = prj_file.read()
  srs = osr.SpatialReference()
  srs.ImportFromESRI([prj_txt])
  print 'Shape prj is: %s' % prj_txt
  print 'WKT is: %s' % srs.ExportToWkt()
  print 'Proj4 is: %s' % srs.ExportToProj4()
  srs.AutoIdentifyEPSG()
  print 'EPSG is: %s' % srs.GetAuthorityCode(None)

esriprj2standards(sys.argv[1])

a questo punto eseguendo lo script dandogli in pasto un file .prj:

$ python esriprj2standards.py /home/pcorti/data/shapefile/country.prj Shape prj is: GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] WKT is: GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] Proj4 is: +proj=longlat +datum=WGS84 +no_defs EPSG is: 4326