Python 和 Shapefile：导入 shapefile 后的坐标非常大

Question

我下载了波士顿的 shapefile 并想使用下面的代码绘制它。但是它给我一个错误 ValueError: lat_0 must be between -90.000000 and 90.000000 degrees

原来 coords 有值 (33869.92130000144, 777617.2998000011, 330800.31099999696, 959741.1853) 为什么这么大？

代码

# Import Boston shapefile
shapefilename = 'ZIPCODES_NT_POLY'
shp = fiona.open(shapefilename + '.shp')
coords = shp.bounds
shp.close()

w, h = coords[2] - coords[0], coords[3] - coords[1]
extra = 0.01

m = Basemap(
    projection='tmerc', ellps='WGS84',
    lon_0 = np.mean([coords[0], coords[2]]), 
    lat_0 = np.mean([coords[1], coords[3]]), 
    llcrnrlon = coords[0] - extra * w,
    llcrnrlat = coords[1] - extra * h,
    urcrnrlon = coords[2] + extra * w,
    urcrnrlat = coords[3] + extra * h,
    resolution = 'i', suppress_ticks = True)

错误

ValueError: lat_0 must be between -90.000000 and 90.000000 degrees

Answer 1

您需要从东距和北距的原生投影重新投影到另一个坐标参考系。如果你想要纬度和经度，通常是 WGS84 or EPSG:4326。重新投影的方法如下：

import fiona
import pyproj
from functools import partial
from shapely.geometry import box
from shapely.ops import transform

shp = fiona.open('ZIPCODES_NT_POLY.shp', 'r')
p_in = pyproj.Proj(shp.crs)
bound_box = box(*shp.bounds)
shp.close()

p_out = pyproj.Proj({'init': 'EPSG:4326'})  # aka WGS84
project = partial(pyproj.transform, p_in, p_out)

bound_box_wgs84 = transform(project, bound_box)

print('native box: ' + str(bound_box))
print('WGS84 box: ' + str(bound_box_wgs84))

native box: POLYGON ((330800.310999997 777617.2998000011, 330800.310999997 959741.1853, 33869.92130000144 959741.1853, 33869.92130000144 777617.2998000011, 330800.310999997 777617.2998000011))

WGS84 box: POLYGON ((-69.93980848410942 41.23787282321487, -69.899038698261 42.87724537285449, -73.53324195423785 42.8704709990465, -73.48147096070339 41.2312695091688, -69.93980848410942 41.23787282321487))

否则，Basemap需要的参数大部分都在shp.crs（看一看）

Answer 2

您必须更改坐标参考系。我推荐 geopandas:

import geopandas as gpd
import matplotlib.pyplot as plt

# load zips with the source projection
shapefilename = 'ZIPCODES_NT_POLY'
zips = gpd.read_file(shapefilename + '.shp')
# convert projection to familiar lat/lon
zips = zips.to_crs('epsg:4326')

zips.plot()
plt.show()

Python 和 Shapefile：导入 shapefile 后的坐标非常大

Python and Shapefile: Very large coordinates after importing shapefile

python

gis

shapefile

python-2.7

matplotlib-basemap