地图投影：如何转换 shapefile

问题描述

我从 Natural Earth 下载了“Admin 0- countries (Download countries 4.67MB v4.1.0)”数据集。如果我使用 Aspose SHP Viewer 查看下载的数据集，我会看到：

但是，我希望能够转换每个坐标，以便我可以创建自己的地图投影。我该怎么做？

我查看了 Fiona 和 pyshp，但我不确定如何才能做到这一点。

import fiona
import matplotlib
import matplotlib.pyplot as plt
from math import log,tan,pi,radians


def mercator_proj(lon,lat):
    lat,lon = map(radians,[lat,lon])

    R = 6371000
    lon_0 = 0
    x = R*(lon-lon_0)
    y = R*log(tan(pi/4+lat/2))
    return (x,y)

shape = fiona.open("ne_10m_admin_0_countries.shp")
# print(shape.schema)

k = iter(shape)

x,y = [],[]
for i in range(150): #first 150 countries
    s = next(k)
    for lis in s["geometry"]["coordinates"]:
        for thing in lis[0]: #it always appears as a nested list of single element,[[]]
            if not isinstance(thing,tuple): #occasionally there is a random float
                continue
            _x,_y = mercator_proj(*thing)
            x.append(_x)
            y.append(_y)
            # print(_x,_y)

colors = ["green"]
color_indices = [0 for i in range(len(x))]
colormap = matplotlib.colors.ListedColormap(colors)
plt.scatter(x,y,c=color_indices,cmap=colormap)
plt.show()

产生以下结果：

相反，我希望能够复制上面的图片。它也应该适用于其他预测（而不仅仅是 mercator_proj）。

编辑： “我希望能够转换每个坐标，以便我可以创建自己的地图投影”，我的意思是我希望能够使用我自己的数学定义公式。大多数地图投影都有公式，例如下面的墨卡托投影：

这需要一个纬度和经度，并返回 2D 地图上的 x,y 值。

解决方法

我不是 100% 确定您要做什么，但我相信您可以这样做：

import geopandas as gpd

# Location of the zip file on disk
shp_file = 'ne_10m_admin_0_countries.zip'

# Loading the data into a GeoDataFrame
my_geodata = gpd.read_file(shp_file)

# Transforming the data's projection/coordinate system 
# to EPSG 4326 (which is the default system used in the 
# Aspose viewer you posted.)
my_geodata = my_geodata.to_crs('epsg:4326')

# Plotting the data
my_geodata.plot()

上面的代码产生下图：

请注意，您需要安装 GeoPandas 库。它在后台使用 Fiona，可以轻松处理从一个投影到另一个投影的转换。

它还使操作数据库的其他部分变得非常容易。