使用 Cartopy 在 Python 中进行地理可视化

最新推荐文章于 2024-10-22 22:17:57 发布

nimabcsd

最新推荐文章于 2024-10-22 22:17:57 发布

阅读量1.1k

点赞数

文章标签： python 大数据

本文链接：https://blog.csdn.net/nimabcsd/article/details/124382707

版权

在这里插入图片描述
Cartopy 是一个制图 Python 库，因此它的命名约定是为地理数据操作和可视化中的应用程序开发的。它是 Basemap Toolkit 的继承者，后者是以前用于地理可视化的 Python 库（可以在我们的网站上找到几个 Basemap 教程：来自 CSV 文件的地理映射、卫星图像分析、GOES-16 纬度/经度投影算法）。Cartopy 可用于在真实地图上绘制卫星数据、可视化城市和国家边界、基于地理定位跟踪和预测运动，以及与地理编码数据系统相关的一系列其他应用。在本教程中，Anaconda 3 将用于安装 Cartopy 和相关的地理库。作为库和地理可视化的介绍，将进行一些简单的测试，以确保 Cartopy 库已成功安装并正常工作。在后续教程中：将使用 shapefile 作为边界，绘制真实的城市街道，并分析卫星数据。

Anaconda 3 被选为 Python 平台是因为它对 Python 库和包的强大处理。它是免费的、开源的，并在安装新软件包时动态管理软件包。这对于需要广泛先决条件的复杂库（例如 Cartopy）尤为重要。另一个建议是**Mac 或 Linux与 Cartopy 一起使用。原因是 -由于环境变量、编程导航和管理权限的问题，一些 Python库在 Windows 平台上遇到了困难。**因此，强烈建议使用 Linux 或 Mac 操作系统。

安装 Anaconda 3

以下安装过程适用于从 Linux 操作系统（此处使用 Ubuntu 19.x）基于终端安装 Anaconda 3：

1. 从（Python 3.x 版本）下载 Anaconda 3：https 😕/www.anaconda.com/distribution/#download-section

2. 下载完成后，导航到下载的文件夹并键入：

user@user1:~ $ cd ./Downloads
user@user1:~/Downloads $ bash Anaconda3-2020.02-Linux-x86_64.sh

3.接下来，安装程序会要求用户同意条款，点击进入并在提示时输入yes。有关完整的下载过程，请转到 Anaconda 页面并按照他们的说明进行操作：https 😕/docs.anaconda.com/anaconda/install/linux/ 。此处将继续使用默认值（同意所有条款和默认初始化）。

4.一旦安装到达“感谢您安装Anaconda3” - 关闭终端

5. 重新打开终端后，键入以下内容：

user@user1:~ $ source ~/.bashrc

6. 接下来，通过键入以下命令查看 Anaconda 3 中的当前环境，确保已安装 conda：

user@user1:~ $ conda 环境列表

