Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am a complete n00b when it comes to GIS, so please bear with me. I started out by posting here: http://stackoverflow.com/questions/6963991/how-can-i-transform-the-coordinates-of-a-shapefile .

In short, I found a shapefile here, and I want to be able to convert that to Latitude/Longitude (Geodetic?) coordinates that I can store in my PostGIS database. I am having trouble working with the Ruby libraries, but that is a separate issue. First, I need to correctly understand the projection.

In the linked zip file, there is a .prj file, which contains the following (It is only one line, but I broke it up for readability):

PROJCS["NAD_1983_StatePlane_California_III_FIPS_0403_Feet",
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",6561666.666666666],PARAMETER["False_Northing",1640416.666666667],
PARAMETER["Central_Meridian",-120.5],PARAMETER["Standard_Parallel_1",37.06666666666667],
PARAMETER["Standard_Parallel_2",38.43333333333333],PARAMETER["Latitude_Of_Origin",36.5],
UNIT["Foot_US",0.3048006096012192]]

I was unable to find tools that translate .prj files into a projection string, so I used this as a reference and came up with:

+proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +units=us-ft +pm=greenwich +no_defs 

I happen to know that the points in this shapefile should correspond to locations in San Francisco, so I could test my work:

cs2cs -f "%.8f" +proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +units=us-ft +pm=greenwich +no_defs +to +proj=lonlat +datum=WGS84 +ellps=WGS84
6011287.4999795845 2100857.2499904726
-164.18117482   17.72087079 0.00000000

That's about 1000km west of Hawaii, so clearly that's not right. I tried many different things, but the closest I got was by stripping the +units=us-ft off.

cs2cs -f "%.8f" +proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +pm=greenwich +no_defs +to +proj=lonlat +datum=WGS84 +ellps=WGS84
6011287.4999795845 2100857.2499904726
-126.98864867   40.47509388 0.00000000

That's about 500 km NW of San Francisco.

What am I missing here? Please educate this desperate newbie.
Worth noting: If I import the shapefile into the online ArcGIS Explorer, it understands it perfectly.

Thanks in advance.

share|improve this question

3 Answers 3

up vote 3 down vote accepted

I believe that the correct string for this projection is:

+proj=lcc +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000002 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

See: http://spatialreference.org/ref/esri/102643/

Regards, Nick.

share|improve this answer
    
Just to clarify, the proj4 addendum ftp.remotesensing.org/proj/proj.4.3.pdf (p. 7) explains that the + parameters always need to be specified in metres, regardless of the units the projection is in. That's why the +x and +y definitions in Nick's answer are different, because Larry's parameters were in us-ft. This is something that has caught me out in the past. –  MerseyViking Aug 10 '11 at 9:28
    
Yes! This seemed to work. Thank you SO much. And thanks for the explanation MerseyViking. –  Larry Aug 10 '11 at 15:55

Coordinate reference systems and transformations are a complicated matter and are a source of trouble in practice.

I'm not aware of any tool that converts WKT to proj4. As a general hint, look up the PROJCS keyworks on the Spatial Reference website (that's how nhopton found his answer) to find hints to your projection.

share|improve this answer
1  
While there isn't a GDAL tool as such, it is four lines of Python to convert from WKT to Proj4. Posting code in a comment isn't easy, but you just need to create a osr.SpatialReference object, call ImportFromWkt() passing in the WKT string, and then call ExportToProj4() and out pops the answer. In this case, the PROJECTION parameter needs to be changed from Lambert_Conformal_Conic to Lambert_Conformal_Conic_2SP, but that is a corner case in an otherwise very robust implementation. –  MerseyViking Aug 10 '11 at 13:56
    
Great! Post it as a separate answer! –  EPSG31468 Aug 10 '11 at 15:36
    
I know I'm probably pushing my luck here, but if anybody has used RGeo for Ruby, please take a look at stackoverflow.com/questions/7006581/… –  Larry Aug 10 '11 at 15:57

As requested by EPSG31468, here's a very simple Python script to convert from any projection format handled by GDAL to Proj4 strings:

#!/usr/bin/python
# This code is released into the public domain, but may be subject to the licence(s) of the third-party libraries that have been used.
# There is absolutely no warranty or guarantee of fitness for purpose.

import sys
import osgeo.osr as osr

projectionString = sys.stdin.read()

sr = osr.SpatialReference()
sr.SetFromUserInput(projectionString)
proj4 = sr.ExportToProj4()

if proj4 == "":
    print("Unable to convert input string.")
else:
    print(proj4)

There's not a great deal of error handling, but it works with the strings I've tried it with.

To use it, type something like:

$ cat larrysWKT.txt | toproj4.py
+proj=lcc +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000001 +ellps=GRS80 +units=us-ft +no_defs

or even:

$ echo EPSG:27700 | toproj4.py
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

Equivalently, Windows users will need to do this:

c:\> type larrysWKT.txt | python toproj4.py
c:\> echo EPSG:27700 | python toproj4.py

Naturally you'll need Python installed, as well as GDAL/OGR and the corresponding Python bindings. Linux users should be able to get this from their package repositories, and I think Windows users can install it from the OSGeo4W installer.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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