使用 Colab 笔记本可视化地理空间分析数据

使用集合让一切井井有条 根据您的偏好保存内容并对其进行分类。

创建 Colab 笔记本

本教程将构建一个 Colab 笔记本，以可视化地理空间分析数据。您可以点击 GitHub 版教程 Colab 中的 BigQuery 地理空间可视化顶部的链接，在 Colab、Colab Enterprise 或 BigQuery Studio 中打开预构建版本的笔记本。

1. 打开 Colab。

   打开 Colab

2. 在打开笔记本对话框中，点击新建笔记本。

3. 点击 Untitled0.ipynb，然后将笔记本的名称更改为 bigquery-geo.ipynb。

4. 依次选择文件 > 保存。

使用 Google Cloud 和 Google 地图进行身份验证

本教程将查询 BigQuery 数据集并使用 Google Maps JavaScript API。如需使用这些资源，您需要使用 Google Cloud 和 Google 地图 API 对 Colab 运行时进行身份验证。

使用 Google Cloud进行身份验证

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需使用您的项目进行身份验证，请输入以下代码：

   # REQUIRED: Authenticate with your project.
   GCP_PROJECT_ID = "PROJECT_ID"  #@param {type:"string"}
   
   from google.colab import auth
   from google.colab import userdata
   
   auth.authenticate_user(project_id=GCP_PROJECT_ID)

   将 PROJECT_ID 替换为您的项目 ID。

3. 点击 play_circle_filled 运行单元。

4. 如果您同意，请在系统提示时点击允许，以向 Colab 授予对您的凭据的访问权限。

5. 在使用 Google 账号登录页面上，选择您的账号。

6. 在登录第三方撰写的笔记编码页面上，点击继续。

7. 在选择第三方撰写的记事代码可以访问的内容中，点击全选，然后点击继续。

   完成授权流程后，您的 Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

可选：使用 Google 地图进行身份验证

如果您将 Google Maps Platform 用作底图的地图提供程序，则必须提供 Google Maps Platform API 密钥。该笔记本会从您的 Colab 密钥中检索密钥。

只有在您使用 Google 地图 API 时，才需要执行此步骤。如果您未使用 Google Maps Platform 进行身份验证，pydeck 将改用 carto 映射。

1. 按照 Google 地图文档中使用 API 密钥页面上的说明获取 Google Maps API 密钥。

2. 切换到您的 Colab 笔记本，然后点击 vpn_key Secrets（密钥）。

3. 点击添加新 Secret。

4. 在名称中，输入 GMP_API_KEY。

5. 在值中，输入您之前生成的 Maps API 密钥值。

6. 关闭密钥面板。

7. 如需插入代码单元，请点击 add Code（代码）。

8. 如需使用 Maps API 进行身份验证，请输入以下代码：

   # Authenticate with the Google Maps JavaScript API.
   GMP_API_SECRET_KEY_NAME = "GMP_API_KEY" #@param {type:"string"}
   
   if GMP_API_SECRET_KEY_NAME:
     GMP_API_KEY = userdata.get(GMP_API_SECRET_KEY_NAME) if GMP_API_SECRET_KEY_NAME else None
   else:
     GMP_API_KEY = None

9. 如果您同意，请在系统提示时点击授予访问权限，以向该记事本授予对您的密钥的访问权限。

10. 点击 play_circle_filled 运行单元。

    完成授权流程后，您的 Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

安装 Python 软件包并导入数据科学库

除了 colabtools (google.colab) Python 模块之外，本教程还使用了其他一些 Python 软件包和数据科学库。

在本部分中，您将安装 pydeck 和 h3 软件包。pydeck 在 Python 中提供高规模空间渲染，由 deck.gl 提供支持。h3-py 以 Python 语言提供 Uber 的 H3 六边形分层地理空间索引系统。

然后，导入 h3 和 pydeck 库以及以下 Python 地理空间库：