此时，Anaconda 3 应该已经安装在本地系统上，上面应该打印出以下内容：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-G64i9Sh5-1650782775886)(https://images.squarespace-cdn.com/content/v1/59b037304c0dbfb092fbe894/1587753538439-IRRGA8YDXFCGKO1AXLLI/anaconda_active_base_environment.png)]

这意味着已经安装了基础 Anaconda 3 环境。

安装 Cartopy 和依赖项

我们现在需要做的是创建一个全新的环境，并根据 Cartopy 库和依赖项对其进行调整。Anaconda 团队有一个备忘单，在 Anaconda 中创建、更新和导航环境时非常有用，可以在此处找到。在那里，用户可以找到快速的终端命令来创建和更改环境。对于 Cartopy，第一个命令将涉及创建一个名为“cartopy_env”的 Cartopy 环境，该环境将包含所有必要的地理映射库：

user@user1:~ $ conda create -n cartopy_env cartopy --yes -c conda-forge

第一个命令安装环境“cartopy_env”并使用“conda-forge”作为安装通道。使用“conda-forge”允许用户以动态方式在 Anaconda 环境中维护库，而无需安装冲突或不兼容的库或版本。

创建新环境后，必须将其激活。在终端中键入并运行以下命令：

user@user1:~ $ conda activate cartopy_env

用户现在应该会看到一个更新的终端行，其中新环境处于活动状态，它应该类似于以下内容：

( cartopy_env )用户@user1:~ $

( cartopy_env ) user@user1:~ $ conda info --envs

如果星号位于新环境“cartopy_env”旁边，则安装过程已完成并且已安装正确的环境。接下来，可以安装 Jupyter Lab，让 Python 编程更加通用。Jupyter Lab 是一个通过 Web 界面工作的 GUI，并连接到 Jupyter Notebook 系列，使其成为实时可视化和编码的绝佳工具。假设“cartopy_env”处于活动状态，可以使用以下命令下载 Jupyter 框架：

( cartopy_env ) user@user1:~ $ conda install -c conda-forge jupyterlab

上述命令将再次使用“conda-forge”来安装 Jupyter Lab 和 Jupyter Notebook 依赖项。在安装使用环境所需的任何新库时，应使用此通用命令：

conda install -c conda-forge ______

这种形式将使所有图书馆保持在同一频道上，并防止任何冲突或不兼容。Jupyter Lab 下载完成后，通过在活动的“cartopy_env”下键入以下内容来启动它：

( cartopy_env ) user@user1:~ $ jupyter lab

Jupyter Lab 现在应该在用户的本地 Web 浏览器中打开，并显示类似于以下内容：

至此，Anaconda 3 和 Cartopy 已准备好进行测试。可以通过输入“import cartopy.crs as ccrs”来进一步验证 Cartopy 的安装——如果脚本运行后没有错误，那么库和环境就得到了验证！下一节将简要介绍 Cartopy。

本节将使用 Jupyter Lab 来输出 Cartopy 生成的可视化。处理新软件时最简单的起点是尝试与已知条件相关的示例。Maker Portal 团队位于纽约市 (NYC)，因此，下面的示例使用世界贸易中心一号楼 (40.713, -74.0135) 的位置，使用开放街道地图 (OSM) 绘制纽约市周边地区的地图) 和 Cartopy：


