如何使等高线图叠加在底图上 [英] How can I get my contour plot superimposed on a basemap

查看:122
本文介绍了如何使等高线图叠加在底图上的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

这是我几个月前问过的一个问题,但仍在努力寻求解决方案.我的代码并排提供了底图和轮廓图(但打印到文件仅给出轮廓图),但是我希望它们叠加在一起.最好的解决方案是方法,它们既更快又更灵活:

来自scipy.interpolate的

 以gd格式导入griddatazi = gd((数据[['projected_lon','projected_lat']]),data.Z.值,(xi,yi),method ='linear') 

如果您使用scipy的 griddata 方法,则还必须确定哪种方法( nearest linear 立方)给出最佳的结果图.

我还要补充一点,上面演示和讨论的插值方法是最简单的方法,不一定对降雨数据的插值有效.本文很好地概述了有效的方法和注意事项用于水文和水文建模.这些的实现(可能使用Scipy)留作练习& c.

This is a question I asked several months ago and am still struggling to come to a solution. My code gives me a basemap and a contour plot side by side (but printing to file only gives the contour plot), but I want them superimposed. The best solution would the one here https://gist.github.com/oblakeobjet/7546272 but this does not show how to introduce the data, and it is difficult when you are learning from scratch online. Without tiring very kind people, I hope the solution is easy as changing a line of code and that someone can help. My code

#!/usr/bin/python
# vim: set fileencoding=UTF8

import numpy as np
import pandas as pd
from matplotlib.mlab import griddata
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

#fig = plt.figure(figsize=(10,8))  #when uncommented draws map with colorbar but no contours

#prepare a basemap

m = Basemap(projection = 'merc',llcrnrlon = 21, llcrnrlat = -18, urcrnrlon = 34, urcrnrlat = -8, resolution='h')

# draw country outlines.

m.drawcountries(linewidth=0.5, linestyle='solid', color='k', antialiased=1, ax=None, zorder=None)
m.drawmapboundary(fill_color = 'white')
m.fillcontinents(color='coral',lake_color='blue')
parallels = np.arange(-18, -8, 2.)
m.drawparallels(np.arange(-18, -8, 2.), color = 'black', linewidth = 0.5)
m.drawparallels(parallels,labels=[True,False,False,False])
meridians = np.arange(22,34, 2.)
m.drawmeridians(np.arange(21,36, 2.), color = '0.25', linewidth = 0.5)
m.drawmeridians(meridians,labels=[False,False,False,True])

fig = plt.figure(figsize=(10,8))       # At this position or commented draws teo figures side by side

#-- Read the data.

data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True)

#-- Now gridding data.  First making a regular grid to interpolate onto

numcols, numrows = 300, 300
xi = np.linspace(data.Lon.min(), data.Lon.max(), numcols)
yi = np.linspace(data.Lat.min(), data.Lat.max(), numrows)
xi, yi = np.meshgrid(xi, yi)

#-- Interpolating at the points in xi, yi

x, y, z = data.Lon.values, data.Lat.values, data.Z.values
zi = griddata(x, y, z, xi, yi)

#-- Display and write the results

m = plt.contourf(xi, yi, zi)
plt.scatter(data.Lon, data.Lat, c=data.Z, s=100,
       vmin=zi.min(), vmax=zi.max())
fig.colorbar(m)
plt.savefig("rainfall.jpg", format="jpg")

The plots I get look like this and

and my data

32.6  -13.6   41
27.1  -16.9   43
32.7  -10.2   46
24.2  -13.6   33
28.5  -14.4   43
28.1  -12.6   33
27.9  -15.8   46
24.8  -14.8   44
31.1  -10.2   35
25.9  -13.5   24
29.1   -9.8   10
25.8  -17.8   39
33.2  -12.3   44
28.3  -15.4   46
27.6  -16.1   47
28.9  -11.1   31
31.3   -8.9   39
31.9  -13.3   45
23.1  -15.3   31
31.4  -11.9   39
27.1  -15.0   42
24.4  -11.8   15
28.6  -13.0   39
31.3  -14.3   44
23.3  -16.1   39
30.2  -13.2   38
24.3  -17.5   32
26.4  -12.2   23
23.1  -13.5   27