• geopandas，用于扩展 pandas 使用的类型，以允许对几何类型执行空间操作。
• shapely，用于操作和分析各个平面几何图形对象。
• branca 用于生成 HTML 和 JavaScript 颜色映射。
• geemap.deck（用于使用 pydeck 和 earthengine-api 进行可视化）。

导入库后，您可以为 Colab 中的 pandas DataFrame 启用交互式表格。

安装 pydeck 和 h3 软件包

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需安装 pydeck 和 h3 软件包，请输入以下代码：

   # Install pydeck and h3.
   !pip install pydeck>=0.9 h3>=4.2

3. 点击 play_circle_filled 运行单元。

   完成安装后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

导入 Python 库

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需导入 Python 库，请输入以下代码：

   # Import data science libraries.
   import branca
   import geemap.deck as gmdk
   import h3
   import pydeck as pdk
   import geopandas as gpd
   import shapely

3. 点击 play_circle_filled 运行单元。

   运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

为 pandas DataFrame 启用交互式表格

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需启用 pandas DataFrame，请输入以下代码：

   # Enable displaying pandas data frames as interactive tables by default.
   from google.colab import data_table
   data_table.enable_dataframe_formatter()

3. 点击 play_circle_filled 运行单元。

   运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

创建共享例程

在本部分中，您将创建两个共享例程：一个用于在底图上渲染图层的共享例程，以及一个用于将 pandas DataFrame 转换为 geopandas DataFrame 的共享例程。

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需创建用于在地图上渲染图层的共享例程，请输入以下代码：

   # Set Google Maps as the base map provider.
   MAP_PROVIDER_GOOGLE = pdk.bindings.base_map_provider.BaseMapProvider.GOOGLE_MAPS.value
   
   # Shared routine for rendering layers on a map using geemap.deck.
   def display_pydeck_map(layers, view_state, **kwargs):
     deck_kwargs = kwargs.copy()
   
     # Use Google Maps as the base map only if the API key is provided.
     if GMP_API_KEY:
       deck_kwargs.update({
         "map_provider": MAP_PROVIDER_GOOGLE,
         "map_style": pdk.bindings.map_styles.GOOGLE_ROAD,
         "api_keys": {MAP_PROVIDER_GOOGLE: GMP_API_KEY},
       })
   
     m = gmdk.Map(initial_view_state=view_state, ee_initialize=False, **deck_kwargs)
   
     for layer in layers:
       m.add_layer(layer)
     return m

3. 点击 play_circle_filled 运行单元。

   运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

4. 如需插入代码单元，请点击 add Code（代码）。

5. 如需创建用于将 pandas DataFrame 转换为 geopandas DataFrame 的共享例程，请输入以下代码：

   # Shared routine for converting a pandas dataframe to a GeoPandas dataframe.
   def pandas_to_geopandas(df, geometry_column='geometry'):
   
     # Use shapely library to parse WKT strings into shapely Geometry based
     # objects
     df[geometry_column] = df[geometry_column].dropna().apply(shapely.from_wkt)
   
     # Convert to a geopandas Dataframe
     return gpd.GeoDataFrame(df, geometry=geometry_column, crs='EPSG:4326')
     

6. 点击 play_circle_filled 运行单元。

   运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

创建散点图

在本部分中，您将从 bigquery-public-data.san_francisco_bikeshare.bikeshare_station_info 表中检索数据，以创建旧金山 Ford GoBike Share 公共数据集中所有共享单车站点的散点图。散点图是使用 deck.gl 框架中的图层和散点图图层创建的。

当您需要查看部分数据点（也称为抽查）时，散点图会很有用。

以下示例演示了如何使用图层和散点图图层将各个点渲染为圆形。如需渲染这些点，您需要从自行车共享数据集中的 station_geom 列中提取经度和纬度作为 x 坐标和 y 坐标。

通过将数据转换为 geopandas DataFrame 来提取坐标，该操作会将 station_geom 列中的每个项转换为 shapely 对象。shapely 提供了 Python 方法和属性，可让您提取 x 和 y 等组件。

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需查询旧金山 Ford GoBike 共享单车公共数据集，请输入以下代码。以下代码使用 %%bigquery 魔法函数运行查询，并以 DataFrame 的形式返回结果：

   # Query the station ID, station name, station short name, and station
   # geometry from the bike share dataset.
   # NOTE: In this tutorial, the denormalized 'lat' and 'lon' columns are
   # ignored. They are decomposed components of the geometry.
   %%bigquery df_sanfrancisco_bike_stations --project {GCP_PROJECT_ID}
   
   SELECT
     station_id,
     name,
     short_name,
     station_geom
   FROM
     `bigquery-public-data.san_francisco_bikeshare.bikeshare_station_info`

3. 点击 play_circle_filled 运行单元。

   输出类似于以下内容：

   Job ID 12345-1234-5678-1234-123456789 successfully executed: 100%

4. 如需插入代码单元，请点击 add Code（代码）。

5. 如需将数据转换为 geopandas DataFrame，请输入以下代码：

   # Convert the data to a geopandas.GeoDataFrame.
   gdf_sf_bikestations = pandas_to_geopandas(df_sanfrancisco_bike_stations, geometry_column='station_geom')
   gdf_sf_bikestations.info()

6. 点击 play_circle_filled 运行单元。

   输出应如下所示：

   <class 'geopandas.geodataframe.GeoDataFrame'>
   RangeIndex: 472 entries, 0 to 471
   Data columns (total 4 columns):
   #   Column        Non-Null Count  Dtype
   ---  ------        --------------  -----
   0   station_id    472 non-null    object
   1   name          472 non-null    object
   2   short_name    472 non-null    object
   3   station_geom  472 non-null    geometry
   dtypes: geometry(1), object(3)
   memory usage: 14.9+ KB
   

7. 如需插入代码单元，请点击 add Code（代码）。

8. 如需预览 DataFrame 的前五行，请输入以下代码：

   # Preview the first five rows
   gdf_sf_bikestations.head()

9. 点击 play_circle_filled 运行单元。

   输出类似于以下内容：

   

10. 如需插入代码单元，请点击 add Code（代码）。

11. 如需从 station_geom 列中提取经度和纬度值，请输入以下代码：

    # Extract the longitude (x) and latitude (y) from station_geom.
    gdf_sf_bikestations["longitude"] = gdf_sf_bikestations["station_geom"].x
    gdf_sf_bikestations["latitude"] = gdf_sf_bikestations["station_geom"].y

12. 点击 play_circle_filled 运行单元。

    运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

13. 如需插入代码单元，请点击 add Code（代码）。

14. 如需根据您之前提取的经纬度值渲染共享自行车站的散点图，请输入以下代码：

    # Render a scatter plot using pydeck with the extracted longitude and
    # latitude columns in the gdf_sf_bikestations geopandas.GeoDataFrame.
    scatterplot_layer = pdk.Layer(
      "ScatterplotLayer",
      id="bike_stations_scatterplot",
      data=gdf_sf_bikestations,
      get_position=['longitude', 'latitude'],
      get_radius=100,
      get_fill_color=[255, 0, 0, 140],  # Adjust color as desired
      pickable=True,
    )
    
    view_state = pdk.ViewState(latitude=37.77613, longitude=-122.42284, zoom=12)
    display_pydeck_map([scatterplot_layer], view_state)

15. 点击 play_circle_filled 运行单元。

    输出类似于以下内容：

    

直观呈现多边形

借助地理空间分析，您可以使用 GEOGRAPHY 数据类型和 GoogleSQL 地理位置函数来分析和直观呈现 BigQuery 中的地理空间数据。

在地理空间分析中，GEOGRAPHY 数据类型是一组点、线串和多边形，表示为一个点集或地球表面的部分区域。GEOGRAPHY 类型可以包含以下对象：

• 积分
• 行数
• 多边形
• 多边形集

如需查看所有受支持对象的列表，请参阅 GEOGRAPHY 类型文档。

如果您获得了地理空间数据，但不知道预期的形状，可以可视化数据来发现形状。您可以将地理数据转换为 GeoJSON 格式，以直观呈现形状。然后，您可以使用 deck.gl 框架中的 GeoJSON 图层直观呈现 GeoJSON 数据。

在本部分中，您将查询“San Francisco Neighborhoods”数据集中的地理数据，然后直观呈现多边形。

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需从 San Francisco Neighborhoods 数据集中的 bigquery-public-data.san_francisco_neighborhoods.boundaries 表中查询地理数据，请输入以下代码。以下代码使用 %%bigquery 魔法函数运行查询，并以 DataFrame 的形式返回结果：

   # Query the neighborhood name and geometry from the San Francisco
   # neighborhoods dataset.
   %%bigquery df_sanfrancisco_neighborhoods --project {GCP_PROJECT_ID}
   
   SELECT
     neighborhood,
     neighborhood_geom AS geometry
   FROM
     `bigquery-public-data.san_francisco_neighborhoods.boundaries`

3. 点击 play_circle_filled 运行单元。

   输出类似于以下内容：

   Job ID 12345-1234-5678-1234-123456789 successfully executed: 100%

4. 如需插入代码单元，请点击 add Code（代码）。

5. 如需将结果转换为 geopandas 格式，请输入以下代码：

   # Convert the query results to geopandas format.
   gdf_sanfrancisco_neighborhoods = pandas_to_geopandas(df_sanfrancisco_neighborhoods)
   gdf_sanfrancisco_neighborhoods.info()

6. 点击 play_circle_filled 运行单元。

   结果应如下所示：

   <class 'geopandas.geodataframe.GeoDataFrame'>
   RangeIndex: 117 entries, 0 to 116
   Data columns (total 2 columns):
   #   Column        Non-Null Count  Dtype
   ---  ------        --------------  -----
   0   neighborhood  117 non-null    object
   1   geometry      117 non-null    geometry
   dtypes: geometry(1), object(1)
   memory usage: 2.0+ KB
   

7. 如需预览 DataFrame 的第一行，请输入以下代码：

   # Preview the first row
   gdf_sanfrancisco_neighborhoods.head(1)

8. 点击 play_circle_filled 运行单元。

   输出类似于以下内容：

   

   请注意，结果中的数据是多边形。

9. 如需插入代码单元，请点击 add Code（代码）。

10. 如需直观呈现多边形，请输入以下代码。pydeck 用于将几何图形列中的每个 shapely 对象实例转换为 GeoJSON 格式：

    # Visualize the polygons.
    geojson_layer = pdk.Layer(
        'GeoJsonLayer',
        id="sf_neighborhoods",
        data=gdf_sanfrancisco_neighborhoods,
        get_line_color=[127, 0, 127, 255],
        get_fill_color=[60, 60, 60, 50],
        get_line_width=100,
        pickable=True,
        stroked=True,
        filled=True,
      )
    view_state = pdk.ViewState(latitude=37.77613, longitude=-122.42284, zoom=12)
    display_pydeck_map([geojson_layer], view_state)

11. 点击 play_circle_filled 运行单元。

    输出类似于以下内容：

    

创建分级统计图

如果您要探索的数据包含难以转换为 GeoJSON 格式的多边形，则可以改用 deck.gl 框架中的多边形图层。多边形图层可以处理特定类型的输入数据，例如点数组。

在本部分中，您将使用多边形图层渲染点数组，并使用结果渲染分级统计图。该 Choropleth 地图通过将“旧金山街区”数据集与“旧金山 Ford GoBike 共享单车”数据集的数据联接起来，按街区显示共享单车站的密度。

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需汇总和统计每个社区的充电站数量，并创建包含点数组的 polygon 列，请输入以下代码：

   # Aggregate and count the number of stations per neighborhood.
   gdf_count_stations = gdf_sanfrancisco_neighborhoods.sjoin(gdf_sf_bikestations, how='left', predicate='contains')
   gdf_count_stations = gdf_count_stations.groupby(by='neighborhood')['station_id'].count().rename('num_stations')
   gdf_stations_x_neighborhood = gdf_sanfrancisco_neighborhoods.join(gdf_count_stations, on='neighborhood', how='inner')
   
   # To simulate non-GeoJSON input data, create a polygon column that contains
   # an array of points by using the pandas.Series.map method.
   gdf_stations_x_neighborhood['polygon'] = gdf_stations_x_neighborhood['geometry'].map(lambda g: list(g.exterior.coords))

3. 点击 play_circle_filled 运行单元。

   运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

4. 如需插入代码单元，请点击 add Code（代码）。

5. 如需为每个多边形添加 fill_color 列，请输入以下代码：

   # Create a color map gradient using the branch library, and add a fill_color
   # column for each of the polygons.
   colormap = branca.colormap.LinearColormap(
     colors=["lightblue", "darkred"],
     vmin=0,
     vmax=gdf_stations_x_neighborhood['num_stations'].max(),
   )
   gdf_stations_x_neighborhood['fill_color'] = gdf_stations_x_neighborhood['num_stations'] \
     .map(lambda c: list(colormap.rgba_bytes_tuple(c)[:3]) + [0.7 * 255])   # force opacity of 0.7

6. 点击 play_circle_filled 运行单元。

   运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

7. 如需插入代码单元，请点击 add Code（代码）。

8. 如需渲染多边形图层，请输入以下代码：

   # Render the polygon layer.
   polygon_layer = pdk.Layer(
     'PolygonLayer',
     id="bike_stations_choropleth",
     data=gdf_stations_x_neighborhood,
     get_polygon='polygon',
     get_fill_color='fill_color',
     get_line_color=[0, 0, 0, 255],
     get_line_width=50,
     pickable=True,
     stroked=True,
     filled=True,
   )
   view_state = pdk.ViewState(latitude=37.77613, longitude=-122.42284, zoom=12)
   display_pydeck_map([polygon_layer], view_state)

9. 点击 play_circle_filled 运行单元。

   输出类似于以下内容：

   

创建热图

如果您有已知的有意义的边界，则可以使用 Choropleth 图。如果您的数据没有已知的有意义的边界，您可以使用热图图层来渲染其连续密度。

在以下示例中，您将查询旧金山警察局 (SFPD) 报告数据集中 bigquery-public-data.san_francisco_sfpd_incidents.sfpd_incidents 表中的数据。这些数据用于直观呈现 2015 年事故的分布情况。

对于热图，建议您在渲染之前对数据进行量化和汇总。在此示例中，系统使用 Carto H3 空间索引对数据进行量化和汇总。热图是使用 deck.gl 框架中的热图图层创建的。

在此示例中，使用 h3 Python 库对入射点进行汇总，将其汇总为六边形。h3.latlng_to_cell 函数用于将事故的位置（纬度和经度）映射到 H3 单元格索引。H3 分辨率为 9 可为热图提供足够的汇总六边形。h3.cell_to_latlng 函数用于确定每个六边形的中心。

1. 如需插入代码单元，请点击 add Code（代码）。

2. 如需查询旧金山警察局 (SFPD) 报告数据集中的数据，请输入以下代码。以下代码使用 %%bigquery 魔法函数运行查询，并在 DataFrame 中返回结果：

   # Query the incident key and location  data from the SFPD reports dataset.
   %%bigquery df_sanfrancisco_incidents_2015  --project {GCP_PROJECT_ID}
   
   SELECT
     unique_key,
     location
   FROM (
     SELECT
       unique_key,
       location, # WKT string
       EXTRACT(YEAR FROM timestamp) AS year,
     FROM `bigquery-public-data.san_francisco_sfpd_incidents.sfpd_incidents` incidents
   )
   WHERE year = 2015

