Converting coordinates from shapefile using python?

Question

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?

• ZrSiO4: If your main goal is to convert a shapefile to GeoJSON, this can be done with QGIS (FOSS) by simply saving the shapefile as GeoJSON or with GeoPandas if you want to use code (geopandas.org/io.html?highlight=geojson).
• Michael Stimson: GeoJSON is an OGR format gdal.org/drv_geojson.html you should be able to convert in less than 10 lines of python if you start with one of the GDAL python bindings. (GDAL and OGR are fairly much the same lib nowdays, it used to be GDAL was raster and OGR was vector but now they're distributed together under the same package). The OSGeo installed package is fairly stable. Tutorial gdal.org/ogr_apitut.html (there is some python there amongst the C/C++).

Answer 1

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.