Convierta shapefile de coordenadas proyectadas usando Python

Novato aquí luchando con SIG. Estoy tratando de mapear los barrios de la ciudad de Milwuakee usando archivos de formas que se encuentran en el sitio web de su condado . Estoy siguiendo el hilo aquí con cierto éxito. Mi código es da:

from pyproj import Proj, transform
# wisconsing EPSG:32054
# epsg:4326 is for the entire world, wgs 84...not obvious
inProj = Proj(init='epsg:32054')
outProj = Proj(init='epsg:4326')
x1,y1 = 2560131.496875003, 406816.434375003
x2,y2 = transform(inProj,outProj,x1,y1)
print(x2,y2)

con salida,

-65.70220967836329 43.08590211722421

El problema es que esto está mal. El lon / lat para Milwaukee es -87.863984 y 42.920816.

En segundo lugar, ¿cómo puedo hacer esto mediante programación para todo el archivo de forma? Me gustaría trazar esto en el mapa base. Cuando trato de seguir este hilo me sale un código de error:

with fiona.open("ward2012/ward.shp") as shp:
    ori = Proj(init='epsg:32054' ),
    dest= Proj(init='EPSG:4326',preserve_units=True)
    with fiona.open('ward2012/MKE_wards_lat_lon.shp', 'w', 'ESRI Shapefile', shp.schema.copy(), crs=from_epsg(4326))as output:
        for point in shp:
            x,y =  point['geometry']['coordinates']
            point['geometry']['coordinates'] = transform(ori, dest,x,y)
            output.write(point)

error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-139-a5079ab39f99> in <module>()
      4     with fiona.open('ward2012/MKE_wards_lat_lon.shp', 'w', 'ESRI Shapefile', shp.schema.copy(), crs=from_epsg(4326))as output:
      5         for point in shp:
----> 6             x,y =  point['geometry']['coordinates']
      7             point['geometry']['coordinates'] = transform(ori, dest,x,y)
      8             output.write(point)

ValueError: not enough values to unpack (expected 2, got 1)

python shapefile fiona pyproj superhéroe
fuente

Respuestas:

En la primera pregunta, el código 'epsg: 32054' tiene unidades de pies. Por esta razón, es necesario usar 'preserve_units = True' como parámetro en inProj = Proj(init='epsg:32054')línea. Ahora, el siguiente código funciona bien:

from pyproj import Proj, transform
# wisconsing EPSG:32054
# epsg:4326 is for the entire world, wgs 84...not obvious
inProj = Proj(init='epsg:32054', preserve_units=True)
outProj = Proj(init='epsg:4326')
x1,y1 = 2560131.496875003, 406816.434375003
x2,y2 = transform(inProj,outProj,x1,y1)
print (x2,y2)
(-87.9028568836077, 43.09691266312185)

En la segunda pregunta, ward.shp es un archivo de forma poligonal; No es un archivo de forma puntual. En este caso, puede usar el módulo de geopandas para la reproyección. Mi código propuesto es (con mi ruta particular):

import geopandas as gpd

tmp = gpd.GeoDataFrame.from_file('/home/zeito/pyqgis_data/ward2012/ward.shp')

tmpWGS84 = tmp.to_crs({'proj':'longlat', 'ellps':'WGS84', 'datum':'WGS84'})

tmpWGS84.to_file('/home/zeito/pyqgis_data/ward2012/wardWGS84.shp')

En la siguiente imagen, puede ver el archivo de forma reproyectado (wardWGS84.shp) y el punto (-87.9028568836077, 43.09691266312185) antes considerado en Map Canvas de QGIS:

xunilk
fuente