3. 点击 play_circle_filled 运行单元。

   输出类似于以下内容：

   Job ID 12345-1234-5678-1234-123456789 successfully executed: 100%

4. 如需插入代码单元，请点击 add Code（代码）。

5. 如需将结果转换为 geopandas DataFrame，请输入以下代码：

   # Convert the results into a geopandas.GeoDataFrame.
   gdf_incidents = pandas_to_geopandas(df_sanfrancisco_incidents_2015, geometry_column='location')

6. 点击 play_circle_filled 运行单元。

   运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

7. 如需插入代码单元，请点击 add Code（代码）。

8. 如需计算每个事故的纬度和经度对应的单元格，汇总每个单元格的事故数据，构建 geopandas DataFrame，并为热图图层添加每个六边形的中心，请输入以下代码：

   # Compute the cell for each incident's latitude and longitude.
   H3_RESOLUTION = 9
   gdf_incidents['h3_cell'] = df_sanfrancisco_incidents_2015['location'].apply(
       lambda location: h3.latlng_to_cell(location.y, location.x, H3_RESOLUTION)
   )
   
   # Aggregate the incidents for each hexagon cell.
   count_incidents = gdf_incidents.groupby(by='h3_cell')['unique_key'].count().rename('num_incidents')
   
   # Construct a new geopandas.GeoDataFrame with the aggregate results.
   # Add the center of each hexagon for the HeatmapLayer to render.
   gdf_incidents_x_cell = gpd.GeoDataFrame(data=count_incidents).reset_index()
   gdf_incidents_x_cell['h3_center'] = gdf_incidents_x_cell['h3_cell'].apply(h3.cell_to_latlng)
   gdf_incidents_x_cell.info()

9. 点击 play_circle_filled 运行单元。

   输出类似于以下内容：

   <class 'geopandas.geodataframe.GeoDataFrame'>
   RangeIndex: 969 entries, 0 to 968
   Data columns (total 3 columns):
   #   Column         Non-Null Count  Dtype
   --  ------         --------------  -----
   0   h3_cell        969 non-null    object
   1   num_incidents  969 non-null    Int64
   2   h3_center      969 non-null    object
   dtypes: Int64(1), object(2)
   memory usage: 23.8+ KB
   

10. 如需插入代码单元，请点击 add Code（代码）。

11. 如需预览 DataFrame 的前五行，请输入以下代码：

    # Preview the first five rows.
    gdf_incidents_x_cell.head()

12. 点击 play_circle_filled 运行单元。

    输出类似于以下内容：

    

13. 如需插入代码单元，请点击 add Code（代码）。

14. 如需将数据转换为 HeatmapLayer 可使用的 JSON 格式，请输入以下代码：

    # Convert to a JSON format recognized by the HeatmapLayer.
    def _make_heatmap_datum(row) -> dict:
      return {
          "latitude": row['h3_center'][0],
          "longitude": row['h3_center'][1],
          "weight": float(row['num_incidents']),
      }
    
    heatmap_data = gdf_incidents_x_cell.apply(_make_heatmap_datum, axis='columns').values.tolist()

15. 点击 play_circle_filled 运行单元。

    运行代码后，Colab 笔记本中不会生成任何输出。单元格旁边的对勾标记表示代码已成功运行。

16. 如需插入代码单元，请点击 add Code（代码）。

17. 如需渲染热图，请输入以下代码：

    # Render the heatmap.
    heatmap_layer = pdk.Layer(
      "HeatmapLayer",
      id="sfpd_heatmap",
      data=heatmap_data,
      get_position=['longitude', 'latitude'],
      get_weight='weight',
      opacity=0.7,
      radius_pixels=99,  # this limitation can introduce artifacts (see above)
      aggregation='MEAN',
    )
    view_state = pdk.ViewState(latitude=37.77613, longitude=-122.42284, zoom=12)
    display_pydeck_map([heatmap_layer], view_state)

18. 点击 play_circle_filled 运行单元。

    输出类似于以下内容：