How To Plot Interpolate Station Data Over A Map (basemap)?

I've got a file with accumulated rainfall over a month in 5 stations. There are lat, lon and rain data in csv file. My file is just like that: Out[18]: lat lon rai

Solution 1:

There is no serious error in your code. (1) I suspect that (Lat, long) is reversed because all points are located in South Atlantic. (2) Content in data file should not have extra spaces on the first line.

This is the data. No extra spaces on the first line.

lat,lon,rain
0, -48.379000, -1.067000,  213.0
1, -48.435548, -1.401513,  157.2
2, -48.482217, -1.449707,  147.0
3, -48.457779, -1.249272,  182.6
4, -48.479847, -1.308735,   49.4

Here is the working code based on yours, and the resulting map. Note that I commented out the parts that plot land features. They are useless with the current data.

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
from math import ceil

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)

data = pd.read_csv('rainfall-2010-12.txt', na_values=['NaN'], sep=',')
norm = Normalize()

#mapextent
lllon = data['lon'].min()
lllat = data['lat'].min()
urlon = data['lon'].max()
urlat = data['lat'].max()

#Basemap
pad = 0.01# padding around map extents
m = Basemap(
    projection = 'merc', \
    llcrnrlon = lllon - pad, \
    llcrnrlat = lllat - pad, \
    urcrnrlon = urlon + pad, \
    urcrnrlat = urlat + pad, \
    resolution='h', \
    ax=ax)

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

#griddata
numcols, numrows = 100, 100
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.rain.values
zi = griddata(x, y, z, xi, yi, interp='linear')

# contours plot
conf = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu')  # filled contour
cont = m.contour(xi, yi, zi, zorder=5)                           # line contour# *** If the location is on land, uncomment this block of code ***# draw some map features#m.drawcoastlines()#m.drawrivers(color='b')#m.fillcontinents(color='lightyellow')# paint the ocean
watercolor=np.array([0.6, 0.8, 0.9])
m.drawlsmask(land_color=watercolor, ocean_color=watercolor, lakes=False, resolution='i')

# draw parallels and meridians; also accompanying labels
m.drawparallels(np.arange(ceil(lllat*100)/100, ceil(urlat*100)/100, 0.05), \
                labels=[1,0,0,0], linewidth=0.2, dashes=[5, 3])   # draw parallels
m.drawmeridians(np.arange(ceil(lllon*100)/100, ceil(urlon*100)/100, 0.05), \
                labels=[0,0,0,1], linewidth=0.2, dashes=[5, 3])   # draw meridians

cbar = plt.colorbar(conf, orientation='horizontal', fraction=.057, pad=0.075)
cbar.set_label("Rainfall - mm")

plt.title("Rainfall")
plt.show()

(Edit 1) With the same code above, but 2 different setups:

(1) Python 2.7.14 / Basemap 1.1.0 (conda-forge)

(2) Python 3.5.4 / Basemap 1.1.0 (conda-forge)

The resulting plots are visually the same as shown below, Left image: setup 1, right image: setup 2.