解决方案

You're almost there, but Basemap can be temperamental, and you have to manage the z-order of plots / map details. Also, you have to transform your lon / lat coordinates to map projection coordinates before you plot them using basemap.

Here's a complete solution, which gives the following output. I've changed some colours and linewidths in order to make the whole thing more legible, YMMV. I've also scaled the size of the scatter points by the normalised 'mean' value (data['Z']) – you can simply remove it and substitute e.g. 50 if you'd prefer a constant size (it'll look like the largest marker).

Please also specify the units of rainfall, and the duration of the measurement which generated the means, if possible:

import numpy as np
import pandas as pd
from matplotlib.mlab import griddata
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
%matplotlib inline

# set up plot
plt.clf()
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, axisbg='w', frame_on=False)

# grab data
data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True)
norm = Normalize()

# define map extent
lllon = 21
lllat = -18
urlon = 34
urlat = -8

# Set up Basemap instance
m = Basemap(
    projection = 'merc',
    llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat,
    resolution='h')

# transform lon / lat coordinates to map projection
data['projected_lon'], data['projected_lat'] = m(*(data.Lon.values, data.Lat.values))

# grid data
numcols, numrows = 1000, 1000
xi = np.linspace(data['projected_lon'].min(), data['projected_lon'].max(), numcols)
yi = np.linspace(data['projected_lat'].min(), data['projected_lat'].max(), numrows)
xi, yi = np.meshgrid(xi, yi)

# interpolate
x, y, z = data['projected_lon'].values, data['projected_lat'].values, data.Z.values
zi = griddata(x, y, z, xi, yi)

# draw map details
m.drawmapboundary(fill_color = 'white')
m.fillcontinents(color='#C0C0C0', lake_color='#7093DB')
m.drawcountries(
    linewidth=.75, linestyle='solid', color='#000073',
    antialiased=True,
    ax=ax, zorder=3)

m.drawparallels(
    np.arange(lllat, urlat, 2.),
    color = 'black', linewidth = 0.5,
    labels=[True, False, False, False])
m.drawmeridians(
    np.arange(lllon, urlon, 2.),
    color = '0.25', linewidth = 0.5,
    labels=[False, False, False, True])

# contour plot
con = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu')
# scatter plot
m.scatter(
    data['projected_lon'],
    data['projected_lat'],
    color='#545454',
    edgecolor='#ffffff',
    alpha=.75,
    s=50 * norm(data['Z']),
    cmap='RdPu',
    ax=ax,
    vmin=zi.min(), vmax=zi.max(), zorder=4)

# add colour bar and title
# add colour bar, title, and scale
cbar = plt.colorbar(conf, orientation='horizontal', fraction=.057, pad=0.05)
cbar.set_label("Mean Rainfall - mm")

m.drawmapscale(
    24., -9., 28., -13,
    100,
    units='km', fontsize=10,
    yoffset=None,
    barstyle='fancy', labelstyle='simple',
    fillcolor1='w', fillcolor2='#000000',
    fontcolor='#000000',
    zorder=5)

plt.title("Mean Rainfall")
plt.savefig("rainfall.png", format="png", dpi=300, transparent=True)
plt.show()

Using matplotlib's griddata method is convenient, but it can also be slow. As an alternative, you can use scipy's griddata methods, which are both faster, and more flexible:

from scipy.interpolate import griddata as gd

zi = gd(
    (data[['projected_lon', 'projected_lat']]),
    data.Z.values,
    (xi, yi),
    method='linear')

If you use scipy's griddata method, you'll also have to determine which of the methods (nearest, linear, cubic) give the best resulting plot.

I should add that the interpolation methods demonstrated and discussed above are the simplest possible, and aren't necessarily valid for the interpolation of rainfall data. This article gives a good overview of valid approaches and considerations for use in hydrology and hydrological modelling. The implementation of these (probably using Scipy) is left as an exercise &c.

这篇关于如何使等高线图叠加在底图上的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