Python-地理空间分析学习手册-四-

布客飞龙
70 阅读51分钟

Python 地理空间分析学习手册（四）

原文：Learning Geospatial Analysis with Python

协议：CC BY-NC-SA 4.0

九、实时数据

地理空间分析师的一句俗语是:地图一创建就过时了。这句话反映了地球和地球上的一切都在不断变化的事实。对于地理空间分析的大部分历史和本书的大部分内容，地理空间产品是相对静态的。原始数据集通常在几个月到几年的时间内更新。地图中地理空间数据的时代被称为数据货币。

由于收集数据所需的时间和费用，数据货币传统上不是主要关注点。网络地图、无线蜂窝调制解调器和低成本的全球定位系统天线改变了这一焦点。现在，监控一个快速变化的对象或系统，并将这些变化在线传播给数百万人，在后勤上是可行的，甚至是非常经济的。这一变化正在彻底改变地理空间技术，并将其带入新的方向。这场革命最直接的证据是使用谷歌地图或 OpenLayers 等系统和网络可访问数据格式的网络地图混搭。每天，越来越多的电子设备上线，广播它们的位置和数据，以实现自动化或远程控制。例子包括恒温器、照相机、汽车等等。你也可以使用廉价的嵌入式计算机，比如流行的树莓皮，将几乎任何东西都变成一个互联的智能设备。这种将设备连接成数据和信息网络的概念被称为物联网 ( 物联网)。

在本章中，我们将查看以下主题:

  • 实时数据的局限性
  • 使用实时数据
  • 跟踪车辆
  • 风暴追逐
  • 来自现场的报告

到最后，您将学会使用实时地理空间数据，并且能够构建一个字段报告工具，该工具可以作为任何类型数据的数据传输源。

技术要求

本章要求以下内容:

实时数据的局限性

术语实时数据通常意味着接近实时。一些跟踪设备捕捉实时数据，并且可能每秒更新几次。但是，广播这些数据的基础设施的局限性可能会将输出限制在每 10 秒或更长时间一次。天气雷达就是一个完美的例子。一台多普勒天气雷达 ( DWR )持续扫描，但数据通常每五分钟在线一次。但是考虑到与传统地理空间数据更新的对比，几分钟的刷新就足够实时了。局限性可以总结如下:

  • 限制数据大小的网络带宽限制
  • 限制数据更新频率的网络延迟
  • 由于电池寿命等限制，数据源的可用性
  • 由于消费者可以立即获得数据，因此缺乏质量控制
  • 由于快速接收未经验证的数据而导致的安全漏洞

实时数据为地理空间应用开辟了更多的机会，因此我们接下来将研究如何使用它。

使用实时数据

Web 混搭通常使用实时数据。网络混搭令人惊叹，改变了许多不同行业的运作方式。但是它们通常是有限的，因为它们通常只在地图上显示预处理过的数据，并允许开发人员访问 JavaScript API。但是如果你想以某种方式处理数据呢？如果您想要过滤、更改，然后将其发送到另一个系统，该怎么办？要将实时数据用于地理空间分析，您需要能够将其作为点数据或地理参考栅格进行访问。

You can find out more about web map mashups here: www.esri.com/arcgis-blog….

与前面几章中的例子一样，脚本尽可能简单，设计为从头到尾阅读，没有太多的思维循环。当使用函数时，首先列出它们，然后是脚本变量声明，最后是主程序执行。

现在让我们看看如何使用 NextBus API 中的车辆访问实时和点位置数据源。

跟踪车辆

对于我们的第一个实时数据源，我们将使用优秀的 NextBus API 。next bus(www.nextbus.com/)是一项商业服务，跟踪城市公共交通，包括公共汽车、无轨电车和火车。乘坐这些公交线路的人们可以追踪下一班公交车的到达时间。

更好的是，在客户的允许下，NextBus 通过表征状态转移 ( REST ) **API 发布跟踪数据。**使用 URL API 调用，开发者可以请求关于车辆的信息，并接收关于其位置的 XML 文档。这个应用编程接口是开始使用实时数据的简单方法。

如果你去 NextBus，你会看到如下截图所示的网络界面，显示了加州洛杉矶市地铁系统的数据:

系统允许您选择几个参数来学习下一站的当前位置和时间预测。在屏幕的右侧，有一个指向谷歌地图混搭的链接，显示了特定路线的交通跟踪数据，如下图所示:

这是一个非常有用的网站，但它不能让我们控制数据如何显示和使用。让我们直接使用 Python 和 NextBus REST API 访问原始数据，开始处理实时数据。

对于本章中的示例，我们将使用在此找到的有文档记录的 next bus API:www.nextbus.com/xmlFeedDocs…。

从这个例子开始，我们需要一个所需总线的列表。

下一个巴士公司名单

NextBus 的客户被称为机构。在我们的例子中，我们将在前往加州洛杉矶的路线上跟踪公共汽车。首先，我们需要获得一些关于该机构的信息。下一个总线应用编程接口由一个名为publicXMLFeed的网络服务组成，在这个服务中你可以设置一个名为command的参数。我们将在浏览器中调用agencyList命令，使用以下 REST URL 获取包含机构信息的 XML 文档:webservices.nextbus.com/service/pub…。

当我们在浏览器中访问该链接时，它会返回一个包含<agency/>标签的 XML 文档。洛杉矶的标签如下所示:

<agency tag="lametro" title="Los Angeles Metro" regionTitle="California-Southern"/>

既然我们有了公共汽车的清单，我们需要知道它们可以行驶的路线。

下一班车路线列表

tag属性是雷霆湾的 ID，我们需要它来执行其他的 NextBus API 命令。其他属性是人类可读的元数据。我们需要的下一条信息是关于路线 2 巴士路线的细节。为了获取这些信息，我们将使用机构标识和routeList REST 命令，通过将网址粘贴到我们的网络浏览器中来获取另一个 XML 文档。

Note that the agency ID is set to the parameter in the REST URL: webservices.nextbus.com/service/pub….

当我们在浏览器中调用这个网址时，我们会得到以下 XML 文档:

<?xml version="1.0" encoding="utf-8" ?>
<body copyright="All data copyright Los Angeles Metro 2015."><route tag="2" title="2 Downtown LA - Pacific Palisades Via"/><route tag="4" title="4 Downtown LA - Santa Monica Via Santa"/>
<route tag="10" title="10 W Hollywood-Dtwn LA -Avalon Sta Via"/>
...
<route tag="901" title="901 Metro Orange Line"/>
<route tag="910" title="910 Metro Silver Line"/>
</body>

我们有公共汽车和路线。我们准备开始追踪他们的位置！

NextBus 车辆位置

所以根据这些结果，存储在tag属性中的主线路线 ID 简单来说就是1。因此，现在，我们已经掌握了沿着洛杉矶地铁 2 号线追踪公交车所需的所有信息。

