5

I am trying to work with a shapefile using shapely, fiona and python. The original coordinate system is in latitude and longitude that I am trying to convert to state plane system. So I want an updated shapefile with coordinates as statePlane. I converted the json object:

fc = fiona.open("sample.shp")   
geoJsonObj = shapefile_record['geometry']
array_coordinates = np.array(geoJsonObj['coordinates'])
array_newCoordinates = copy.deepcopy(array_coordinates)
for counter in range(0,len(array_coordinates[0])):
    long,lat = p1(array_coordinates[0][counter][0],array_coordinates[0][counter][1])
    #where p1 is a function to do the conversion
    array_newCoordinates[0][counter][0] = long
    array_newCoordinates[0][counter][1] = lat

geoJsonObj['coordinates'] = array_newCoordinates.tolist()

when i check the coordinates, i get state plane coordinates according to the transformation.

n. However, when i open the individual shapes/polygons within the shapefile, the coordinates are still latitude and longitude. why does this happen?

----EDIT 1-----

from pyproj import Proj, transform
import fiona
from fiona.crs import from_epsg

dest_crs = from_epsg(4269)
shape = fiona.open("Sample.shp")
original = Proj(shape.crs) # EPSG:4326 in your case
destination = Proj('+proj=lcc +lat_1=36.41666666666666 +lat_2=35.25      +lat_0=34.33333333333334 +lon_0=-86 +x_0=600000 +y_0=0 +ellps=GRS80 +datum=NAD83  +no_defs')
with fiona.open('new.shp', 'w', 'ESRI Shapefile',shape.schema.copy(),crs = dest_crs) as output:
    for feat in shape:       
        print feat['geometry']['coordinates']
        for counter in range(0,len(feat['geometry']['coordinates'][0])):
            long,lat =  feat['geometry']['coordinates'][0][counter]
            x,y = transform(original, destination,long,lat)
            feat['geometry']['coordinates'][0][counter] = (x,y)
        output.write(feat)
  • Although I made an answer, I don't see the point to use Python here. Why don't you use ogr2ogr? Do you know the EPSG, ESRI code for your expected state plane coordinates ? – ThomasG77 Nov 8 '14 at 22:26
  • I didnt know about that. I am very new to this field and therefore, didnt know. – statBeginner Nov 9 '14 at 1:08
  • I had skipped defining the destination crs while opening the new shp file. When i updated the code, the program crashes as soon as i try to open the new.shp file. – statBeginner Nov 9 '14 at 1:20
  • If your file is Sample.shp, the Sample.prj file content (open it with a text editor) associated will help us deduce the exact projection. Can you add it content in your post? – ThomasG77 Nov 9 '14 at 3:03
  • The content is : 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]] – statBeginner Nov 9 '14 at 3:37
13

To change projections with Fiona, use the pyproj module.

Example with a point shapefile (you can simplify the algorithm):

from pyproj import Proj, transform
import fiona
from fiona.crs import from_epsg
shape = fiona.open('sample.shp')
original = Proj(shape.crs) # EPSG:4326 in your case
destination = Proj(init='EPSG:...') # your new EPSG
with fiona.open('new.shp', 'w', 'ESRI Shapefile', shape.schema.copy(), crs=from_epsg(...)) as output:
    for feat in shape:
        long,lat =  feat['geometry']['coordinates']
        x,y = transform(original, destination,long,lat)
        # change only the coordinates of the feature
        feat['geometry']['coordinates'] = (x,y)
        output.write(feat)

If you want to work with the GeoJSON format:

with fiona.open('sample.shp') as source:
    records = list(source)
geojson = {"type": "FeatureCollection","features": records}

And use the same method with the GeoJSON format (with pyproj)

For a polygon, do the same thing with the points of the polygon

for feat in shape:
   for point in feat['geometry']['coordinates'][0]: # = LinearRing of the polygon
       x,y = pt
       ....

new

In your script, why counter?
- feat['geometry']['coordinates'] is a nested list (Polygon)
- feat['geometry']['coordinates'][0] also (LinearRing of the Polygon)
- feat['geometry']['coordinates'][0][x] is a Point of the LinearRing

So:

with fiona.open('new.shp', 'w', 'ESRI Shapefile', shape.schema.copy(), crs=from_epsg(...)) as output:
    for feat in shape: # feat = one polygon of the shapefile
        out_linearRing = [] # empty list for the LinearRing of transformed coordinates
        for point in feat['geometry']['coordinates'][0]: # LinearRing of the Polygon
            long,lat =  point  # one point of the LinearRing
            x,y = transform(original, destination,long,lat) # transform the point
            out_linearRing.append((x,y)) # add all the points to the new LinearRing
        # transform the resulting LinearRing to  a Polygon and write it
        feat['geometry']['coordinates'] = [out_linearRing]
        output.write(feat)
  • Thanks. I was using pyproj. The p in my case was the destination projection (should have mentioned that). Here is what i am using based on your suggestion. But as soon as I execute the statement output.write(feat), the program crashes. please See Edit 1 for updated code – statBeginner Nov 9 '14 at 1:16
  • Thanks for the edit. But some of the problems still persist. If I include crs=from_epsg(4269) while opening the the new.shp, the program crashes. If I omit it, the program crashes while executing output.write(feat) – statBeginner Nov 10 '14 at 23:47
  • Here are the details: Problem Event Name: APPCRASH Application Name: python.exe Application Version: 0.0.0.0 Application Timestamp: 53b1ecd8 Fault Module Name: ntdll.dll Fault Module Version: 6.1.7601.18247 Fault Module Timestamp: 521ea8e7 Exception Code: c0000005 Exception Offset: 00022289 OS Version: 6.1.7601.2.1.0.768.2 Locale ID: 16393 Additional Information 1: 0a9e Additional Information 2: 0a9e372d3b4ad19135b953a78882e789 Additional Information 3: 0a9e Additional Information 4: 0a9e372d3b4ad19135b953a78882e789 – statBeginner Nov 10 '14 at 23:47
  • Can you add a print feat right before output.write(feat) and write back here with the value of the last feature that is printed? – sgillies Nov 11 '14 at 16:35
3

Didn't try but the most obvious is where do you save your new file? Or why if you want to do 'in-place' change, you only open the shp file without changing the mode to be able to write it.

Also more in general, use the with statement like in official docs because for example in your code, you do not close the file whereas with with it's automatic.

  • I will use the with statement. Thank you. How do I save the file? More importantly, does changing the geojsonObj's coordinates in my case ensure that all the polygon's will have their coordinates transformed as well? – statBeginner Nov 8 '14 at 22:02
  • Read the docs, please, writing is already documented toblerity.org/fiona/manual.html#writing-vector-data – ThomasG77 Nov 8 '14 at 22:32
  • In practice, you don't need withto close a shapefile with Fiona. It is recommended but not necessary. – gene Nov 8 '14 at 22:48
  • You're right, IMO, just more Pythonic and a good way to not forget closing file – ThomasG77 Nov 8 '14 at 23:33

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.