Convertir mientras lee un shapefile

Question

Tengo un conjunto de funciones que estoy usando para leer los shapefiles y producir algunos geojson archivos para su análisis:

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)

que funciona bien, pero a veces estoy recibiendo los archivos en X\Y en lugar de lat long. Veo Cómo convertir coordenadas proyectadas de lat/lon usando Python? cual es útil para la conversión de puntos pero ¿hay alguna manera dentro de read_shapefile fácilmente se podría aplicar la conversión a la totalidad de la geometría del objeto en vez de analizar a través de mi resultado final de la línea por línea?

Answer 1

El uso de ogr2ogr

La forma más sencilla de hacer esto es probablemente sólo uso ogr2ogr, sin necesidad de código. Algo como esto:-

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

Reemplazar XXXX con el CRS del shapefile.

En Python

Si desea utilizar Python, usted puede hacer esto utilizando shapely / fiona / pyproj. Esto será más lento que ogr2ogr, pero más flexible.

Este código se reproyecta un shapefile de polígonos y escribir los resultados como geoJSON. No funciona con los puntos en el momento.

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")

Esto es todo geojson en la memoria, por lo que usted puede necesitar para refactorizar es muy grande shapefiles.

En GeoJSON, las coordenadas son siempre en X,Y (formato de lon,lat) en lugar de Y,X formato (lat,lon).