我有一个超过500个点的列表,给出了纬度和经度。这些点表示环形山,我想绘制这些环形山的热图。例如,我希望有很多环形山的区域被认为是“热”的,而较少的环形山是“冷”的。我看过使用SciPy的KDE,也尝试过在Mathematica中使用ListSliceDensityPlot3D,但我一直无法创建一个足够的图形。
我将每个点从纬度/经度转换为笛卡尔x,y,z坐标,并将它们绘制在球体的表面上,但我不知道如何获取点列表并计算给定区域的密度,然后将其绘制在3D表面上。
我的想法是得到一个类似于this image of Ceres的图!
提前感谢,如果需要的话,请提问,如果我一开始没有张贴足够的信息,很抱歉。
发布于 2018-10-24 17:01:14
这是一种蛮力方法,但它在一定程度上是有效的。这将是有问题的,如果你使网格非常精细或有数千个环形山。如果容器大小足够小,表面上的距离和3D距离之间没有太大的差异,所以我采用了3D距离,因为它更容易计算,但有人可能想要改变这一点。
代码看起来像这样:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
def random_point( r=1 ):
ct = 2*np.random.rand() - 1
st = np.sqrt( 1 - ct**2 )
phi = 2* np.pi * np.random.rand()
x = r * st * np.cos( phi)
y = r * st * np.sin( phi)
z = r * ct
return np.array( [x, y, z ] )
def near( p, pntList, d0 ):
cnt=0
for pj in pntList:
dist=np.linalg.norm( p - pj )
if dist < d0:
cnt += 1 - dist/d0
return cnt
"""
https://stackoverflow.com/questions/22128909/plotting-the-temperature-distribution-on-a-sphere-with-python
"""
pointList = np.array([ random_point( 10.05 ) for i in range( 65 ) ] )
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1, projection='3d')
u = np.linspace( 0, 2 * np.pi, 120)
v = np.linspace( 0, np.pi, 60 )
# create the sphere surface
XX = 10 * np.outer( np.cos( u ), np.sin( v ) )
YY = 10 * np.outer( np.sin( u ), np.sin( v ) )
ZZ = 10 * np.outer( np.ones( np.size( u ) ), np.cos( v ) )
WW = XX.copy()
for i in range( len( XX ) ):
for j in range( len( XX[0] ) ):
x = XX[ i, j ]
y = YY[ i, j ]
z = ZZ[ i, j ]
WW[ i, j ] = near(n p.array( [x, y, z ] ), pointList, 3)
WW = WW / np.amax( WW )
myheatmap = WW
# ~ ax.scatter( *zip( *pointList ), color='#dd00dd' )
ax.plot_surface( XX, YY, ZZ, cstride=1, rstride=1, facecolors=cm.jet( myheatmap ) )
plt.show() 结果是这样的:

您还可以修改距离函数,以考虑到环形山的大小。
发布于 2018-11-15 01:42:04
这里有两个独立的问题:在球体上定义合适的bin大小和计数函数,以便您可以构建适当的颜色函数,然后在3D球体上绘制该颜色函数。我将为这两个问题提供一个数学解决方案。
1.构建示例数据
首先,这里有一些示例数据:
data = Map[
Apply@Function[{x, y, z},
{
Mod[ArcTan[x, y], 2 π],
ArcSin[z/Sqrt[x^2 + y^2 + z^2]]
}
],
Map[
#/Norm[#]&,
Select[
RandomReal[{-1, 1}, {200000, 3}],
And[
Norm[#] > 0.01,
Or[
Norm[# - {0, 0.3, 0.2}] < 0.6,
Norm[# - {-0.3, -0.15, -0.3}] < 0.3
]
] &
]
]
]我把它弄得有点凹凸不平,这样在绘图时间上就会有更多有趣的特性。
2.构建颜色函数
要在Mathematica中构建颜色函数,最干净的解决方案是使用HistogramList,但这需要进行修改,以考虑到高纬度的存储箱将具有不同的区域,因此需要调整密度。
尽管如此,内置的直方图构建工具还是相当不错的:
DensityHistogram[
data,
{5°}
, AspectRatio -> Automatic
, PlotRangePadding -> None
, ImageSize -> 700
]

您可以通过以下方式获取原始数据
{{ϕbins, θbins}, counts} = HistogramList[data, {15°}]为了方便起见,让我们定义
ϕcenters = 1/2 (Most[ϕbins] + Rest[ϕbins])
θcenters = 1/2 (Most[θbins] + Rest[θbins])使用以下命令计算bin面积:
SectorArea[ϕmin_, ϕmax_, θmin_, θmax_] = (Abs[ϕmax - ϕmin]/(4 π)) *
Integrate[Sin[θ], {θ, θmin, θmax}]然后,您可以将自己的颜色函数定义为
function[ϕ_, θ_] := With[{
iϕ = First[Nearest[ϕcenters -> Range[Length[ϕcenters]], ϕ]],
iθ = First[Nearest[θcenters -> Range[Length[θcenters]], θ]]
},
(N@counts[[iϕ, iθ]]/
SectorArea[ϕbins[[iϕ]], ϕbins[[iϕ + 1]], θbins[[iθ]], θbins[[iθ + 1]]])/max
]所以,下面是实际运行的函数:
texture = ListDensityPlot[
Flatten[
Table[
{
ϕcenters[[iϕ]],
θcenters[[iθ]],
function[ϕcenters[[iϕ]], θcenters[[iθ]]]
}
, {iϕ, Length[ϕbins] - 1}
, {iθ, Length[θbins] - 1}
], 1]
, InterpolationOrder -> 0
, AspectRatio -> Automatic
, ColorFunction -> ColorData["GreenBrownTerrain"]
, Frame -> None
, PlotRangePadding -> None
]

3.绘图
要在球体上绘制数据,我看到了两个主要选项:您可以创建一个曲面图,然后将其作为Texture包装在参数图周围,如
ParametricPlot3D[
{Cos[ϕ] Sin[θ], Sin[ϕ] Sin[θ],
Cos[θ]}
, {ϕ, 0, 2 π}, {θ, 0, π}
, Mesh -> None
, Lighting -> "Neutral"
, PlotStyle -> Directive[
Specularity[White, 30],
Texture[texture]
]
]

或者,您可以将其定义为同一参数图中的显式ColorFunction:
ParametricPlot3D[
{Cos[ϕ] Sin[θ], Sin[ϕ] Sin[θ],
Cos[θ]}
, {ϕ, 0, 2 π}, {θ, 0, π}
, ColorFunctionScaling -> False
, ColorFunction -> Function[{x, y, z, ϕ, θ},
ColorData["GreenBrownTerrain"][function[ϕ, θ]]
]
]

当然,以上所有功能都是非常模块化的,因此您可以根据自己的优势自由地进行混合和匹配。
https://stackoverflow.com/questions/52937846
复制相似问题