Converting while reading a shapefile

Question

I have a set of functions I am using to read shapefiles and produce some geojson files for parsing:

def read_shapefile(path, buffer = []):
    # read the shapefile
    reader = shapefile.Reader(path)
    fields = reader.fields[1:]
    field_names = [field[0] for field in fields]
    #buffer = []
    for sr in reader.shapeRecords():
        for i in range(len(sr.record)):
            if type(sr.record[i]) == bytes:
                sr.record[i] = str(sr.record[i])
        atr = dict(zip(field_names, sr.record))
        geom = sr.shape.__geo_interface__
        buffer.append(dict(type="Feature", \
        geometry=geom, properties=atr))
    return buffer

def write_geojson(buffer):
    # write the GeoJSON file
    #print(buffer)
    print(type(buffer))
    geojson = open("pyshp-demo.json", "w")
    geojson.write(dumps({"type": "FeatureCollection",\
    "features": buffer}, indent=2) + "\n")
    geojson.close()

def read_shapes():
    """ read all of the shapes in our shape directory and write them to geojson"""
    files = get_file_list(SHAPE_DIRECTORY, [".shp"])
    print(files)
    buffer = []
    for file in files:
        buffer = read_shapefile(file, buffer)
    write_geojson(buffer)

which works fine but sometimes I am getting files in X\Y instead of lat long. I see How to convert projected coordinates to lat/lon using Python? which is useful for converting single points but is there any way inside read_shapefile I could easily apply the conversion to the whole geometry object instead of parsing through my end result line by line?

could you clarify, do you want to reproject from one CRS to another, or to swap the x and y values? — Steven Kay, Feb 28 '17 at 21:17
I have values in a projection/coordinate system other then lat long and I am trying to convert the values to lat long — lathomas64, Mar 01 '17 at 18:43

score 7 · Answer 1 · answered Mar 06 '17 at 23:39

Using ogr2ogr

The simplest way to do this is probably to just use ogr2ogr, no need for code. Something like this:-

ogr2ogr -f GeoJSON -s_srs XXXX -t_srs 4326 output.geojson input.shp

Replace XXXX with the CRS of your shapefile.

In Python

If you want to use Python, you can do this using shapely / fiona / pyproj. This will be slower than ogr2ogr, but more flexible.

This code will reproject a polygon shapefile and write the results as geoJSON. It doesn't work with points at the moment.

import json
import fiona
import pyproj
from shapely.geometry import asShape, mapping
from shapely.ops import transform
from functools import partial

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:XXXX'), # source coordinate system, you need to set this
    pyproj.Proj(init='epsg:4326'), # destination coordinate system
)

output = {"type": "FeatureCollection", "features": []}

with open("/path/to/output.geojson","w") as file_out:
    with fiona.open("/path/to/input.shp", "r") as source:
        for feature in source:
            try:
                geom = asShape(feature["geometry"])
                feature["geometry"] = mapping(transform(project, geom))
                output["features"].append(feature)
            except:
                pass # i'll leave this an exercise for you ;)
    file_out.write("{}".format(json.dumps(output, indent=4, sort_keys=True)))

print("Finished")

This holds the whole geojson in memory, so you may need to refactor it for very large shapefiles.

In GeoJSON, coordinates are always in X,Y format (lon,lat) rather than Y,X format (lat,lon).

Converting while reading a shapefile

1 Answers1