将球形地形映射到 Geopandas 和 Cartopy蛋白质表面制图映射

如何解决将球形地形映射到 Geopandas 和 Cartopy蛋白质表面制图映射

我一直在研究这个问题，但似乎无法解决这个问题。本质上，我想使用制图软件将任意 3D 表面映射到球体（我正在尝试使用 geopandas 和 cartopy）在下面链接的论文的 SI 中有一个代码示例，但他们使用了 Basemap，也是他们的代码片段对我来说也不是很好。

link to paper

我有一个顶点输入（附加 np.array[x,y,z,value]） 以下代码将表面映射到一个球体（地球的大小）

def cart2pol(coords):
    x = coords[0]
    y = coords[1]
    z = coords[2]
    rho = np.sqrt(x**2 + y**2 + z**2)
    theta = np.arctan2(y,x)
    phi = np.arccos(z/rho)
    return([rho,theta,phi])

def pol2cart(coords):
    rho = coords[0]
    theta = coords[1]
    phi = coords[2]
    x = rho * np.cos(theta) * np.sin(phi)
    y = rho * np.sin(theta) * np.sin(phi)
    z = rho * np.cos(phi)
    return([np.round(x,decimals=3),z])

def cart2geo(coords):
    x = coords[0]
    y = coords[1]
    z = coords[2]
    rho = np.sqrt(x**2 + y**2 + z**2)
    lat = np.rad2deg(np.arcsin(z / rho))
    lon = np.rad2deg(np.arctan2(y,x))
    return(lon,lat)

sph_coords = []
for vert in vertices:
    sph_coords.append(cart2pol(vert))
sphc_pol = np.array(sph_coords)
sphc_pol[:,0] = 6378100

sph_coords = []
for vert in sphc_pol:
    sph_coords.append(pol2cart(vert))
sphc_cart = np.array(sph_coords)

这段代码会给我我想要映射到的经度和纬度（在正射投影中）

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
import cartopy.crs as ccrs               
from pyproj import CRS
from scipy.interpolate import griddata
from matplotlib.colors import LinearSegmentedColormap 
from matplotlib import pyplot,transforms

crs = CRS("EPSG:4326")

lons = np.arange(-180,181,1)
lats = np.arange(-90,91,1)
x,y = np.meshgrid(lons,lats)
positions = np.vstack([x.ravel(),y.ravel()])
geo = gpd.points_from_xy(positions[0],positions[1],crs=crs)
geogs = gpd.GeoSeries(geo)
geogs.crs = crs

mol3 = ccrs.Orthographic()
crs3 = mol3.proj4_init

geogs2 = geogs.to_crs(crs3)

最终结果是插入一个可以映射到“地球”的表面，然后访问 cartopy 的投影。

surf = griddata(spherical_points,values,(lons[None,:],lats[:,None]),method = 'cubic')

sph_coords = []
for vert in vertices:
    sph_coords.append(cart2geo(vert))
sphc = np.array(sph_coords)
surf = griddata(sphc,df['value'],method = 'cubic')

colors = ['darkred','red','white','lightblue','blue']
cm = LinearSegmentedColormap.from_list('clist',colors)
base = ccrs.PlateCarree()

fig = plt.figure(figsize=(11,8.5))
levels = [-4,-3,-0.6,0.6,3,4]
ax=plt.axes(projection=ccrs.Mollweide())
ax.set_global()
ax.set_extent([-100,100,-50,50])
ax.gridlines(draw_labels=False,dms=True,x_inline=False,y_inline=False)
ax.contourf(lons,lats,surf,cmap=cm,levels=levels,extend='both',transform=base)

但是这个最终结果的几何值都搞砸了，并没有真正与原始球面投影对齐（上图）

我知道这是一篇很长的文章。但也许有一些制图大师可以轻松完成。我试图对我的购物车进行故障排除 -> lat long 转换，但似乎无处可去。我附上了我的顶点文件。谢谢！

numpy array