Convertir shapefile de coordenadas proyectadas utilizando python

Question

Novato en apuros con el SIG. Estoy tratando de trazar el mapa de los barrios de la ciudad de Milwuakee el uso de los shapefiles de encontrar en su sitio web del condado de condado de sitio web. Estoy siguiendo el hilo aquí con cierto éxito. Mi código es:

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. La 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 parcela esta en el mapa base. Cuando intento seguir este hilo me sale un error Código:

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)

Answer 1

En la primera pregunta, "epsg:32054' código pies de unidades. Por esta razón, es necesario el uso de '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, de barrio.shp es un polígono shapefile; no es un punto de shapefile. En este caso, puede utilizar geopandas módulo para la reproyección. Mi propuesta de código es (con mi particular ruta de acceso):

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, se puede ver que se vuelve a proyectar shapefile (wardWGS84.shp) y el punto (-87.9028568836077, 43.09691266312185) antes de considerarse en el Mapa de Lona de QGIS: