首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >在给定纬度和经度列表的情况下,如何在球体上绘制热图?

在给定纬度和经度列表的情况下,如何在球体上绘制热图?
EN

Stack Overflow用户
提问于 2018-10-23 05:14:07
回答 2查看 1.9K关注 0票数 8

我有一个超过500个点的列表,给出了纬度和经度。这些点表示环形山,我想绘制这些环形山的热图。例如,我希望有很多环形山的区域被认为是“热”的,而较少的环形山是“冷”的。我看过使用SciPy的KDE,也尝试过在Mathematica中使用ListSliceDensityPlot3D,但我一直无法创建一个足够的图形。

我将每个点从纬度/经度转换为笛卡尔x,y,z坐标,并将它们绘制在球体的表面上,但我不知道如何获取点列表并计算给定区域的密度,然后将其绘制在3D表面上。

我的想法是得到一个类似于this image of Ceres的图!

提前感谢,如果需要的话,请提问,如果我一开始没有张贴足够的信息,很抱歉。

EN

回答 2

Stack Overflow用户

发布于 2018-10-24 17:01:14

这是一种蛮力方法,但它在一定程度上是有效的。这将是有问题的,如果你使网格非常精细或有数千个环形山。如果容器大小足够小,表面上的距离和3D距离之间没有太大的差异,所以我采用了3D距离,因为它更容易计算,但有人可能想要改变这一点。

代码看起来像这样:

代码语言:javascript
复制
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() 

结果是这样的:

您还可以修改距离函数,以考虑到环形山的大小。

票数 6
EN

Stack Overflow用户

发布于 2018-11-15 01:42:04

这里有两个独立的问题:在球体上定义合适的bin大小和计数函数,以便您可以构建适当的颜色函数,然后在3D球体上绘制该颜色函数。我将为这两个问题提供一个数学解决方案。

1.构建示例数据

首先,这里有一些示例数据:

代码语言:javascript
复制
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,但这需要进行修改,以考虑到高纬度的存储箱将具有不同的区域,因此需要调整密度。

尽管如此,内置的直方图构建工具还是相当不错的:

代码语言:javascript
复制
DensityHistogram[
 data,
 {5°}
 , AspectRatio -> Automatic
 , PlotRangePadding -> None
 , ImageSize -> 700
 ]

您可以通过以下方式获取原始数据

代码语言:javascript
复制
{{ϕbins, θbins}, counts} =  HistogramList[data, {15°}]

为了方便起见,让我们定义

代码语言:javascript
复制
ϕcenters = 1/2 (Most[ϕbins] + Rest[ϕbins])
θcenters = 1/2 (Most[θbins] + Rest[θbins])

使用以下命令计算bin面积:

代码语言:javascript
复制
SectorArea[ϕmin_, ϕmax_, θmin_, θmax_] = (Abs[ϕmax - ϕmin]/(4 π)) * 
                                         Integrate[Sin[θ], {θ, θmin, θmax}]

然后,您可以将自己的颜色函数定义为

代码语言:javascript
复制
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
  ]

所以,下面是实际运行的函数:

代码语言:javascript
复制
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包装在参数图周围,如

代码语言:javascript
复制
ParametricPlot3D[
 {Cos[ϕ] Sin[θ], Sin[ϕ] Sin[θ], 
  Cos[θ]}
 , {ϕ, 0, 2 π}, {θ, 0, π}
 , Mesh -> None
 , Lighting -> "Neutral"
 , PlotStyle -> Directive[
   Specularity[White, 30],
   Texture[texture]
   ]
 ]

或者,您可以将其定义为同一参数图中的显式ColorFunction

代码语言:javascript
复制
ParametricPlot3D[
 {Cos[ϕ] Sin[θ], Sin[ϕ] Sin[θ], 
  Cos[θ]}
 , {ϕ, 0, 2 π}, {θ, 0, π}
 , ColorFunctionScaling -> False
 , ColorFunction -> Function[{x, y, z, ϕ, θ},
   ColorData["GreenBrownTerrain"][function[ϕ, θ]]
   ]
 ]

当然,以上所有功能都是非常模块化的,因此您可以根据自己的优势自由地进行混合和匹配。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/52937846

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档