1

I am using Python to reproject a shapefile (called "SCL" in my project) into a new coordinate system that will allow me to plot it on the same map as another shapefile (called "PUMA" in my project). After the transformation, my resulting plot looks like this:

enter image description here

This issue is that red area should be overlapping with the blue area around -122 long and -47 lat, but it's shifted right and seems elongated as well. This makes me think something is going wrong with my projection process of converting the red shape ("SCL") into the same system as the blue shape ("PUMA").

I've tried for hours to no avail... Can someone take a look at my code and see what might be going wrong?

1. LIBRARIES LOADED:

    import matplotlib.pyplot as plt
    from shapely.geometry import Polygon, MultiPolygon, shape, mapping
    import fiona
    from descartes import PolygonPatch
    from pyproj import Proj, transform
    import shapefile

2. GETTING THE TRANSFORMATION SYSTEMS:

    puma_shape = fiona.open("PUMA_2016_boundaries\cb_2016_53_puma10_500k.shp")
    print(puma_shape.crs) 

Gives {'init': 'epsg:4269'}

    scl_shape = fiona.open("SCL_Service_Territory_Polygon\SCL_Service_Territory_Polygon.shp")
    print(scl_shape.crs) 

Gives {'proj': 'lcc', 'lat_1': 47.5, 'lat_2': 48.73333333333333, 'lat_0': 47, 'lon_0': -120.8333333333333, 'x_0': 500000.0000000001, 'y_0': 0, 'ellps': 'GRS80', 'towgs84': '0,0,0,0,0,0,0', 'units': 'us-ft', 'no_defs': True}

    original = Proj(scl_shape.crs)
    destination = Proj(puma_shape.crs)

3. TRANSFORMING THE COORDINATES

Based entirely off of: https://glenbambrick.com/2016/01/24/reproject-shapefile/

I'm pretty sure is correct / doesn't need revising.

    shpf = shapefile.Reader("SCL_Service_Territory_Polygon\SCL_Service_Territory_Polygon.shp")
    wgs_shp = shapefile.Writer(shapefile.POLYGON)
    fields = shpf.fields
    wgs_fields = wgs_shp.fields
    for name in fields: 
        if type(name) == "tuple":
        continue
    else: 
        args = name
        wgs_shp.field(*args)

    records = shpf.records()
    for row in records:
        args = row
        wgs_shp.record(*args)
    geom = shpf.shapes()
    for feature in geom:
        if len(feature.parts) == 1: 
                poly_list = []
                for coords in feature.points: 
                    long,lat=coords[0],coords[1]
                    new_x,new_y = transform(original, destination,long,lat) # TRANSFORMATION! 
                    poly_coord = [float(new_x),float(new_y)]
                    poly_list.append(poly_coord)
                wgs_shp.poly(parts=[poly_list])
            else:
                print('multipart!')
                feature.parts.append(len(feature.points))
                poly_list = []
                parts_counter = 0
                while parts_counter <len(feature.parts)-1:
                    coord_count = feature.parts[parts_counter]
                    no_of_points = abs(feature.parts[parts_counter]-feature.parts[parts_counter+1])
                part_list = []
                end_point = coord_count + no_of_points
                while coord_count < end_point:
                    for coords in feature.points[coord_count:end_point]:
                        long,lat = coords[0],coords[1]
                        new_x,new_y = transform(original, destination,long,lat) # TRANSFORMATION! 
                        poly_coord = [float(new_x),float(new_y)]
                        part_list.append(poly_coord)
                        coord_count = coord_count + 1
                    poly_list.append(part_list)
                parts_counter = parts_counter + 1
            wgs_shp.poly(parts=poly_list)
    wgs_shp.save("SCL_Service_Territory3.shp")

4. PLOTTING THE TWO SHAPEFILES:

    scl_shape2 = fiona.open("SCL_Service_Territory3.shp")
    fig = plt.figure()
    ax = fig.gca()
    BLUE = '#6699cc'

    for poly in all_polygons:
        ax.add_patch(PolygonPatch(poly['geometry'], fc=BLUE, ec=BLUE, alpha=0.5,zorder=2))

    for poly in all_polygons2:
        ax.add_patch(PolygonPatch(poly['geometry'], fc='RED', ec="RED", alpha=0.5,zorder=2))
    ax.axis('scaled')
    plt.show()

This produces the image shared above. Any idea what's going wrong?

marked as duplicate by AndreJ, Andre Silva, xunilk, Jochen Schwarze, mkennedy Sep 3 '18 at 22:39

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

1

Figured it out. It was a simple issue with units!

Instead of:

original = Proj(scl_shape.crs)

I should have been doing:

original = Proj(scl_shape.crs,preserve_units=True)

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