使用 Python 库 Manim 创建动态地理空间可视化

https://readmedium.com/zh/create-geo-spatial-visualization-using-manim-2d179b2c21b9

使用 Manim 创建动态地理空间可视化

开始之前,我们需要安装两个主要软件包。以下是软件包名称和安装文件:

Manim

Geopandas

现在,让我们收集一些数据来进行可视化,我将以 2017 年美国各州人口数据为例,但你也可以使用任何数据。

import geopandas as gpd
import matplotlib.pyplot as plt
from manim import *

us_population_2017 = gpd.read_file("https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/USA_States_Generalized/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson")
us_population_2017.head()

img

让我们使用 Matplotlib 将数据可视化。当然,也可以使用 matplotlib 库制作动画,但我们将在另一篇文章中讨论这个问题。

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 7))
cax = fig.add_axes([0.93, 0.2, 0.02, 0.6])
us_population_2017.plot(ax=ax, 
                        column='POP_SQMI', 
                        legend=True, 
                        cax=cax,
                        vmax=500,
                        legend_kwds = {'label': 'per sqmi'}
                       )

custom_ticklabels = []
for i in cax.get_yticklabels():
    ticklabel = i.get_text()
    if ticklabel=='500':
        custom_ticklabels.append(f'>{ticklabel}')
    else:
        custom_ticklabels.append(ticklabel)
        
cax.set_yticklabels(custom_ticklabels)

ax.set_title('Population per Square Mile')
plt.show()

下面是静态地图的样子

img

从地理空间数据集创建 Manim 动画的步骤可分为 3 个主要步骤:

1) 将 GIS 几何数据转换为 Manim 对象;
2) 设置轴线并应用适当的比例;
3) 根据属性应用颜色;
4) 自定义动画。

将 GIS 几何数据转换为 Manim 对象:要将多边形 GIS 几何数据转换为 Manim 线条和多边形对象,以实现纵横图可视化,请按照以下步骤操作:

提取坐标:从 GIS 数据集中提取多边形的顶点坐标。这些坐标定义了每个多边形的边界。

创建直线和多边形对象使用提取的坐标创建马尼姆直线和多边形对象。每个多边形对应一个由其顶点定义的封闭形状。此外,如有必要,还可创建线条对象来表示多边形的边界。

设置坐标轴并应用适当的比例:确定要用于可视化的坐标系。这可以是笛卡尔坐标系、极坐标系或任何其他适合数据的几何投影。我们可以通过创建轴对象来实现这一点:使用 Manim 的轴类来创建基于所选坐标系的轴对象。为每个轴指定范围、标签和刻度线等参数

下面是一个使用美国边界的示例。为简单起见,我将阿拉斯加和夏威夷排除在外。

%%manim -qm -v WARNING Animate_population_density
# the above magic command is needed for running Manim in jupyter notebook

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

exploded_usa = world[world['iso_a3']=='USA'].explode() # to convert multipolygon to polygon
exploded_usa['part'] = [f'part_{i}' for i in range(len(exploded_usa))]
# exploded_usa.explore(column='part')

polygon = exploded_usa[exploded_usa['part']=='part_0'].geometry.values[0] # excluding Alaska and Hawai

#extract the coordinates
polygon_border_xy = np.array(polygon.boundary.coords)

class Animate_population_density(Scene):
   def construct(self):
        # set up Manim Axis
        axes = Axes(
            x_range=[min(polygon_border_xy[:,0])-2.5, max(polygon_border_xy[:,0])+2.5],
            y_range=[min(polygon_border_xy[:,1]), max(polygon_border_xy[:,1])],
            axis_config={"color": BLUE},
        )
        
        boundary_line = axes.plot_line_graph(polygon_border_xy[:,0], polygon_border_xy[:,1], 
                                       add_vertex_dots=False, 
                                       line_color=WHITE,
                                       stroke_width=5 )
#         bd_line.move_to(LEFT*3)
        self.play( Create(boundary_line, run_time=2))
        self.wait(1)

运行这段代码后,我们在 “media/videos “目录下创建了名为 “Animate_population_density.mp4 “的动画。

img

以下是为各州生成人口密度 Choropleth 动画的完整代码。

