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().points 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:
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:
If I use the crs.proj4 for the input projection I get:
+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:
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.