[GIS] Converting coordinates from shapefile using python

pyprojpyshppythonshapefile

I am really new to working with this kind of data.

I have been trying to work with a shapefile that I wish to convert to geojson. The source has both shapefiles and geojson downloads available, but the geojson files they have available are truncated so are unusable, but do contain enough information for me to compare my script's output to expected values.

I am using a combination of pyshp to load the shapes/records, pycrs to load the .prj file, and pyproj to preform the transformations.

The simplest check I run is something like:

sf  = shapefile.Reader(fileName)
crs = pycrs.loader.from_file(fileName + ".prj")
input_projection = Proj(crs.to_proj4()) # or init="ESRI:102718"
output_projection = Proj(init="epsg:4326")
x,y     =  sf.shapes()[0].points[0]
lo,la   =  transform(input_projection,output_projection,x,y)

I have tried both using the proj4 string, and the ESRI code (found via spatialreference.org), both seem to be producing different incorrect results.

Example: The first point of the first shape in the file has points:

(994133.507019,214848.897583)

This should, according to the truncated geojson (which does line up correctly on google maps so I believe its points are correct) produce something like:

(-73.964326872592096, 40.756389797390298)

If I use the crs.proj4 for the input projection I get:

(-73.8805457169,42.1010909572)

Exact string:

+proj=lcc +datum=NAD83 +ellps=GRS80 +a=6378137.0 +f=298.257222101 +pm=0.0  +x_0=984250.0 +y_0=0.0 +lon_0=-74.0 +lat_1=40.6666666667 +lat_2=41.0333333333 +lat_0=40.1666666667 +units=us-ft +axis=enu +no_defs

While if I use the ESRI code for the input projection I end up with:

(-65.6360519953,41.8026275179)

When I plot the whole set, the shapes end up 'above' their correct positions and several times larger than they should be, but otherwise undistorted relative to each other.

So I am clearly doing something wrong or not understanding something important, but am at a loss for where I am going so horribly wrong.

Any ideas?

Best Answer

Looks like pyproj has some issue working in imperial units, so a shapefile set in feet does not convert quite right.

The issue can be fixed by something like:

input_projection = Proj(init="ESRI:102718", preserve_units=True)

Curiously enough neither of the proj4 strings produced by:

pycrs.parser.from_esri_code(102718).to_proj4()

+proj=lcc +datum=NAD83 +ellps=GRS80 +a=6378137.0 +f=298.257222101 +pm=0.0  +x_0=984250.0 +y_0=0.0 +lon_0=-74.0 +lat_1=40.6666666667 +lat_2=41.0333333333 +lat_0=40.1666666667 +units=us-ft +axis=enu +no_defs

nor:

pycrs.loader.from_file(fileName + ".prj").crs.to_proj4())

+proj=lcc +datum=NAD83 +ellps=GRS80 +a=6378137.0 +f=298.257222101 +pm=0  +lon_0=-74 +x_0=300000 +y_0=0 +lat_0=40.16666666666666 +to_meter=0.3048006096012192 +axis=enu +no_defs

manage to produce correct results. So it is possible pycrs is somehow borked.

Related Question