#Full code 
import geopandas as gpd
import matplotlib.pyplot as plt
from manim import *
from colormap import rgb2hex

us_population_2017 = gpd.read_file("https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/USA_States_Generalized/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson")
relevant_columns = ['STATE_NAME', 'POPULATION', 'POP_SQMI', 'geometry']
# us_population_2017[relevant_columns].head()

us_population_2017 = us_population_2017.explode()

us_population_2017 = us_population_2017[~us_population_2017['STATE_NAME'].isin(['Alaska', 'Hawaii'])]

us_population_2017 = us_population_2017.sort_values(by='POP_SQMI', ascending=False)

cmap = plt.get_cmap('viridis')
norm = plt.Normalize(us_population_2017['POP_SQMI'].min(), 500)


def get_line_coord(linestring):

    border_xy = np.array(linestring.boundary.coords)
    
    return border_xy



def get_line_and_area(axes, border_xy, pop):

    color = cmap(norm(pop))
    hex_color = rgb2hex(*color)

    line = axes.plot_line_graph(border_xy[:,0], border_xy[:,1], 
                                add_vertex_dots=False,
                                line_color=BLACK, 
                                stroke_width=1 )
    points = axes.coords_to_point(border_xy)
    area = Polygon(*points,fill_opacity=0.5, color=hex_color, stroke_width=1)
    dist_center = area.get_center()

    # line.to_edge(UR)
    # area.to_edge(UR)

    return line, area, dist_center


world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

exploded_usa = world[world['iso_a3']=='USA'].explode() # to convert multipolygon to polygon
exploded_usa['part'] = [f'part_{i}' for i in range(len(exploded_usa))]
# exploded_usa.explore(column='part')

polygon = exploded_usa[exploded_usa['part']=='part_0'].geometry.values[0] # excluding Alaska and Hawai

#extract the coordinates
polygon_border_xy = np.array(polygon.boundary.coords)

class Animate_population_density(Scene):
   def construct(self):
    
        self.camera.background_color = WHITE
        # set up Manim Axis
        axes = Axes(
            x_range=[min(polygon_border_xy[:,0])-2.5, max(polygon_border_xy[:,0])+2.5],
            y_range=[min(polygon_border_xy[:,1]), max(polygon_border_xy[:,1])],
            axis_config={"color": BLUE},
        )
        
        boundary_line = axes.plot_line_graph(polygon_border_xy[:,0], polygon_border_xy[:,1], 
                                       add_vertex_dots=False, 
                                       line_color=BLACK,
                                       stroke_width=5 )
#         bd_line.move_to(LEFT*3)
        self.play( Create(boundary_line, run_time=1))
        self.wait(0.5)

        
        
        for state in us_population_2017['STATE_NAME'].unique():
            state_df = us_population_2017[us_population_2017['STATE_NAME']==state]
            
            if len(state_df)>1:
                
                states_line_group = VGroup()
                state_area_group = VGroup()
                
                for row in state_df.itertuples():
                
                    geometry = row.geometry
                    density = row.POP_SQMI

                    state_line, state_area, center = get_line_and_area(axes, 
                                                                       get_line_coord(geometry), 
                                                                       density
                                                                      )

                    states_line_group.add(state_line)
                    state_area_group.add(state_area)





                state_label = Text(f"{state}: {density:.2f} per square mile ",
                                   font_size=25, color=BLACK)

                state_label.move_to(UP + 2.3 * UP)

                self.play( Create(states_line_group) , Write(state_label),  run_time=2)

                self.play( Create(state_area_group), run_time=2)

                self.play(FadeOut(state_label))
               
                
            else:
                geometry = state_df.geometry.values[0]
                density = state_df.POP_SQMI.values[0]

                state_line, state_area, center = get_line_and_area(axes, 
                                                                   get_line_coord(geometry), 
                                                                   density
                                                                  )

                
                state_label = Text(f"{state}: {density:.2f} per square mile ",
                                   font_size=25, color=BLACK)
                
                state_label.move_to(UP + 2.3 * UP)
                
                self.play( Create(state_line) , Write(state_label),  run_time=2)
                
                self.play( Create(state_area), run_time=2)
                
                self.play(FadeOut(state_label))

这段代码将生成这段视频:

https://youtu.be/lAWEr46F_Kw