3

I'm working with las/laz files which often have GeoTIFF information attached to them. However, for further processing I need proj4 definitions. http://spatialreference.org/ has helped me so far to convert between both manualy but ideally I'd like to do that automatically. Is this possible somehow?

Here 2 examples of the GeoTIFF information I have and the resulting proj4 string according to spatialreference.org

GeoTIFF:

  reserved             43707
  user ID              'LASF_Projection'
  record ID            34735
  length after header  72
  description          'GeoTIFF GeoKeyDirectoryTag'
    GeoKeyDirectoryTag version 1.1.0 number of keys 8
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 1025 tiff_tag_location 0 count 1 value_offset 1 - GTRasterTypeGeoKey: RasterPixelIsArea
      key 1026 tiff_tag_location 34737 count 14 value_offset 0 - GTCitationGeoKey: CH1903 / LV03
      key 2049 tiff_tag_location 34737 count 7 value_offset 14 - GeogCitationGeoKey: CH1903
      key 2054 tiff_tag_location 0 count 1 value_offset 9102 - GeogAngularUnitsGeoKey: Angular_Degree
      key 2062 tiff_tag_location 34736 count 3 value_offset 0 - key ID 2062 not implemented
      key 3072 tiff_tag_location 0 count 1 value_offset 21781 - ProjectedCSTypeGeoKey: CH1903 / LV03
      key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
variable length header record 2 of 4:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34736
  length after header  24
  description          'GeoTIFF GeoDoubleParamsTag'
    GeoDoubleParamsTag (number of doubles 3)
      674.4 15.1 405.3
variable length header record 3 of 4:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34737
  length after header  22
  description          'GeoTIFF GeoAsciiParamsTag'
    GeoAsciiParamsTag (number of characters 22)
      CH1903 / LV03|CH1903|
variable length header record 4 of 4:
  reserved             43707
  user ID              'liblas'
  record ID            2112
  length after header  715
  description          'OGR variant of OpenGIS WKT SRS'

proj4:

"+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +units=m +no_defs"

GeoTIFF:

  description          'GeoTIFF GeoKeyDirectoryTag'
    GeoKeyDirectoryTag version 1.1.0 number of keys 11
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 1025 tiff_tag_location 0 count 1 value_offset 1 - GTRasterTypeGeoKey: RasterPixelIsArea
      key 1026 tiff_tag_location 34737 count 31 value_offset 0 - GTCitationGeoKey: NAD83 / UTM zone 10N + VERT_CS
      key 2049 tiff_tag_location 34737 count 6 value_offset 31 - GeogCitationGeoKey: NAD83
      key 2054 tiff_tag_location 0 count 1 value_offset 9102 - GeogAngularUnitsGeoKey: Angular_Degree
      key 3072 tiff_tag_location 0 count 1 value_offset 26910 - ProjectedCSTypeGeoKey: UTM 10 northern hemisphere
      key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
      key 4096 tiff_tag_location 0 count 1 value_offset 5703 - VerticalCSTypeGeoKey: NAVD88 height (Reserved EPSG)
      key 4097 tiff_tag_location 34737 count 14 value_offset 37 - VerticalCitationGeoKey: NAVD88 height
      key 4098 tiff_tag_location 0 count 1 value_offset 5103 - VerticalDatumGeoKey: Vertical Datum Codes 5103
      key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
variable length header record 2 of 3:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34737
  length after header  51
  description          'GeoTIFF GeoAsciiParamsTag'
    GeoAsciiParamsTag (number of characters 51)
      NAD83 / UTM zone 10N + VERT_CS|NAD83|NAVD88 height|
variable length header record 3 of 3:
  reserved             43707
  user ID              'liblas'
  record ID            2112
  length after header  951
  description          'OGR variant of OpenGIS WKT SRS'

proj4:

"+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
3

I'm assuming you want to get the proj4 string by just specifying the EPSG code.

You can get the proj4 string with Python.

import urllib2
EPSG = 21781
proj4_url = 'http://spatialreference.org/ref/epsg/{0}/proj4/'
proj4_str = urllib2.urlopen(proj4_url.format(EPSG)).read()
print proj4_str

The result is

+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs

Just change the EPSG code to get another proj4 string.

3

You can also directly use pyepsg, A simple interface to epsg.io

import pyepsg
projection = pyepsg.get(21781)
pro4_str = projection.as_proj4()
print pro4_str
+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +units=m +no_defs 
wkt_str = projection.as_wkt()
print wkt_str
PROJCS["CH1903 / LV03",GEOGCS["CH1903",DATUM["CH1903",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[674.4,15.1,405.3,0,0,0,0],AUTHORITY["EPSG","6149"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4149"]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",600000],PARAMETER["false_northing",200000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Y",EAST],AXIS["X",NORTH],AUTHORITY["EPSG","21781"]]
2

Python is probably best for automatic system but for temporary needs I suggest to try gdalsrsinfo http://www.gdal.org/gdalsrsinfo.html because with that you can read the the projection info directly from your images. Program supports seven formats for output:

  • proj4 PROJ.4 string
  • wkt OGC WKT format (full)
  • wkt_simple OGC WKT (simplified)
  • wkt_noct OGC WKT (without OGC CT params)
  • wkt_esri ESRI WKT format
  • mapinfo Mapinfo style CoordSys format
  • xml XML format (GML based)

Example: Find proj4 definitions for image "32bit.tif"

gdalsrsinfo -o proj4 32bit.tif
'+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0 =-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs '

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.