Konvertieren Sie Shapefiles mit Python aus projizierten Koordinaten

Neuling hier kämpft mit GIS. Ich versuche, Schutzzauber für die Stadt Milwuakee mithilfe von Shapefiles zu kartieren, die auf der Website des Landkreises zu finden sind . Ich verfolge den Thread hier mit einigem Erfolg. Mein Code lautet:

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)

mit Ausgabe,

-65.70220967836329 43.08590211722421

Das Problem ist, dass dies falsch ist. Der Lon / Lat für Milwaukee beträgt -87,863984 und 42,920816.

Zweitens, wie kann ich dies programmgesteuert für das gesamte Shapefile tun? Ich möchte dies in eine Grundkarte eintragen. Wenn ich versuche, diesem Thread zu folgen, erhalte ich einen Fehlercode:

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 Superheld
quelle

Antworten:

Bei der ersten Frage hat der Code 'epsg: 32054' Fußeinheiten. Aus diesem Grund muss'erve_units = True 'als Parameter in der inProj = Proj(init='epsg:32054')Zeile verwendet werden. Der nächste Code funktioniert nun gut:

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)

Bei der zweiten Frage ist ward.shp ein Polygon-Shapefile. kein Punkt-Shapefile. In diesem Fall können Sie das Geopandas-Modul für die Neuprojektion verwenden. Mein vorgeschlagener Code lautet (mit meinem speziellen Pfad):

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

Im nächsten Bild sehen Sie das neu projizierte Shapefile (wardWGS84.shp) und den Punkt (-87.9028568836077, 43.09691266312185), bevor Sie ihn auf Map Canvas von QGIS betrachten:

xunilk
quelle