只有一个必需的参数(称为t)表示自 1970 年纪元日期(1970 年 1 月 1 日，世界协调时午夜以来的毫秒数。纪元日期只是机器用来跟踪时间的计算机标准。在 NextBus API 内最容易做的事情就是为这个值指定0，它返回最后 15 分钟的数据。

有一个可选的direction标签，允许您在一条路线上有多辆相反方向行驶的公共汽车的情况下指定一个终止公共汽车站。但是，如果我们不指定，API 将返回第一个，这符合我们的需要。获取洛杉矶地铁主线路线的 REST 网址如下:webservices.nextbus.com/service/pub… 。

在浏览器中调用此 REST URL 会返回以下 XML 文档:

<?xml version="1.0" encoding="utf-8" ?>
<body copyright="All data copyright Los Angeles Metro 2015.">
<vehicle id="7582" routeTag="2" dirTag="2_758_0" lat="34.097992" lon="-118.350365" secsSinceReport="44" predictable="true" heading="90" speedKmHr="0"/>
<vehicle id="7583" routeTag="2" dirTag="2_779_0" lat="34.098076" lon="-118.301399" secsSinceReport="104" predictable="true" heading="90" speedKmHr="37"/>
. . .
</body >

每个vehicle标记代表最近 15 分钟内的一个位置。last标签是最近的位置(尽管 XML 在技术上是无序的)。

这些公共交通系统不会一直运行。许多商店在当地时间晚上 10 点(22 点)关门。如果您在脚本中遇到错误，请使用 NextBus 网站找到正在运行的系统，并更改代理和路由变量到该系统。

我们现在可以编写一个 Python 脚本，返回给定路线上的公共汽车位置。如果我们不指定direction标签，NextBus 将返回第一个标签。在本例中，我们将通过使用前面章节中演示的内置 Python urllib库调用 REST URL 来轮询 NextBus 跟踪 API。

我们将使用简单的内置minidom模块解析返回的 XML 文档，也显示在迷你文档模块部分、第 4 章T5、地理空间 Python 工具箱中。这个脚本只是输出路线 2 公共汽车的最新纬度和经度。您将在顶部附近看到代理和路由变量。为此，我们需要遵循以下步骤:

  1. 首先，我们导入我们需要的库:
import urllib.request
import urllib.parse
import urllib.error
from xml.dom import minidom
  1. 现在，我们为应用编程接口模式以及要查询的客户和路线设置变量:
# Nextbus API command mode
command = "vehicleLocations"

# Nextbus customer to query
agency = "lametro"

# Bus we want to query
route = "2"
  1. 我们将时间值设置为0，这将抓取最后15分钟的数据:
# Time in milliseconds since the
# 1970 epoch time. All tracks
# after this time will be returned.
# 0 only returns data for the last
# 15 minutes
epoch = "0"
  1. 现在，我们需要构建用于访问应用编程接口的查询网址:
# Build our query url
# webservices base url
url = "http://webservices.nextbus.com"

# web service path
url += "/service/publicXMLFeed?"

# service command/mode
url += "command=" + command

# agency
url += "&a=" + agency
url += "&r=" + route
url += "&t=" + epoch
  1. 接下来，我们可以使用urllib调用 API:
# Access the REST URL
feed = urllib.request.urlopen(url)
if feed:
 # Parse the xml feed
 xml = minidom.parse(feed)
 # Get the vehicle tags
 vehicles = xml.getElementsByTagName("vehicle")
 # Get the most recent one. Normally there will
 # be only one.
  1. 最后，我们可以访问结果并打印出每辆公共汽车的位置:
 if vehicles:
   bus = vehicles.pop()
   # Print the bus latitude and longitude
   att = bus.attributes
   print(att["lon"].value, ",", att["lat"].value)
 else:
   print("No vehicles found.")

这个脚本的输出只是一个纬度和经度值，这意味着我们现在可以控制并理解 API。输出应该是纬度和经度的坐标值。

现在我们准备使用这些位置值来创建我们自己的地图。

映射下一个总线位置

免费获取街道地图数据的最佳来源是开放街道地图 ( OSM )项目:www.openstreetmap.org。OSM 还有一个可公开获取的 REST API，用于创建静态地图图像，名为 StaticMapLite : 。

OSM 静态地图应用编程接口提供了一个基于谷歌静态地图应用编程接口的GET应用编程接口，以创建具有有限数量的点标记和线的简单地图图像。一个GET API，与 REST 相反，API 允许你在网址上的问号后面附加名称/值参数对。REST 应用编程接口将参数作为网址路径的一部分。我们将使用该应用编程接口按需创建我们自己的下一个总线应用编程接口地图，带有一个红色的图钉图标用于总线位置。

在下一个例子中，我们已经将之前的脚本压缩为一个名为nextbus()的紧凑函数。nextbus()函数接受代理、路线、命令和纪元作为参数。该命令默认为vehicleLocations，纪元默认为0，以获取最后 15 分钟的数据。在这个脚本中，我们将传入 LA route-2 路由信息，并使用默认命令返回最近的总线纬度/经度。

我们有第二个名为nextmap()的函数，它会创建一个地图，在每次调用公共汽车时的当前位置都有一个紫色标记。地图是通过为 OSM StaticMapLite应用编程接口建立一个GET网址创建的，该网址以公共汽车的位置为中心，并使用 1-18 和地图大小之间的缩放级别来确定地图范围。

You can access the API directly in a browser to see an example of what the nextmap() function does. You will need a free MapQuest Developer API key available by registering here: developer.mapquest.com/plan_purcha…. Once you have the key, insert it in the key parameter where it says YOUR_API_KEY_HERE. Then, you can test the following example URL: https://www.mapquestapi.com/staticmap/v4/getmap?size=865,512&type=map&pois=mcenter,40.702147,-74.015794|&zoom=14&center=40.714728,-73.998672&imagetype=JPEG&key=YOUR_API_KEY_HERE.

静态地图如下所示:

nextmap()函数接受地图基础图像名称的下一个公交机构标识、路线标识和字符串。该函数调用nextbus()函数获取纬度/经度对。该程序的执行以定时间隔循环进行，在第一次通过时创建一个地图，然后在后续通过时覆盖该地图。该程序还会在每次保存地图时输出时间戳。requests变量指定通过次数，freq变量表示每次循环之间的时间(秒)。让我们检查下面的代码，看看这个例子是如何工作的:

  1. 首先，我们导入我们需要的库:
import urllib.request
import urllib.parse
import urllib.error
from xml.dom import minidom
import time
  1. 接下来，我们创建一个函数，它可以获取给定路线上公共汽车的最新位置:
def nextbus(a, r, c="vehicleLocations", e=0):
 """Returns the most recent latitude and
 longitude of the selected bus line using
 the NextBus API (nbapi)
 Arguments: a=agency, r=route, c=command,
 e=epoch timestamp for start date of track,
 0 = the last 15 minutes"""
 nbapi = "http://webservices.nextbus.com"
 nbapi += "/service/publicXMLFeed?"
 nbapi += "command={}&a={}&r={}&t={}".format(c, a, r, e)
 xml = minidom.parse(urllib.request.urlopen(nbapi))
 # If more than one vehicle, just get the first 
 bus = xml.getElementsByTagName("vehicle")[0]
 if bus: 
 at = bus.attributes
 return(at["lat"].value, at["lon"].value)
 else:
 return (False, False)
  1. 现在，我们有了一个在地图图像上绘制公交位置的功能:
def nextmap(a, r, mapimg):
 """Plots a nextbus location on a map image
 and saves it to disk using the MapQuest OpenStreetMap Static Map 
 API (osmapi)"""
 # Fetch the latest bus location
 lat, lon = nextbus(a, r)
 if not lat:
   return False
 # Base url + service path
  1. 在该函数中，我们在 URL 中设置了 API 参数:
 osmapi = "https://www.mapquestapi.com/staticmap/v4/getmap?
 type=map&"
# Use a red, pushpin marker to pin point the bus
 osmapi += "mcenter={},{}|&".format(lat, lon)
# Set the zoom level (between 1-18, higher=lower scale)
 osmapi += "zoom=18&"
# Center the map around the bus location
 osmapi += "center={},{}&".format(lat, lon)
# Set the map image size
 osmapi += "&size=1500,1000"
# Add our API Key
 osmapi += "&key=YOUR_API_KEY_HERE"
  1. 现在，我们可以通过调用网址来创建图像并保存它:
 # Create a PNG image
 osmapi += "imagetype=png&"
 img = urllib.request.urlopen(osmapi)

# Save the map image
 with open("{}.png".format(mapimg), "wb") as f:
   f.write(img.read())
return True
  1. 现在，在我们的主程序中，我们可以设置要跟踪的公共汽车的变量:
# Nextbus API agency and bus line variables
agency = "lametro"
route = "2"
# Name of map image to save as PNG
nextimg = "nextmap"
  1. 然后，我们可以指定我们想要的跟踪点的数量和频率:
# Number of updates we want to make
requests = 1
# How often we want to update (seconds)
freq = 5
  1. 最后，我们可以开始跟踪和更新我们的地图图像:
# Map the bus location every few seconds
for i in range(requests):
 success = nextmap(agency, route, nextimg)
 if not success:
   print("No data available.")
   continue
 print("Saved map {} at {}".format(i, time.asctime()))
 time.sleep(freq)
  1. 当脚本运行时，您将看到类似以下内容的输出，显示脚本保存每个地图的时间:
Saved map 0 at Sun Nov 1 22:35:17 2015
Saved map 1 at Sun Nov 1 22:35:24 2015
Saved map 2 at Sun Nov 1 22:35:32 2015

该脚本会保存一个类似于以下内容的地图图像，具体取决于运行该脚本时公共汽车的位置:

这张地图是使用应用编程接口创建自定义地图产品的一个很好的例子。但这是一个非常基本的跟踪应用。为了开始将它开发成一个更有趣的地理空间产品，我们需要将它与其他一些实时数据源结合起来，让我们能够更好地了解情况。

现在我们可以跟踪公共汽车了，让我们在地图上添加一些对乘坐公共汽车的乘客有用的附加信息。让我们添加一些天气数据。

风暴追逐

到目前为止，我们已经创建了一个更简单版本的 NextBus 网站。但我们这样做的方式最终让我们完全控制了产量。现在，我们希望使用这个控件来超越 NextBus 谷歌地图混搭的功能。我们将添加另一个对旅行者和公交运营商都非常重要的实时数据源:天气。

爱荷华州立大学的 Mesonet 项目为应用程序提供免费且经过打磨的天气数据。我们使用这些数据为我们的公交位置图创建实时天气图。我们可以使用开放地理空间联盟 ( OGC ) 网络地图服务 ( WMS )标准来请求我们感兴趣区域的单个图像。WMS 是通过网络提供地理参考地图图像的 OGC 标准；它们是由地图服务器通过 HTTP 请求生成的。

Mesonet 系统提供了一个出色的网络地图服务，它根据一个格式正确的 WMS 请求，从全球降水镶嵌图中返回一个子集图像。这种请求的一个例子是以下查询: mesonet.agron .传道论. edu/cgi-bin/wms…

因为本章中的示例依赖于实时数据，所以如果感兴趣的区域没有活动，列出的特定请求可能会生成空白天气图像。您可以访问此链接(radar.weather.gov/ridge/Conus…)来查找正在发生风暴的区域。该页面包含谷歌地球或 QGIS 的 KML 链接。这些 WMS 图像是透明的 PNG 图像，类似于以下示例:

另一方面，OSM 网站不再通过 WMS 提供街道地图——只提供瓷砖。但是，它们允许其他组织下载切片或原始数据来扩展免费服务。美国国家海洋和大气管理局 ( 美国国家海洋和大气管理局)已经做到了这一点，并为他们的 OSM 数据提供了一个 WMS 接口，允许请求检索我们的公交路线所需的单个底图图像:

我们现在有了获取底图和天气数据的数据源。我们想把这些图像结合起来，画出公共汽车的当前位置。这次，我们将不再使用简单的点，而是变得更加复杂，并添加以下总线图标:

您需要从这里将此图标busicon.png下载到您的工作目录中:https://github . com/GeospatialPython/Learn/blob/master/busicon . png？原始=真实。

现在，我们将结合以前的脚本和新的数据源来创建实时天气巴士地图。因为我们要混合街道图和天气图，我们需要前几章使用的 Python 图像库 ( PIL )。我们将用一个简单的wms()函数替换前面例子中的nextmap()函数，该函数可以通过任何 WMS 服务的边界框抓取地图图像。我们还将添加一个将十进制度数转换为米的函数，名为ll2m()。

该脚本获取公交车位置，将该位置转换为米，在该位置周围创建一个 2 英里(3.2 公里)的矩形，然后下载街道和天气地图。然后使用 PIL 将地图图像混合在一起。然后，PIL 将公交图标图像缩小到 30 x 30 像素，并将其粘贴到地图的中心，也就是公交位置。让我们看看下面的代码是如何工作的:

  1. 首先，我们将导入我们需要的库:
import sys
import urllib.request
import urllib.parse
import urllib.error
from xml.dom import minidom
import math
try:
 import Image
except:
 from PIL import Image
  1. 现在我们将重用前面例子中的nextbus函数来获取总线跟踪数据:
def nextbus(a, r, c="vehicleLocations", e=0):
 """Returns the most recent latitude and
 longitude of the selected bus line using
 the NextBus API (nbapi)"""
 nbapi = "http://webservices.nextbus.com"
 nbapi += "/service/publicXMLFeed?"
 nbapi += "command=%s&a=%s&r=%s&t=%s" % (c, a, r, e)
 xml = minidom.parse(urllib.request.urlopen(nbapi))
 # If more than one vehicle, just get the first
 bus = xml.getElementsByTagName("vehicle")[0]
 if bus:
 at = bus.attributes
 return(at["lat"].value, at["lon"].value)
 else:
 return (False, False)
  1. 我们还需要一个将纬度和经度转换为米的函数:
def ll2m(lon, lat):
 """Lat/lon to meters"""
 x = lon * 20037508.34 / 180.0
 y = math.log(math.tan((90.0 + lat) *
 math.pi / 360.0)) / (math.pi / 180.0)
 y = y * 20037508.34 / 180
 return (x, y)
  1. 现在，我们需要一个函数来检索 WMS 地图图像，这将用于我们的天气图像:
def wms(minx, miny, maxx, maxy, service, lyr, epsg, style, img, w, 
        h):
    """Retrieve a wms map image from
    the specified service and saves it as a JPEG."""
    wms = service
    wms += "?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&"
    wms += "LAYERS={}".format(lyr)
    wms += "&STYLES={}&".format(style)
    wms += "SRS=EPSG:{}&".format(epsg)
    wms += "BBOX={},{},{},{}&".format(minx, miny, maxx, maxy)
    wms += "WIDTH={}&".format(w)
    wms += "HEIGHT={}&".format(h)
    wms += "FORMAT=image/jpeg"
    wmsmap = urllib.request.urlopen(wms)
    with open(img + ".jpg", "wb") as f:
        f.write(wmsmap.read())
  1. 现在我们可以在主程序中设置所有变量来使用我们的函数:
# Nextbus agency and route ids
agency = "roosevelt"
route = "shuttle"
# OpenStreetMap WMS service
basemap = "http://ows.mundialis.de/services/service"
# Name of the WMS street layer
streets = "TOPO-OSM-WMS"
# Name of the basemap image to save
mapimg = "basemap"
# OpenWeatherMap.org WMS Service
weather = "https://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0q.cgi?"
# If the sky is clear over New York,
# use the following url which contains
# a notional precipitation sample:
# weather = "http://git.io/vl4r1"
# WMS weather layer
weather_layer = "nexrad-n0q-900913"
# Name of the weather image to save
skyimg = "weather"
# Name of the finished map to save
final = "next-weather"
# Transparency level for weather layer
# when we blend it with the basemap.
# 0 = invisible, 1 = no transparency
opacity = .5
# Pixel width and height of the
# output map images
w = 600
h = 600
# Pixel width/height of the the
# bus marker icon
icon = 30
  1. 现在我们准备好了我们的巴士位置:
# Get the bus location
lat, lon = nextbus(agency, route)
if not lat:
 print("No bus data available.")
 print("Please try again later")
 sys.exit()
# Convert strings to floats
lat = float(lat)
lon = float(lon)
# Convert the degrees to Web Mercator
# to match the NOAA OSM WMS map
x, y = ll2m(lon, lat)
# Create a bounding box 1600 meters
# in each direction around the bus
minx = x - 1600
maxx = x + 1600
miny = y - 1600
maxy = y + 1600
  1. 然后，我们可以下载我们的街道地图:
# Download the street map
wms(minx, miny, maxx, maxy, basemap, streets, mapimg, w, h)
  1. 然后，我们可以下载天气图:
# Download the weather map
wms(minx, miny, maxx, maxy, weather, weather_layer, skyimg, w, h)
  1. 现在我们可以将天气数据叠加在公交地图上:
# Open the basemap image in PIL
im1 = Image.open("basemap.png").convert('RGBA')
# Open the weather image in PIL
im2 = Image.open("weather.png").convert('RGBA')
# Convert the weather image mode
# to "RGB" from an indexed PNG
# so it matches the basemap image
im2 = im2.convert(im1.mode)
# Create a blended image combining
# the basemap with the weather map
im3 = Image.blend(im1, im2, opacity)
  1. 接下来，我们需要将巴士图标添加到我们的组合地图中，以显示巴士的位置:
# Open up the bus icon image to
# use as a location marker.
# http://git.io/vlgHl 
im4 = Image.open("busicon.png")
# Shrink the icon to the desired
# size
im4.thumbnail((icon, icon))
# Use the blended map image
# and icon sizes to place
# the icon in the center of
# the image since the map
# is centered on the bus
# location.
w, h = im3.size
w2, h2 = im4.size
# Paste the icon in the center of the image
center_width = int((w/2)-(w2/2))
center_height = int((h/2)-(h2/2))
im3.paste(im4, (center_width, center_height), im4)
  1. 最后，我们可以保存完成的地图:
# Save the finished map
im3.save(final + ".png")

该脚本将生成类似于以下内容的地图:

地图向我们显示，公共汽车在当前位置正经历中度降水。如之前的 Mesonet 网站截图所示，颜色渐变的范围从淡蓝色代表轻度降水，然后是绿色、黄色、橙色，随着雨越下越大，颜色渐变的范围会变成红色(或者从浅灰色变成黑色和白色的深灰色)。因此，在创建这张地图时，公交线路运营者可以使用这张图片告诉他们的司机开得慢一点，乘客会知道他们可能想在去公交车站之前拿把伞。

因为我们想在低层次上学习 NextBus API，所以我们直接使用内置的 Python 模块来使用该 API。但是该应用编程接口有几个第三方 Python 模块，包括 PyPI 上的一个，简称为nextbus，它允许您为所有 NextBus 命令处理更高级别的对象，并提供本章简单示例中没有的更健壮的错误处理。

现在我们已经学会了如何检查天气，让我们使用 Python、HTML 和 JavaScript 将离散的实时数据源组合成更有意义的产品。

来自现场的报告

在这一章的最后一个例子中，我们将下车到野外。现代智能手机、平板电脑和笔记本电脑允许我们更新地理信息系统，并从任何地方查看这些更新。我们将使用 HTML、GeoJSON、小叶 JavaScript 库和名为 leaf 的纯 Python 库来创建一个客户端-服务器应用程序，该应用程序允许我们将地理空间信息发布到服务器，然后创建一个交互式网络地图来查看这些数据更新。

首先，我们需要一个 web 表单，显示您的当前位置，并在您提交表单时更新服务器，其中包含关于您的位置的注释。你可以在这里找到表格:geospatialpython.github.io/Learn/field…。

下面的截图显示了表单:

您可以查看该表单的源代码，了解其工作原理。映射是使用传单库完成的，并将地理信息发布到 myjson.com 的一个唯一的网址上。您可以在移动设备上使用此页面，将其移动到任何 web 服务器，甚至可以在本地硬盘上使用。

该表单在myjson.com上公开发布到以下网址:api.myjson.com/bins/467pm。您可以在浏览器中访问该网址来查看原始地理信息。

接下来，您需要从 PyPI 安装叶库。leaf 为创建 leaf 网络地图提供了一个简单的 Python 应用编程接口。你可以在这里找到更多关于叶的信息:github.com/python-visu…。

leaf 使得制作一张 leaf 地图变得非常简单。这个脚本只有几行，将输出一个名为map.html的网页。我们将 GeoJSON URL 传递给map对象，该对象将在地图上绘制位置:

import folium
m = folium.Map()
m.geo_json(geo_path="https://api.myjson.com/bins/467pm")
m.create_map(path="map.html")

生成的交互式地图会将点显示为标记。单击标记时，将显示表单中的信息。你可以在任何浏览器中打开 HTML 文件。

摘要

实时数据是进行新类型地理空间分析的一种令人兴奋的方式，只是最近几项不同技术的进步才使其成为可能，包括网络地图、全球定位系统和无线通信。在本章中，您学习了如何访问实时位置数据的原始源，如何获取实时栅格数据的子集，如何仅使用 Python 将不同类型的实时数据组合到自定义地图分析产品中，以及如何构建客户端-服务器地理空间应用程序来实时更新地理信息系统。

与前几章一样，这些示例包含构建块，允许您使用 Python 构建新类型的应用程序，这些应用程序远远超出了典型的流行且无处不在的基于 JavaScript 的混搭。

在下一章中，我们将把迄今为止所学的一切结合成一个完整的地理空间应用程序，在现实场景中应用算法和概念。

十、把这一切放在一起

在整本书中，我们已经触及了地理空间分析的所有重要方面，并在 Python 中使用了各种不同的技术来分析不同类型的地理空间数据。在这最后一章中，我们将利用几乎所有我们已经讨论过的主题来制作一个非常受欢迎的现实世界产品:全球定位系统路线分析报告。

这些报告常见于数十种移动应用服务、全球定位系统手表、车载导航系统和其他基于全球定位系统的工具。全球定位系统通常记录位置、时间和高度。从这些值中，我们可以得出大量的辅助信息，这些信息与记录数据的路线沿线发生的情况有关。包括 RunKeeper、MapMyRun、Strava 和 Nike Plus 在内的健身应用程序都使用类似的报告来呈现来自跑步、徒步旅行、骑自行车和步行的全球定位系统跟踪的锻炼数据。

我们将使用 Python 创建其中一个报告。这个程序有将近 500 行代码，是我们迄今为止最长的，所以我们将一点一点地来看。我们将结合以下技术:

  • 理解典型的全球定位系统报告
  • 构建全球定位系统报告工具

当我们逐步完成这个程序时，所有使用的技术都将是熟悉的，但是我们将以新的方式使用它们。

技术要求

本章我们需要以下内容:

  • Python 3.6 或更高版本
  • 内存:最小–6gb(Windows)，8gb(macOS)；推荐 8 GB
  • 存储:最低 7200 转/分的 SATA，20 GB 可用空间，推荐的固态硬盘，40 GB 可用空间
  • 处理器:最低英特尔酷睿 i3 2.5 GHz，推荐英特尔酷睿 i5
  • PIL:Python 图像库
  • 多维数组处理库
  • pygooglechart:优秀的谷歌图表 API 的 Python 包装器
  • FPDF:一个简单而纯粹的 Python PDF 作者

理解典型的全球定位系统报告

典型的全球定位系统报告具有共同的元素，包括路线图、高程剖面图和速度剖面图。以下截图是通过 RunKeeper(runkeeper.com/index)记录的典型路线的报告:

我们的报告将是相似的，但我们将增加一个转折。我们将像这项服务一样包括路线图和海拔剖面图，但我们还将添加记录路线时发生的天气情况以及在路线上拍摄的地理位置照片。

现在我们已经知道了什么是全球定位系统报告，让我们学习如何构建它。

构建全球定位系统报告工具

我们节目的名字是GPX-Reporter.py。如果你还记得第 2 章学习地理空间数据中的标记和基于标记的格式一节， GPX 格式是存储全球定位系统路线信息的最常见方式。几乎每一个依靠全球定位系统数据的程序和设备都可以在 GPX 之间进行转换。

对于本例，您可以从:git.io/vl7qi下载一个示例 GPX 文件。此外，您还需要安装一些来自 PyPI 的 Python 库。

你应该简单地使用easy_install或pip来安装这些工具。我们还将使用一个名为SRTM.py的模块。该模块用于处理奋进号航天飞机在 2000 年为期 11 天的航天飞机雷达地形任务 ( SRTM )中收集的近全球高程数据。使用pip安装 SRTM 模块:

pip install srtm.py

或者，您也可以下载压缩文件，将其解压缩，并将srtm文件夹复制到您的 Python site-packages目录或您的工作目录:git.io/vl5Ls。

您还需要注册一个免费的黑暗天空应用编程接口。这项免费服务提供独特的工具。这是唯一一项为几乎任何地点提供全球历史天气数据的服务，每天最多可免费提供 1000 个请求:darksky.net/dev。

黑暗天空将为您提供一个文本键，您可以在运行 GPX-记者程序之前将其插入到名为api_key的变量中。最后，根据黑暗天空的服务条款，您需要下载一个徽标图像以插入报告中:https://raw . githubusercontent . com/GeospatialPython/Learn/master/darksky . png。

You can review the Dark Sky Terms of Service here: darksky.net/dev/docs/te….

现在，我们准备通过 GPX-记者项目开展工作。像本书中的其他脚本一样，这个程序试图最小化函数，这样您就可以更好地跟踪程序，并更轻松地修改它。以下列表包含程序中的主要步骤:

  1. 设置 Python logging模块
  2. 建立我们的助手函数
  3. 解析 GPX 数据文件
  4. 计算路线边界框
  5. 缓冲边界框
  6. 将盒子转换为仪表
  7. 下载底图
  8. 下载高程数据
  9. 对高程数据进行山体阴影
  10. 增加山体阴影对比度
  11. 混合山体阴影和底图
  12. 在单独的图像上绘制 GPX 轨迹
  13. 混合轨迹图像和底图
  14. 绘制起点和终点
  15. 保存地图图像
  16. 计算路线英里标记
  17. 构建高程剖面图
  18. 获取路线时间段的天气数据
  19. 生成 PDF 报告

下一小节将带您完成第一步。

初始设置

程序的开头是import语句，后面是 Python logging模块。与简单的print语句相比，logging模块提供了一种更强大的方法来跟踪和记录程序状态。在程序的这一部分，我们按照以下步骤进行配置:

  1. 我们首先需要安装我们需要的所有库，如下面的代码所示:
from xml.dom import minidom
import json
import urllib.request
import urllib.parse
import urllib.error
import math
import time
import logging
import numpy as np
import srtm # Python 3 version: http://git.io/vl5Ls
import sys
from pygooglechart import SimpleLineChart
from pygooglechart import Axis
import fpdf
import glob
import os
try:
 import Image
 import ImageFilter
 import ImageEnhance
 import ImageDraw
except:
 from PIL import Image
 from PIL import ImageFilter
 from PIL import ImageEnhance
 from PIL import ImageDraw
 from PIL.ExifTags import TAGS
  1. 现在我们可以配置 Python logging模块来告诉我们整个过程中发生了什么，如下所示:
# Python logging module.
# Provides a more advanced way
# to track and log program progress.
# Logging level - everything at or below
# this level will output. INFO is below.
level = logging.DEBUG
# The formatter formats the log message.
# In this case we print the local time, logger name, and message
formatter = logging.Formatter("%(asctime)s - %(name)s - %(message)s")
# Establish a logging object and name it
log = logging.getLogger("GPX-Reporter")
# Configure our logger
log.setLevel(level)
# Print to the command line
console = logging.StreamHandler()
console.setLevel(level)
console.setFormatter(formatter)
log.addHandler(console)

这个记录器打印到控制台，但是只需做一些简单的修改，您就可以将它打印到一个文件，甚至一个数据库，只需更改这一部分中的配置。该模块内置于 Python 中，并记录在这里:docs.python.org/3/howto/log…。

接下来，我们有几个在整个程序中多次使用的实用函数。

使用实用功能

除了与时间相关的功能之外，以下所有功能都以某种形式在前面的章节中使用过。让我们看看如何在我们的例子中使用实用函数:

  1. 首先，ll2m()功能将经纬度转换为米:
def ll2m(lat, lon):
 """Lat/lon to meters"""
 x = lon * 20037508.34 / 180.0
 y = math.log(math.tan((90.0 + lat) *
 math.pi / 360.0)) / (math.pi / 180.0)
 y = y * 20037508.34 / 180
 return (x, y)
  1. world2pixel()功能将地理空间坐标转换为输出地图图像上的像素坐标:
def world2pixel(x, y, w, h, bbox):
 """Converts world coordinates
 to image pixel coordinates"""
 # Bounding box of the map
 minx, miny, maxx, maxy = bbox
 # world x distance
 xdist = maxx - minx
 # world y distance
 ydist = maxy - miny
 # scaling factors for x, y
 xratio = w/xdist
 yratio = h/ydist
 # Calculate x, y pixel coordinate
 px = w - ((maxx - x) * xratio)
 py = (maxy-y) * yratio
 return int(px), int(py)
  1. 然后，我们有get_utc_epoch()和get_local_time()将存储在 GPX 文件中的 UTC 时间转换为沿线的当地时间:
def get_utc_epoch(timestr):
 """Converts a GPX timestamp to Unix epoch seconds
 in Greenwich Mean Time to make time math easier"""
 # Get time object from ISO time string
 utctime = time.strptime(timestr, '%Y-%m-%dT%H:%M:%S.000Z')
 # Convert to seconds since epoch
 secs = int(time.mktime(utctime))
 return secs
  1. 现在我们有了哈弗斯距离函数和简单的wms函数来检索地图图像:
def haversine(x1, y1, x2, y2):
 """Haversine distance formula"""
 x_dist = math.radians(x1 - x2)
 y_dist = math.radians(y1 - y2)
 y1_rad = math.radians(y1)
 y2_rad = math.radians(y2)
 a = math.sin(y_dist/2)**2 + math.sin(x_dist/2)**2 \
 * math.cos(y1_rad) * math.cos(y2_rad)
 c = 2 * math.asin(math.sqrt(a))
 # Distance in miles. Just use c * 6371
 # for kilometers
 distance = c * (6371/1.609344) # Miles
 return distance
  1. wms()功能通过以下代码检索地图图像:
def wms(minx, miny, maxx, maxy, service, lyr, epsg, style, img, w, h):
 """Retrieve a wms map image from
 the specified service and saves it as a JPEG."""
 wms = service
 wms += "?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&"
 wms += "LAYERS={}".format(lyr)
 wms += "&STYLES={}&".format(style)
 wms += "SRS=EPSG:{}&".format(epsg)
 wms += "BBOX={},{},{},{}&".format(minx, miny, maxx, maxy)
 wms += "WIDTH={}&".format(w)
 wms += "HEIGHT={}&".format(h)
 wms += "FORMAT=image/jpeg"
 wmsmap = urllib.request.urlopen(wms)
 with open(img + ".jpg", "wb") as f:
   f.write(wmsmap.read())
  1. 接下来，我们有一个exif()函数从照片中提取元数据:
def exif(img):
 """Return EXIF metatdata from image"""
 exif_data = {}
 try: 
 i = Image.open(img)
 tags = i._getexif()
 for tag, value in tags.items():
 decoded = TAGS.get(tag, tag)
 exif_data[decoded] = value
 except:
 pass
 return exif_data
  1. 然后我们有一个dms2dd()函数将度/分/秒坐标转换为十进制度，因为照片坐标就是这样存储的:
def dms2dd(d, m, s, i):
 """Convert degrees/minutes/seconds to
 decimal degrees"""
 s *= .01
 sec = float((m * 60.0) + s)
 dec = float(sec / 3600.0)
 deg = float(d + dec)
 if i.upper() == 'W':
 deg = deg * -1.0
 elif i.upper() == 'S':
 deg = deg * -1.0
 return float(deg)
  1. 最后，我们有一个gps()函数从照片元数据中提取坐标:
def gps(exif):
 """Extract GPS info from EXIF metadat"""
 lat = None
 lon = None
 if exif['GPSInfo']: 
 # Lat
 coords = exif['GPSInfo']
 i = coords[1]
 d = coords[2][0][0]
 m = coords[2][1][0]
 s = coords[2][2][0]
 lat = dms2dd(d, m ,s, i)
 # Lon
 i = coords[3]
 d = coords[4][0][0]
 m = coords[4][1][0]
 s = coords[4][2][0]
 lon = dms2dd(d, m ,s, i)
 return lat, lon
  1. 接下来，我们有我们的程序变量。我们将访问由一家名为的公司免费提供的WMS 开放街道地图服务，以及美国宇航局提供的 SRTM 数据。

为了简单起见，我们使用 Python 的urllib库来访问本书中的 WMS 服务，但是如果您计划频繁使用 OGC web 服务，您应该使用通过 PyPI 提供的 Python 包 OWSLib:pypi.python.org/pypi/OWSLib。

现在，让我们执行以下步骤来设置 WMS web 服务:

  1. 我们将输出几个中间产品和图像。这些变量在这些步骤中使用。route.gpx文件在本节中定义为gpx变量。首先，我们用下面的代码为度数到弧度的转换设置一些转换常数:
# Needed for numpy conversions in hillshading
deg2rad = 3.141592653589793 / 180.0
rad2deg = 180.0 / 3.141592653589793
  1. 接下来，我们将.gpx文件的名称设置如下:
# Program Variables

# Name of the gpx file containing a route.
# https://git.io/fjwHW
gpx = "route.gpx"
  1. 现在，我们开始设置 WMS 网络服务，它将检索地图:
# NOAA OpenStreetMap Basemap

# OSM WMS service
osm_WMS = "http://ows.mundialis.de/services/service"

# Name of the WMS street layer
# streets = "osm"
osm_lyr = "OSM-WMS"

# Name of the basemap image to save
osm_img = "basemap"

# OSM EPSG code (spatial reference system)
osm_epsg = 3857

# Optional WMS parameter
osm_style = ""
  1. 接下来，我们设置山体阴影参数，它将决定人造太阳的角度和方向:
# Shaded elevation parameters
#
# Sun direction
azimuth = 315.0

# Sun angle
altitude = 45.0

# Elevation exageration
z = 5.0

# Resolution
scale = 1.0
  1. 然后我们在没有高程信息的地方设置no_data值:
# No data value for output
no_data = 0
  1. 接下来，我们将输出图像的名称设置如下:
# Output elevation image name
elv_img = "elevation"
  1. 现在，我们使用以下代码为最小和最大高程值创建颜色:
# RGBA color of the SRTM minimum elevation
min_clr = (255, 255, 255, 0)

# RGBA color of the SRTM maximum elevation
max_clr = (0, 0, 0, 0)

# No data color
zero_clr = (255, 255, 255, 255)
  1. 然后我们设置输出图像大小，如下所示:
# Pixel width and height of the

# output images
w = 800
h = 800

现在我们理解了函数是如何工作的，让我们解析一下 GPX。

解析 GPX

现在，我们将使用built-in xml.dom.minidom模块解析 GPX 文件，它只是 XML。我们将提取纬度、经度、海拔和时间戳。我们将把它们存储在一个列表中以备后用。使用 Python 的time模块将时间戳转换为struct_time对象，这使得处理更加容易。

解析需要执行以下步骤:

  1. 首先，我们使用minidom模块解析gpx文件:
# Parse the gpx file and extract the coordinates
log.info("Parsing GPX file: {}".format(gpx))
xml = minidom.parse(gpx)
  1. 接下来，我们得到所有包含高程信息的"trkpt"标签:
# Grab all of the "trkpt" elements
trkpts = xml.getElementsByTagName("trkpt")
  1. 现在，我们设置列表来存储解析后的位置和高程值:
# Latitude list
lats = []
# Longitude list
lons = []
# Elevation list
elvs = []
# GPX timestamp list
times = []
  1. 然后，我们遍历 GPX 的全球定位系统条目并解析这些值:
# Parse lat/long, elevation and times
for trkpt in trkpts:
 # Latitude
 lat = float(trkpt.attributes["lat"].value)
 # Longitude
 lon = float(trkpt.attributes["lon"].value)
 lats.append(lat)
 lons.append(lon)
 # Elevation
 elv = trkpt.childNodes[0].firstChild.nodeValue
 elv = float(elv)
 elvs.append(elv)

时间戳需要一点额外的工作，因为我们必须从格林尼治时间转换为当地时间:

 # Times
 t = trkpt.childNodes[1].firstChild.nodeValue
 # Convert to local time epoch seconds
 t = get_local_time(t)
 times.append(t)

解析 GPX 之后，我们需要路径的边界框来从其他地理空间服务下载数据。

获取边界框

当我们下载数据时，我们希望数据集比路线覆盖更多的区域，这样地图就不会在路线边缘附近被裁剪得太紧。所以我们将缓冲边界框的每一边 20%。最后，我们需要东部和北部的数据来配合 WMS 服务。东距和北距是笛卡尔坐标系中点的 x 和 y 坐标，单位为米。它们通常用于 UTM 坐标系:

  1. 首先，我们从坐标列表中获取范围，如下所示:
# Find Lat/Long bounding box of the route
minx = min(lons)
miny = min(lats)
maxx = max(lons)
maxy = max(lats)
  1. 接下来，我们缓冲边界框，以确保轨迹不会靠近边缘:
# Buffer the GPX bounding box by 20%
# so the track isn't too close to
# the edge of the image.
xdist = maxx - minx
ydist = maxy - miny
x20 = xdist * .2
y20 = ydist * .2

# 10% expansion on each side
minx -= x20
miny -= y20
maxx += x20
maxy += y20
  1. 最后，我们在一个变量中设置我们的边界框，并将我们的坐标转换为米，这是 web 服务所需要的:
# Store the bounding box in a single
# variable to streamline function calls
bbox = [minx, miny, maxx, maxy]

# We need the bounding box in meters
# for the OSM WMS service. We will
# download it in degrees though to
# match the SRTM file. The WMS spec
# says the input SRS should match the
# output but this custom service just
# doesn't work that way
mminx, mminy = ll2m(miny, minx)
mmaxx, mmaxy = ll2m(maxy, maxx)

有了这个，我们现在将下载我们的地图和海拔图像。

下载地图和高程图像

我们将首先下载 OSM 底图作为我们的底图，该底图包含街道和标签:

  1. 首先，我们将使用log.info下载 OSM 底图:
# Download the OSM basemap
log.info("Downloading basemap")
wms(mminx, mminy, mmaxx, mmaxy, osm_WMS, osm_lyr,
 osm_epsg, osm_style, osm_img, w, h)

该部分将生成一个中间图像，如下图所示:

  1. 接下来，我们将从 SRTM 数据集下载一些高程数据。SRTM 几乎是全球性的，提供 30-90 米的分辨率。SRTM.py Python 模块使处理这些数据变得容易。SRTM.py下载数据并设置其需要进行请求。因此，如果您从不同的区域下载数据，您可能需要清理位于主目录(~/.srtm)中的缓存。脚本的这一部分可能需要 2-3 分钟才能完成，具体取决于您的计算机和互联网连接速度:
# Download the SRTM image

# srtm.py downloader
log.info("Retrieving SRTM elevation data")
# The SRTM module will try to use a local cache

# first and if needed download it.
srt = srtm.get_data()
# Get the image and return a PIL Image object
image = srt.get_image((w, h), (miny, maxy), (minx, maxx),
 300, zero_color=zero_clr, min_color=min_clr,
 max_color=max_clr)
# Save the image
image.save(elv_img + ".png")

脚本的这一部分还输出了一个中间高程图像，如下图所示:

现在我们有了海拔图像，我们可以把它变成山体阴影。

创建山体阴影

我们可以通过在第 7 章、 Python 和高程数据中创建阴影浮雕*部分使用的相同山体阴影算法来运行该数据。*为此，让我们遵循以下步骤:

  1. 首先，我们打开我们的高程图像，并将其读入一个numpy数组:
# Hillshade the elevation image
log.info("Hillshading elevation data")
im = Image.open(elv_img + ".png").convert("L")
dem = np.asarray(im)
  1. 现在，我们设置我们的处理窗口在网格中移动，并在小部分中进行分析以提高效率:
# Set up structure for a 3x3 windows to
# process the slope throughout the grid
window = []
# x, y resolutions
xres = (maxx-minx)/w
yres = (maxy-miny)/h
  1. 然后，我们将立面图像分成如下窗口:
# Create the windows
for row in range(3):
 for col in range(3):
 window.append(dem[row:(row + dem.shape[0]-2),
 col:(col + dem.shape[1]-2)])
  1. 我们将为我们的处理窗口创建数组，如下所示:
# Process each cell
x = ((z * window[0] + z * window[3] + z * window[3] + z * window[6]) -
 (z * window[2] + z * window[5] + z * window[5] + z * window[8])) \
 / (8.0 * xres * scale)

y = ((z * window[6] + z * window[7] + z * window[7] + z * window[8]) -
 (z * window[0] + z * window[1] + z * window[1] + z * window[2])) \
 / (8.0 * yres * scale)
  1. 最后，我们可以通过numpy一次性处理它们:
# Calculate slope
slope = 90.0 - np.arctan(np.sqrt(x*x + y*y)) * rad2deg

# Calculate aspect
aspect = np.arctan2(x, y)

# Calculate the shaded relief
shaded = np.sin(altitude * deg2rad) * np.sin(slope * deg2rad) \
 + np.cos(altitude * deg2rad) * np.cos(slope * deg2rad) \
 * np.cos((azimuth - 90.0) * deg2rad - aspect)

shaded = shaded * 255

现在我们有了山体阴影图层，我们可以开始创建地图了。

创建地图

我们有了开始构建报告地图所需的数据。我们的方法如下:

  • 使用过滤器增强高程和底图图像
  • 将图像混合在一起，以提供带有山坡的 OSM 地图
  • 创建半透明层以绘制街道路线
  • 将路线图层与山坡地图混合

这些任务都将使用 PIL Image和ImageDraw模块来完成，如以下步骤所示:

  1. 首先，我们将我们的阴影浮雕numpy数组转换回图像并平滑它:
# Convert the numpy array back to an image
relief = Image.fromarray(shaded).convert("L")

# Smooth the image several times so it's not pixelated
for i in range(10):
 relief = relief.filter(ImageFilter.SMOOTH_MORE)

log.info("Creating map image")
  1. 现在，我们将增加图像的对比度，使其更加突出:
# Increase the hillshade contrast to make
# it stand out more
e = ImageEnhance.Contrast(relief)
relief = e.enhance(2)
  1. 接下来，我们将地图图像裁剪为与高程图像相同的大小:
# Crop the image to match the SRTM image. We lose
# 2 pixels during the hillshade process
base = Image.open(osm_img + ".jpg").crop((0, 0, w-2, h-2))
  1. 然后，我们还增加了地图图像的对比度，并将其与山体阴影图像混合:
# Enhance basemap contrast before blending
e = ImageEnhance.Contrast(base)
base = e.enhance(1)

# Blend the the map and hillshade at 90% opacity
topo = Image.blend(relief.convert("RGB"), base, .9)
  1. 现在，我们准备好在混合地图上绘制全球定位系统轨迹，首先将我们的点转换为像素:
# Draw the GPX tracks
# Convert the coordinates to pixels
points = []
for x, y in zip(lons, lats):
 px, py = world2pixel(x, y, w, h, bbox)
 points.append((px, py))
  1. 我们还需要从将要创建的轨迹图像的边缘缓冲区中减去缓冲区:
# Crop the image size values to match the map
w -= 2
h -= 2
  1. 接下来，我们创建一个透明图像，并将我们的轨迹绘制为红线:
# Set up a translucent image to draw the route.
# This technique allows us to see the streets
# and street names under the route line.

track = Image.new('RGBA', (w, h))

track_draw = ImageDraw.Draw(track)

# Route line will be red at 50% transparency (255/2=127)
track_draw.line(points, fill=(255, 0, 0, 127), width=4)
  1. 现在，我们可以使用以下代码将轨道粘贴到图像上:
# Paste onto the basemap using the drawing layer itself
# as a mask.
topo.paste(track, mask=track)
  1. 现在我们将在路线上画一个起点，如下所示:
# Now we'll draw start and end points directly on top
# of our map - no need for transparency
topo_draw = ImageDraw.Draw(topo)

# Starting circle
start_lon, start_lat = (lons[0], lats[0])
start_x, start_y = world2pixel(start_lon, start_lat, w, h, bbox)
start_point = [start_x-10, start_y-10, start_x+10, start_y+10]
topo_draw.ellipse(start_point, fill="lightgreen", outline="black")
start_marker = [start_x-4, start_y-4, start_x+4, start_y+4]
topo_draw.ellipse(start_marker, fill="black", outline="white")
  1. 下面是结束点的代码片段:
# Ending circle
end_lon, end_lat = (lons[-1], lats[-1])
end_x, end_y = world2pixel(end_lon, end_lat, w, h, bbox)
end_point = [end_x-10, end_y-10, end_x+10, end_y+10]
topo_draw.ellipse(end_point, fill="red", outline="black")
end_marker = [end_x-4, end_y-4, end_x+4, end_y+4]
topo_draw.ellipse(end_marker, fill="black", outline="white")

现在我们已经画出了我们的轨迹，我们准备放置我们的地理标记的照片。

定位照片

我们将使用添加了全球定位系统位置坐标的手机拍摄的照片。你可以从 下载。

将图像放在与脚本相同级别的名为photos的目录中。我们将只使用一张照片，但脚本可以处理你想要的男人图像。我们将在地图上绘制并放置一个照片图标，然后保存完成的底图，如以下步骤所示:

  1. 首先，我们获得一个包含以下代码的图像列表:
# Photo icon
images = glob.glob("photos/*.jpg")
  1. 接下来，我们遍历每张图像并获取其全球定位系统信息:
for i in images:
 e = exif(i)
  1. 然后，我们使用全球定位系统函数解析位置信息，如下所示:
 photo_lat, photo_lon = gps(e)
 #photo_lat, photo_lon = 30.311364, -89.324786
  1. 现在，我们可以将照片坐标转换为图像像素坐标:
 photo_x, photo_y = world2pixel(photo_lon, photo_lat, w, h, bbox)
  1. 然后，我们将使用以下代码为照片的位置绘制一个图标:
 topo_draw.rectangle([photo_x - 12, photo_y - 10, photo_x + 12, \
 photo_y + 10], fill="black", outline="black")
 topo_draw.rectangle([photo_x - 9, photo_y - 8, photo_x + 9, \
 photo_y + 8], fill="white", outline="white")
 topo_draw.polygon([(photo_x-8,photo_y+7), (photo_x-3,photo_y-1), (photo_x+2,photo_y+7)], fill = "black")
 topo_draw.polygon([(photo_x+2,photo_y+7), (photo_x+7,photo_y+3), (photo_x+8,photo_y+7)], fill = "black")
  1. 最后，我们将像这样保存我们的地图:
# Save the topo map
topo.save("{}_topo.jpg".format(osm_img))

虽然没有保存到文件系统中，但山坡高程如下所示:

混合地形图如下图所示:

虽然山体阴影制图让我们了解了海拔高度，但它没有给我们任何定量数据。为了获得更详细的信息，我们将创建一个简单的立面图。

测量高度

使用出色的谷歌地图应用编程接口，我们可以快速构建一个漂亮的高程剖面图，显示路线上的高程变化:

  1. 首先，我们将为我们的高程剖面创建chart对象:
# Build the elevation chart using the Google Charts API
log.info("Creating elevation profile chart")
chart = SimpleLineChart(600, 300, y_range=[min(elvs), max(elvs)])
  1. 现在，我们需要为最小值创建一条线，如下所示:
# API quirk - you need 3 lines of data to color
# in the plot so we add a line at the minimum value
# twice.
chart.add_data([min(elvs)]*2)
chart.add_data(elvs)
chart.add_data([min(elvs)]*2)

# Black lines
chart.set_colours(['000000'])
  1. 接下来，我们可以按如下方式填写我们的高程剖面:
# fill in the elevation area with a hex color
chart.add_fill_range('80C65A', 1, 2)
  1. 然后，我们可以按如下方式设置高程标签，并将它们指定给一个轴:
# Set up labels for the minimum elevation, halfway value, and max value
elv_labels = int(round(min(elvs))), int(min(elvs)+((max(elvs)-min(elvs))/2))

# Assign the labels to an axis
elv_label = chart.set_axis_labels(Axis.LEFT, elv_labels)
  1. 接下来，我们可以用以下代码标记轴本身:
# Label the axis
elv_text = chart.set_axis_labels(Axis.LEFT, ["FEET"])
# Place the label at 30% the distance of the line
chart.set_axis_positions(elv_text, [30])
  1. 现在我们可以计算轨迹点之间的距离:
# Calculate distances between track segments
distances = []
measurements = []
coords = list(zip(lons, lats))
for i in range(len(coords)-1):
 x1, y1 = coords[i]
 x2, y2 = coords[i+1]
 d = haversine(x1, y1, x2, y2)
 distances.append(d)
total = sum(distances)
distances.append(0)
j = -1

我们有高程剖面，但是我们需要沿着 x 轴添加距离标记，这样我们就知道剖面沿着路线在哪里发生了变化。

测量距离

为了理解高程数据图，我们需要沿着 x 轴的参考点来帮助我们确定路线沿线的高程。我们将计算路线沿线的英里分割，并将其放置在图表 x 轴上的适当位置。让我们看看以下步骤:

  1. 首先，我们沿着轴定位英里标记，如下所示:
# Locate the mile markers
for i in range(1, int(round(total))):
 mile = 0
 while mile < i:
 j += 1
 mile += distances[j]
 measurements.append((int(mile), j))
 j = -1
  1. 接下来，我们为英里标记设置标签:
# Set up labels for the mile points

positions = []
miles = []
for m, i in measurements:
 pos = ((i*1.0)/len(elvs)) * 100
 positions.append(pos)
 miles.append(m)

# Position the mile marker labels along the x axis
miles_label = chart.set_axis_labels(Axis.BOTTOM, miles)
chart.set_axis_positions(miles_label, positions)
  1. 现在，我们可以将英里标记标记如下:
# Label the x axis as "Miles"
miles_text = chart.set_axis_labels(Axis.BOTTOM, ["MILES", ])
chart.set_axis_positions(miles_text, [50, ])

# Save the chart
chart.download('{}_profile.png'.format(elv_img))

我们的图表现在应该如下图所示:

我们的第一张图表完成了。现在，让我们看看沿途的天气数据。

检索天气数据

在本节中，我们将检索我们的最终数据元素:天气。如前所述，我们将使用黑暗天空服务，该服务允许我们收集世界上任何地方的历史天气报告。天气 API 是基于 REST 和 JSON 的，所以我们将使用urllib模块请求数据，使用json库解析数据。本节中值得注意的是，我们在本地缓存数据，因此如果需要，您可以离线运行脚本进行测试。在这一部分的早期，您可以放置由YOUR KEY HERE文本标记的黑暗天空应用编程接口密钥。让我们看看以下步骤:

  1. 首先，我们需要我们感兴趣区域的中心:
log.info("Creating weather summary")

# Get the bounding box centroid for georeferencing weather data
centx = minx + ((maxx-minx)/2)
centy = miny + ((maxy-miny)/2)
  1. 现在，我们按如下方式设置免费的黑暗应用编程接口密钥，这样我们就可以检索天气数据:
# DarkSky API key
# You must register for free at DarkSky.net
# to get a key to insert here.
api_key = "YOUR API KEY GOES HERE"
  1. 然后，我们从数据中获取最新的时间戳，用于天气查询:
# Grab the latest route time stamp to query weather history
t = times[-1]
  1. 现在，我们准备进行天气数据查询，如下所示:
history_req = "https://api.darksky.net/forecast/{}/".format(api_key)
#name_info = [t.tm_year, t.tm_mon, t.tm_mday, route_url.split(".")[0]]
#history_req += "/history_{0}{1:02d}{2:02d}/q/{3}.json" .format(*name_info)
history_req += "{},{},{}".format(centy, centx, t)
history_req += "?exclude=currently,minutely,hourly,alerts,flags"
request = urllib.request.urlopen(history_req)
weather_data = request.read()
  1. 我们将像这样缓存天气数据，以防以后需要查看:
# Cache weather data for testing
with open("weather.json", "w") as f:
 f.write(weather_data.decode("utf-8"))
  1. 然后，我们解析天气 JSON 数据如下:
# Retrieve weather data
js = json.loads(open("weather.json").read())
history = js["daily"]
  1. 我们只需要天气摘要，这是列表中的第一项:
# Grab the weather summary data.
# First item in a list.
daily = history["data"][0]
  1. 现在，我们将获得具体的天气属性如下:
# Max temperature in Imperial units (Farenheit).
# Celsius would be metric: "maxtempm"
maxtemp = daily["temperatureMax"]

# Minimum temperature
mintemp = daily["temperatureMin"]

# Maximum humidity
maxhum = daily["humidity"]

# Precipitation in inches (cm = precipm)
if "precipAccumulation" in daily:
 precip = daily["precipAccumulation"]
else:
 precip = "0.0"
  1. 现在我们已经将天气数据存储在变量中，我们可以完成最后一步:将其全部添加到 PDF 报告中。

在某些情况下，除了 PIL 之外，库没有依赖关系。就我们的目的而言，它会很好地工作。我们将继续下一页并添加元素。fpdf.ln()分隔行，而fpdf.cells包含文本并允许更精确的布局。

我们终于可以通过以下步骤创建 PDF 报告了:

  1. 首先，我们如下设置我们的pdf对象:
# Simple fpdf.py library for our report.
# New pdf, portrait mode, inches, letter size
# (8.5 in. x 11 in.)
pdf = fpdf.FPDF("P", "in", "Letter")
  1. 然后，我们将为报告添加一个页面，并设置我们的字体首选项:
# Add our one report page
pdf.add_page()

# Set up the title
pdf.set_font('Arial', 'B', 20)
  1. 我们将用下面的代码为我们的报告创建一个标题:
# Cells contain text or space items horizontally
pdf.cell(6.25, 1, 'GPX Report', border=0, align="C")

# Lines space items vertically (units are in inches)
pdf.ln(h=1)
pdf.cell(1.75)

# Create a horizontal rule line
pdf.cell(4, border="T")
pdf.ln(h=0)
pdf.set_font('Arial', style='B', size=14)
  1. 现在，我们可以像这样添加路线图:
# Set up the route map
pdf.cell(w=1.2, h=1, txt="Route Map", border=0, align="C")
pdf.image("{}_topo.jpg".format(osm_img), 1, 2, 4, 4)
pdf.ln(h=4.35)
  1. 接下来，我们添加如下立面图:
# Add the elevation chart
pdf.set_font('Arial', style='B', size=14)
pdf.cell(w=1.2, h=1, txt="Elevation Profile", border=0, align="C")
pdf.image("{}_profile.png".format(elv_img), 1, 6.5, 4, 2)
pdf.ln(h=2.4)
  1. 然后，我们可以用以下代码编写天气数据摘要:
# Write the weather summary
pdf.set_font('Arial', style='B', size=14)
pdf.cell(1.2, 1, "Weather Summary", align="C")
pdf.ln(h=.25)
pdf.set_font('Arial', style='', size=12)
pdf.cell(1.8, 1, "Min. Temp.: {}".format(mintemp), align="L")
pdf.cell(1.2, 1, "Max. Hum.: {}".format(maxhum), align="L")
pdf.ln(h=.25)
pdf.cell(1.8, 1, "Max. Temp.: {}".format(maxtemp), align="L")
pdf.cell(1.2, 1, "Precip.: {}".format(precip), align="L")
pdf.ln(h=.25)
  1. “黑暗天空”条款要求我们在报告中添加一个徽标，这归功于出色的数据源:
# Give Dark Sky credit for a great service (https://git.io/fjwHl)
pdf.image("darksky.png", 3.3, 9, 1.75, .25)
  1. 现在，我们可以使用以下代码添加地理定位图像:
# Add the images for any geolocated photos
pdf.ln(h=2.4)
pdf.set_font('Arial', style='B', size=14)
pdf.cell(1.2, 1, "Photos", align="C")
pdf.ln(h=.25)
for i in images:
 pdf.image(i, 1.2, 1, 3, 3)
 pdf.ln(h=.25)
  1. 最后，我们可以保存报告并查看它:
# Save the report
log.info("Saving report pdf")
pdf.output('report.pdf', 'F')

您的工作目录中应该有一个名为report.pdf的 PDF 文档，其中包含您的成品。它应该看起来像下面截图中显示的图像:

有了这个，我们使用了我们在这本书中学到的所有技术，并建立了一个全球定位系统报告工具。

摘要

恭喜你！在这本书里，你汇集了成为一名现代地理空间分析师所需的最基本的工具和技能。无论您是偶尔使用地理空间数据还是一直使用它，您都将能够更好地充分利用地理空间分析。这本书的重点是使用几乎完全在 PyPI 目录中找到的开源工具，以便于安装和集成。但是，即使您使用 Python 作为商业地理信息系统包或流行库(如 GDAL)的驱动程序，用纯 Python 测试新概念的能力也总是会派上用场。

进一步阅读

Python 为可视化数据提供了一套丰富的库。其中最突出的是 Matplotlib ，可以生成无数类型的图表和地图，并保存为 PDF。Packt 有几本关于 Matplotlib 的书，包括 Matplotlib 30 食谱:https://www . packtpub . com/大数据与商业智能/matplotlib-30 食谱。

第一部分：行业的历史和现状

本节首先使用插图、基本公式、简单代码和 Python 演示常见的地理空间分析过程。在此基础上，您将学习如何使用地理空间数据——获取数据并为各种分析做准备。之后，您将了解地理空间技术生态系统中使用的各种软件包和库。在本节的最后，您将学习如何评估任何地理空间工具。

本节包括以下章节:

第二部分：地理空间分析概念

本节代表了本书的主要构建模块，在这里您将通过不同的代码示例和数据编辑概念了解 Python 在行业中的作用。您将了解地理空间产品以及如何应用它们来解决问题。接下来，您将看到如何使用 Python 实际处理遥感数据。在本节的最后，您将了解如何以任何地理空间格式使用高程数据来分析三维要素。

本节包括以下章节:

第三部分：实用地理空间处理技术

这部分是高级水平，它将需要你以前学习的所有技能。它从学习如何创建地理空间模型来回答特定的问题开始。此外，它将向您展示一些构建地理空间模型的技术，以及如何使用可视化概念帮助预测未来。我们将继续访问和处理实时数据。在这一节的最后，我们将结合前面几节所学的知识，实现一个系统，根据全球定位系统数据和地理标记的照片创建一个户外跑步或徒步报告。

本节包括以下章节:

avatar
文章
阅读
粉丝