In [1]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
from numpy import linspace
from numpy import meshgrid
from matplotlib.colors import LinearSegmentedColormap
# from helper import get_tmp_file
cmap1 = LinearSegmentedColormap.from_list("my_colormap", ((0, 0, 0), (1, 1, 1)), N=6, gamma=1.0)
- 创建contourf绘图的方法是从contourf示例中获取的
- 使用LinearSegmentedColormap.from_list静态方法创建颜色映射。此方法具有以下参数：
- 1 给予颜色图名称
- 2 颜色列表。这是一个列表或序列，每个元素包含0到1范围内的三个浮点，它们是颜色的红色，绿色和蓝色值
- 3 N是要创建的颜色级别的数量。如果数字低于列表的颜色长度，列表将被截断。如果它的时间更长，一些颜色会重复
In [2]:
map = Basemap(projection='tmerc',
lat_0=0, lon_0=3,
llcrnrlon=1.819757266426611,
llcrnrlat=41.583851612359275,
urcrnrlon=1.841589961763497,
urcrnrlat=41.598674173123)
ds = gdal.Open("/gdata/sample_files/dem.tiff")
data = ds.ReadAsArray()
x = linspace(0, map.urcrnrx, data.shape[1])
y = linspace(0, map.urcrnry, data.shape[0])
xx, yy = meshgrid(x, y)
map.contourf(xx, yy, data, (400, 600, 800, 1000, 1200), cmap=cmap1)
plt.show()
- 在我们的例子中，六个级别值将从黑色到白色。
- 轮廓线被创建以强制在序列的每个元素中的级别。
- GDAL包括一个名为gdaldem的实用程序，它使文件中定义的值对栅格进行分类。文件格式最初是由GRASS r.colors函数使用的。格式非常好，因为是真的很容易理解，可以从不同的软件使用。这里是如何阅读它，并与Basemap一起使用：
In [3]:
# from osgeo import gdal
# from numpy import linspace
# from numpy import meshgrid
from os.path import exists
# from matplotlib.colors import LinearSegmentedColormap
def read_color_table(color_file):
'''
The method for reading the color file.
'''
colors = []
levels = []
if exists(color_file) is False:
raise Exception("Color file " + color_file + " does not exist")
fp = open(color_file, "r")
for line in fp:
if line.find('#') == -1 and line.find('/') == -1:
entry = line.split()
levels.append(eval(entry[0]))
colors.append((int(entry[1]) / 255., int(entry[2]) / 255., int(entry[3]) / 255.))
fp.close()
cmap = LinearSegmentedColormap.from_list("my_colormap", colors, N=len(levels), gamma=1.0)
return levels, cmap
levels, cmap = read_color_table("/gdata/sample_files/colorfile.clr")
函数read_color_table打开并读取颜色文件，返回文件中定义的级别，以及具有其颜色的颜色图
检查某行是否已注释
原始值在0-255范围内，必须转换为0-1
使用LinearSegmentedColormap.from_list静态方法创建颜色映射
In [4]:
map = Basemap(projection='tmerc',
lat_0=0, lon_0=3,
llcrnrlon=1.819757266426611,
llcrnrlat=41.583851612359275,
urcrnrlon=1.841589961763497,
urcrnrlat=41.598674173123)
ds = gdal.Open("/gdata/sample_files/dem.tiff")
data = ds.ReadAsArray()
x = linspace(0, map.urcrnrx, data.shape[1])
y = linspace(0, map.urcrnry, data.shape[0])
xx, yy = meshgrid(x, y)
map.contourf(xx, yy, data, levels, cmap=cmap)
plt.show()
- 该函数可以返回级别和颜色映射，并且它们可以一起使用或单独使用
- contourf与文件中定义的值一起使用
In [5]:
fig = plt.figure()
# 在调用Basemap构造函数之前，将调用fig.add_subplot方法。这三个数字是：
# 最终图中的行数
# 最终图中的列数
# *使用哪个轴（子图），从左上角的一个轴计数，如StackOverflow问题所述
# 一旦创建了轴，以后创建的映射将自动使用它（虽然可以传递ax参数以使用所选择的轴）
# 可以使用set_title（）将标题添加到每个子图中
ax = fig.add_subplot(211)
ax.set_title("Hammer projection")
map = Basemap(projection='hammer', lon_0 = 10, lat_0 = 50)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
ax = fig.add_subplot(212)
ax.set_title("Robinson projection")
map = Basemap(projection='robin', lon_0 = 10, lat_0 = 50)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
plt.show()
In [6]:
fig, axes = plt.subplots(2, 1)
axes[0].set_title("Hammer projection")
map = Basemap(projection='hammer', lon_0 = 10, lat_0 = 50, ax=axes[0])
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
axes[1].set_title("Robinson projection")
map = Basemap(projection='robin', lon_0 = 10, lat_0 = 50, ax=axes[1])
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
plt.show()
- 传递给子图的参数是要创建行数和列数的。
- 子图方法返回figure对象，以及创建的轴（子图）的列表，其中第一个元素是左上角位置的元素。
- 在创建Basemap实例时，必须使用创建的轴传递ax参数。
- 结果与前面的例子相同。
- 在示例中，布局用于显示不同的缩放级别。每个级别都用上一个地图上的多边形表示。要创建这些指标需注意以下几点：
- 计算边界框的每个角的位置。
- 可以使用Basemap实例的xmin，xmax，ymin，ymax字段检索拐角（下一个地图，因为这将指示缩放区域）。
- 由于字段是投影单位，因此使用带有反向参数的Basemap实例。请参阅使用Basemap实例转换单位部分，看看它是如何工作的。
- 将每个角的计算的经度和纬度传递到当前地图实例以获得当前地图投影中的坐标。因为我们有一个序列中的值，所以使用星号来解包它们。
- 一旦知道点，就使用Path类创建多边形。
In [7]:
import helper; helper.info()