```python
# Mapping New York City Open Street Map (OSM) with Cartopy
# This code uses a spoofing algorithm to avoid bounceback from OSM servers
#
import matplotlib.pyplot as plt 
import numpy as np
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt
import io
from urllib.request import urlopen, Request
from PIL import Image

def image_spoof(self, tile): # this function pretends not to be a Python script
    url = self._image_url(tile) # get the url of the street map API
    req = Request(url) # start request
    req.add_header('User-agent','Anaconda 3') # add user agent to request
    fh = urlopen(req) 
    im_data = io.BytesIO(fh.read()) # get image
    fh.close() # close url
    img = Image.open(im_data) # open image with PIL
    img = img.convert(self.desired_tile_form) # set image format
    return img, self.tileextent(tile), 'lower' # reformat for cartopy

cimgt.OSM.get_image = image_spoof # reformat web request for street map spoofing
osm_img = cimgt.OSM() # spoofed, downloaded street map

fig = plt.figure(figsize=(12,9)) # open matplotlib figure
ax1 = plt.axes(projection=osm_img.crs) # project using coordinate reference system (CRS) of street map
center_pt = [40.713, -74.0135] # lat/lon of One World Trade Center in NYC
zoom = 0.00075 # for zooming out of center point
extent = [center_pt[1]-(zoom*2.0),center_pt[1]+(zoom*2.0),center_pt[0]-zoom,center_pt[0]+zoom] # adjust to zoom
ax1.set_extent(extent) # set extents

scale = np.ceil(-np.sqrt(2)*np.log(np.divide(zoom,350.0))) # empirical solve for scale based on zoom
scale = (scale<20) and scale or 19 # scale cannot be larger than 19
ax1.add_image(osm_img, int(scale)) # add OSM with zoom specification
# NOTE: zoom specifications should be selected based on extent:
# -- 2     = coarse image, select for worldwide or continental scales
# -- 4-6   = medium coarseness, select for countries and larger states
# -- 6-10  = medium fineness, select for smaller states, regions, and cities
# -- 10-12 = fine image, select for city boundaries and zip codes
# -- 14+   = extremely fine image, select for roads, blocks, buildings
plt.show() # show the plot

Jupyter Lab Notebook 中的结果输出应与以下内容相同：

在这里插入图片描述
上面输出的地图是开放街道地图 (OSM) 可能的最高分辨率。如果我们从该点缩小一点，我们会看到不同的比例集。可以自定义比例，但随着比例的增加，绘制地图所需的时间也会增加。这就是添加自动缩放功能的原因。

让我们尝试另一种称为“QuadtreeTiles”的地图 - 这是对表面的照片般逼真的描绘。例如，如果我们在上面的代码中用“QuadtreeTiles”代替“OSM”，我们会得到以下世界贸易中心一号大楼的快照：

在这里插入图片描述
上面成功的城市可视化表明 Anaconda 3 正在运行，Cartopy 已成功安装，并且我们的测试准确且可重复。对世界贸易中心一号楼的谷歌快速搜索表明，上面的快照是准确的，并且在一定程度上反映了谷歌地图上的视觉效果（不包括多年来不一致的任何更新）。在下一节中，气象站的散点将绘制在美国地图上，以进一步探索 Cartopy 的功能。

与自动地面观测系统 (ASOS) 相关的气象站数据库可在以下链接中找到：

ftp://ftp.ncdc.noaa.gov/pub/data/ASOS_Station_Photos/asos-stations.txt

ASOS 网络是气象站的分布，每分钟更新温度、湿度、风向、降水和其他天气信息。ASOS 网络主要用于飞机和机场在起飞和着陆期间报告天气和安全情况。此外，ASOS 站已用于与卫星算法、数值天气预报和基于网络的天气应用相关的研究。ASOS .txt 文件包含分散在美国及其领土上的 915 个站点。从第 5 行开始，解析站。每个站有 14 列，可以从第 10 列和第 11 列中提取经纬度坐标。下载“asos-stations.txt”文件后，数据可以解析成 Python 并用作查找当地天气详细信息的工具。这将通过 Cartopy 和其他 Python 库进行探索。

下面的代码解析 ASOS .txt 文件并将站点映射到毗邻的美国（有时称为 CONUS）：

# Mapping ASOS Weather Stations Across the Contiguous United States
# This code uses a spoofing algorithm to avoid bounceback from OSM servers
# -- This code also parses lat/lon information from the ASOS station .txt
# -- file located here: ftp://ftp.ncdc.noaa.gov/pub/data/ASOS_Station_Photos/asos-stations.txt
#
import csv
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import io
from urllib.request import urlopen, Request
from PIL import Image

def image_spoof(self, tile): # this function pretends not to be a Python script
    url = self._image_url(tile) # get the url of the street map API
    req = Request(url) # start request
    req.add_header('User-agent','Anaconda 3') # add user agent to request
    fh = urlopen(req) 
    im_data = io.BytesIO(fh.read()) # get image
    fh.close() # close url
    img = Image.open(im_data) # open image with PIL
    img = img.convert(self.desired_tile_form) # set image format
    return img, self.tileextent(tile), 'lower' # reformat for cartopy

################################
# parsing the ASOS coordinates
################################
#
asos_data = []
with open('asos-stations.txt','r') as dat_file:
    reader = csv.reader(dat_file)
    for row in reader:
        asos_data.append(row)

row_delin = asos_data[3][0].split(' ')[:-1]
col_sizes = [len(ii) for ii in row_delin]

col_header = []; iter_ii = 0
for ii,jj in enumerate(col_sizes):
    col_header.append(asos_data[2][0][iter_ii:iter_ii+col_sizes[ii]].replace(' ',''))
    iter_ii+=col_sizes[ii]+1
    
call,names,lats,lons,elevs = [],[],[],[],[]
for row in asos_data[4:]:
    data = []; iter_cc = 0
    for cc in range(0,len(col_header)):
        data.append(row[0][iter_cc:iter_cc+col_sizes[cc]].replace('  ',''))
        iter_cc+=col_sizes[cc]+1
    call.append(data[3])
    names.append(data[4])
    lats.append(float(data[9]))
    lons.append(float(data[10]))
    elevs.append(float(data[11]))

#######################################
# Formatting the Cartopy plot
#######################################
#
cimgt.Stamen.get_image = image_spoof # reformat web request for street map spoofing
osm_img = cimgt.Stamen('terrain-background') # spoofed, downloaded street map

fig = plt.figure(figsize=(12,9)) # open matplotlib figure
ax1 = plt.axes(projection=osm_img.crs) # project using coordinate reference system (CRS) of street map
ax1.set_title('ASOS Station Map',fontsize=16)
extent = [-124.7844079,-66.9513812, 24.7433195, 49.3457868] # Contiguous US bounds
# extent = [-74.257159,-73.699215,40.495992,40.915568] # NYC bounds
ax1.set_extent(extent) # set extents
ax1.set_xticks(np.linspace(extent[0],extent[1],7),crs=ccrs.PlateCarree()) # set longitude indicators
ax1.set_yticks(np.linspace(extent[2],extent[3],7)[1:],crs=ccrs.PlateCarree()) # set latitude indicators
lon_formatter = LongitudeFormatter(number_format='0.1f',degree_symbol='',dateline_direction_label=True) # format lons
lat_formatter = LatitudeFormatter(number_format='0.1f',degree_symbol='') # format lats
ax1.xaxis.set_major_formatter(lon_formatter) # set lons
ax1.yaxis.set_major_formatter(lat_formatter) # set lats
ax1.xaxis.set_tick_params(labelsize=14)
ax1.yaxis.set_tick_params(labelsize=14)

scale = np.ceil(-np.sqrt(2)*np.log(np.divide((extent[1]-extent[0])/2.0,350.0))) # empirical solve for scale based on zoom
scale = (scale<20) and scale or 19 # scale cannot be larger than 19
ax1.add_image(osm_img, int(scale)) # add OSM with zoom specification

#######################################
# Plot the ASOS stations as points
#######################################
#
ax1.plot(lons, lats, markersize=5,marker='o',linestyle='',color='#3b3b3b',transform=ccrs.PlateCarree())

plt.show()

根据查看地图的屏幕大小，比例可能会显得模糊。如果需要更高分辨率的地图，可以使用更大的比例因子。此处使用比例因子 4 来保持此网页的加载时间最短。对于这里的应用程序，上图就足够了。如果散点图用于发布目的，则应将其放大到 6 甚至 7（这会导致更长的加载时间和更大的图像尺寸）。

要再次查看特定区域（例如纽约市）中的站点，可以更改范围（取消注释上面代码中的 NYC 范围）并且可以观察给定边界的 ASOS 站点。而不是在“Stamen()”函数中使用“地形背景”映射，而是将其简单地更改为：“地形”。对于 NYC 案例，结果如下：

由于大多数 ASOS 站点位于机场，因此很容易验证站点在纽约市的位置。三个大型机场：拉瓜迪亚机场、肯尼迪国际机场、纽瓦克机场 - 正如预期的那样，它们的上方都有红点。第四个点位于另一个成熟的位置：中央公园。第五个也是最后一个 ASOS 站位于泰特波罗机场。

如果位置特别重要，可以将文本注释添加到地图中以识别每个站点。这是在下图中完成的：

# Mapping ASOS Weather Stations Across New York City and Annotating Them
# This code uses a spoofing algorithm to avoid bounceback from OSM servers
# -- This code also parses lat/lon information from the ASOS station .txt
# -- file located here: ftp://ftp.ncdc.noaa.gov/pub/data/ASOS_Station_Photos/asos-stations.txt
#
import csv
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import io
from urllib.request import urlopen, Request
from PIL import Image

def image_spoof(self, tile): # this function pretends not to be a Python script
    url = self._image_url(tile) # get the url of the street map API
    req = Request(url) # start request
    req.add_header('User-agent','Anaconda 3') # add user agent to request
    fh = urlopen(req) 
    im_data = io.BytesIO(fh.read()) # get image
    fh.close() # close url
    img = Image.open(im_data) # open image with PIL
    img = img.convert(self.desired_tile_form) # set image format
    return img, self.tileextent(tile), 'lower' # reformat for cartopy

################################
# parsing the ASOS coordinates
################################
#
asos_data = []
with open('asos-stations.txt','r') as dat_file:
    reader = csv.reader(dat_file)
    for row in reader:
        asos_data.append(row)

row_delin = asos_data[3][0].split(' ')[:-1]
col_sizes = [len(ii) for ii in row_delin]

col_header = []; iter_ii = 0
for ii,jj in enumerate(col_sizes):
    col_header.append(asos_data[2][0][iter_ii:iter_ii+col_sizes[ii]].replace(' ',''))
    iter_ii+=col_sizes[ii]+1
    
call,names,lats,lons,elevs = [],[],[],[],[]
for row in asos_data[4:]:
    data = []; iter_cc = 0
    for cc in range(0,len(col_header)):
        data.append(row[0][iter_cc:iter_cc+col_sizes[cc]].replace('  ',''))
        iter_cc+=col_sizes[cc]+1
    call.append(data[3])
    names.append(data[4])
    lats.append(float(data[9]))
    lons.append(float(data[10]))
    elevs.append(float(data[11]))

#######################################
# Formatting the Cartopy plot
#######################################
#
cimgt.Stamen.get_image = image_spoof # reformat web request for street map spoofing
osm_img = cimgt.Stamen('terrain') # spoofed, downloaded street map

fig = plt.figure(figsize=(12,9)) # open matplotlib figure
ax1 = plt.axes(projection=osm_img.crs) # project using coordinate reference system (CRS) of street map
ax1.set_title('ASOS Station Map',fontsize=16)
# extent = [-124.7844079,-66.9513812, 24.7433195, 49.3457868] # Contiguous US bounds
extent = [-74.257159,-73.699215,40.495992,40.915568] # NYC bounds
ax1.set_extent(extent) # set extents
ax1.set_xticks(np.linspace(extent[0],extent[1],7),crs=ccrs.PlateCarree()) # set longitude indicators
ax1.set_yticks(np.linspace(extent[2],extent[3],7)[1:],crs=ccrs.PlateCarree()) # set latitude indicators
lon_formatter = LongitudeFormatter(number_format='0.1f',degree_symbol='',dateline_direction_label=True) # format lons
lat_formatter = LatitudeFormatter(number_format='0.1f',degree_symbol='') # format lats
ax1.xaxis.set_major_formatter(lon_formatter) # set lons
ax1.yaxis.set_major_formatter(lat_formatter) # set lats
ax1.xaxis.set_tick_params(labelsize=14)
ax1.yaxis.set_tick_params(labelsize=14)

scale = np.ceil(-np.sqrt(2)*np.log(np.divide((extent[1]-extent[0])/2.0,350.0))) # empirical solve for scale based on zoom
scale = (scale<20) and scale or 19 # scale cannot be larger than 19
ax1.add_image(osm_img, int(scale)) # add OSM with zoom specification

#######################################
# Plot the ASOS stations as points
#######################################
#
ax1.plot(lons, lats, markersize=10,marker='o',linestyle='',color='#b30909',transform=ccrs.PlateCarree(),label='ASOS Station')
transform = ccrs.PlateCarree()._as_mpl_transform(ax1) # set transform for annotations
# text annotations
for lat_s,lon_s,name_s in zip(lats,lons,names): # loop through lats/lons/names
    if lon_s>extent[0] and lon_s<extent[1] and lat_s>extent[2] and lat_s<extent[3]:
        print(name_s)
        # annotations, with some random placement to avoid overlap
        ax1.text(lon_s, lat_s+(np.random.randn(1)*0.01),name_s, {'color': 'k', 'fontsize': 8},
                 horizontalalignment='center', verticalalignment='bottom',
                 clip_on=False,transform=transform,bbox=dict(boxstyle="round",
                   ec='#121212', fc='#fadede'))

ax1.legend() # add a legend 
plt.show() # show plot

Cartopy 是一个强大的制图库，作为 Python 的许多可视化工具箱的一部分而存在。在本教程中，给出了 Anaconda 3 和 Cartopy 的安装过程，针对基于 Linux 的系统上地理映射的具体应用。在上面的第一个应用程序部分中，街道地图是针对给定的地理坐标绘制的。世贸中心一号地图绘制在两张不同的街道地图上：一张是绘制的，一张是基于图像的。在第二个应用程序中，自动地面观测系统 (ASOS) 坐标点被映射到整个美国，以将真实数据应用于 Cartopy 功能。最后，通过大纽约地区的文本注释进一步探索了 ASOS 坐标。通过查看与每个坐标相关的名称并将每个坐标与相应的附近机场或气象站相关联，在纽约市验证了五个 ASOS 坐标。在接下来的教程中，将继续探索 Cartopy，特别是对于涉及城市 shapefile 和卫星数据的实际应用。