Python遥感可视化 — Basemap作Albers投影

欢迎关注博主的微信公众号：“智能遥感”。

该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。

本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

---

“ Python + Basemap对遥感地图可视化。”

今天的遥感之美封面图—圣克莱尔湖“绽放”的藻类。澳大利亚塔斯马尼亚州中西部湖泊，位于克雷德尔山-圣克莱尔湖国家公园南部边缘，面积为28平方公里。原本为山谷，因冰川刨蚀而加深，后为冰川堰塞湖。为澳大利亚最深的湖，西部最深点约215米。

2015年7月28日，Landsat 8卫星上的陆地成像仪（OLI）拍摄了大湖周围爆发的藻类图像，在圣克莱尔湖和西伊利湖的图像中可见绿色“漩涡”。

7月早些时候，美国国家海洋和大气管理局（NOAA）的科学家们预测，在伊利湖西部，2015年有害藻类大量繁殖将会非常严重。他们认为，伊利湖西部的藻类生长可以与2011年的藻类爆发事件相媲美。当有大量营养物质（许多来自农业径流）和阳光以及适宜水温时，该盆地中的藻类会茁壮成长。整个夏季为藻类生长期，九月达到高峰。

研究证实，2011年来自农田径流的磷与有利的天气、湖水条件相结合，产生的爆发比以前观察到的面积大三倍。研究人员指出，如果土地开采活动和气候变化趋势持续下去，可能会看到比2011年更大面积的爆发。

有害藻类的大量繁殖会导致鱼类死亡，它们还可以影响用于娱乐和消费的水体安全。截至2015年7月30日，据报道这些地区的饮用水是安全的。

2015年4月，美国宇航局（NASA）联合几个单位开发了一项新的用于早期预警淡水中有害藻类大量繁殖的系统。该系统有望使环境和水质管理者更容易获得海洋彩色卫星图像（译自LandSat官网）。

---

Albers投影

遥感数据预处理部分到此结束，下面给大家分享遥感地图可视化。遥感可视化主要基于Basemap模块实现，不同的投影会使得遥感图像/矢量地图的形状不一样，这是给人最直观的感受。不同的投影适用于不同的目的，例如网络映射，精确测量距离或计算面积。我国常用的投影如墨卡托投影、高斯-克吕格投影以及UTM投影。

为了在二维地图上表示三维球面，需要对地图进行投影。由于地图投影均有失真情况，因此有许多种地图投影，每个投影都有自己的优缺点。Basemap模块提供了24种不同的地图投影。有些是全球范围，有些只代表全球某个地方。创建类Basemap实例时，必须指定所需的地图投影，以及有关地图投影将描述的地球表面部分信息。这有两种基本方法：一种是提供矩形地图投影区域的四角点的经纬度信息；另一种是提供地图投影区域的中心经纬度以及地图投影坐标中的区域宽度和高度。下面给大家展示的是Albers投影，主要代码参考Basemap官网。

阿尔伯斯投影，又名“正轴等积割圆锥投影”，“双标准纬线等积圆锥投影”。圆锥投影的一种，为阿伯斯（Albers）拟定。在等面积投影上，面积保持不变但形状发生变化。具体代码实现如下：

 

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2019/1/9 11:40'

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
import matplotlib

# 中文字体设置,防止中文乱码
# 指定默认字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']

# width为投影坐标中所需地图域的宽度(米)
# height投影坐标中所需地图域的高度(米)
# l:low低分辨率; aea:Albers Equal Area
# lat_1为第一标准线
# lat_2为第二标准线
# lon_0,lat_0为地图中心点
mp = Basemap(width=8000000, height=7000000,
            resolution='l', projection='aea',
            lat_1=40., lat_2=60., lon_0=100, lat_0=35)
mp.drawcoastlines()
mp.drawcountries()
# 内陆湖泊采用蓝色
mp.fillcontinents(color='red', lake_color='b')

# [1, 1, 0, 0]左右纬度显示, [0, 0, 0, 1]只显示下方经度
mp.drawparallels(np.arange(-80., 81, 20.), labels=[1, 1, 0, 0], fontsize=12)
mp.drawmeridians(np.arange(-180., 181., 20.), labels=[0, 0, 0, 1], fontsize=12)
# 海洋采用水青色
mp.drawmapboundary(fill_color='aqua')

# 获取当前坐标轴对象
ax = plt.gca()
# 添加地图标题
plt.title(u"Albers等面积投影")
# 将文本叠加到地图上,但是Basemap没有annotate方法,所以使用pyplot接口
# xy=(0.42, 0.47)表示本文在图中位置
# xycoords:坐标轴比例
# color:字体颜色
# fontsize:字体大小
plt.annotate(u"中华人民共和国", xy=(0.42, 0.47), xycoords='axes fraction',
             color='k', fontsize=12)
plt.show()

              上述代码运行结果如下：