I have been working on this problem for sometime and just cant seem to get it right. Essentially I want to map an arbitrary 3D surface to a sphere using cartography software (I am trying to do it with geopandas and cartopy) There is a code example in the SI of the paper linked below, but they used Basemap, also their code snippet didn't work for me very well either.
I have a input of vertices (attached np.array[x,y,z,value]) The following code will map the surface to a sphere (size of the globe)
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), y, 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)
and this code will give me the longitudes and latitudes (in ortho graphic projection) that I want to map to
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)
The end result is to interpolate a surface that can be mapped to the "globe" and then access the projections of 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'], (lons[None, :], lats[:, None]), method = 'cubic')
colors = ['darkred', 'red', 'white', 'white', '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, 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)
But the geometry of the values of this end result are all screwed up and don't really align with the original spherical projection (middle picture above)
I know this is a long post. But maybe there are some cartography gurus out there for which this is a snap. I have tried to troubleshoot my cart -> lat long conversions but cant seem to get anywhere. I have attached my vertices file. Thanks!
