How to interpolate scattered data to a regular grid in Python?

Question

I have three txt files for longitude, latitude and temperature (or let's say three lists lon, lat, temp) from scattered weather station in the UK. I would like firstly to interpolate these data in order to get a nice colourful map of temperatures. Then, I would like to plot this interpolated temperature layer only over the land mask (thus over the british isles and not over the sea). Is that possible with Python and how?

Answer 1

It is straightforward to do so with numpy, scipy.interpolate.griddata, and matplotlib. Here is an example:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata

# data coordinates and values
x = np.random.random(100)
y = np.random.random(100)
z = np.random.random(100)

# target grid to interpolate to
xi = yi = np.arange(0,1.01,0.01)
xi,yi = np.meshgrid(xi,yi)

# set mask
mask = (xi > 0.5) & (xi < 0.6) & (yi > 0.5) & (yi < 0.6)

# interpolate
zi = griddata((x,y),z,(xi,yi),method='linear')

# mask out the field
zi[mask] = np.nan

# plot
fig = plt.figure()
ax = fig.add_subplot(111)
plt.contourf(xi,yi,zi,np.arange(0,1.01,0.01))
plt.plot(x,y,'k.')
plt.xlabel('xi',fontsize=16)
plt.ylabel('yi',fontsize=16)
plt.savefig('interpolated.png',dpi=100)
plt.close(fig)

Result:

How to use this:

• x and y are locations of points - these correspond to lon and lat values of your stations;
• z are the values of points - this corresponds to your temperature observations from stations;
• xi and yi are target grid axes - these will be your target longitude and latitude coordinates, which must match your landmask field;
• zi is the result;
• This example includes a simple way to mask the field. You should replace this mask with the landmask on your grid.

Notice also the method argument to griddata. Besides linear, this can also be cubic or nearest. I suggest you play with each to see what yields the best result for your dataset.

Answer 2

You have a lot of options.

The easiest solution for this simple task would be to use a GIS software, e.g. the free QGIS. Add delimited text layer and try raster interpolation. Download a free coastline vector and clip your raster with the coastline. A few searches at GIS SE can help you out if you get stuck. With a GIS option, it is easy to also plot e.g. cities or extract the interpolated temperature for a location.

Alternatively (according to your updated question), you can use Python. This will somehow give you more control of your workflow. Basemap is a useful package, see e.g. this tutorial for a start. Python is also free and there is a great community at SE and elsewhere. numpy and scipy are good packages for interpolation and all array processes. For more complicated spatial processes (clip a raster from a vector polygon e.g.) GDAL is a great library.

You can also use R, that might be a smart solution if you intend to do some more demanding statistical analysis later. There are some tutorials that can put you on the right track.

GMT should also able to make what you need and there is a python interface, at least under development.

Probably, you'd like to spend some effort on picking the right interpolation method and make sure that your grid is the best estimate for the actual values.

Enjoy your map-making!

---

Update:

I think that GIS would be the first approach, but as you asked for some Python commands, here is a sloppy example of how to use Python, basemap and scipy for your application. It can be greatly improved by creating a mask from a shapefile and, as mentioned, a sensitive use of interpolation method.

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

#Define mapframe
lllon = -11
lllat = 49
urlon = 2
urlat = 61

# Make some toy data, random points + corners
n = 10 # no of stations
lat = np.random.uniform(low=lllat+2, high=urlat-2, size=n)
lat = np.append(lat, [lllat, urlat, urlat, lllat]) 
lon = np.random.uniform(low=lllon+2, high=urlon-2, size=n)
lon = np.append(lon, [lllon, urlon, lllon, urlon])
temp = np.random.randn(n+4) + 8 # British summer?

# set up basemap chose projection!
m = Basemap(projection = 'merc', resolution='i', 
    llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat) 

# transform coordinates to map projection m
m_lon, m_lat = m(*(lon, lat))

# generate grid data
numcols, numrows = 240, 240
xi = np.linspace(m_lon.min(), m_lon.max(), numcols)
yi = np.linspace(m_lat.min(), m_lat.max(), numrows)
xi, yi = np.meshgrid(xi, yi)

# interpolate, there are better methods, especially if you have many datapoints
zi = griddata((m_lon,m_lat),temp,(xi,yi),method='cubic')

fig, ax = plt.subplots(figsize=(12, 12))

# draw map details
m.drawmapboundary(fill_color = 'skyblue', zorder = 1)

# Plot interpolated temperatures
m.contourf(xi, yi, zi, 500, cmap='magma', zorder = 2)

m.drawlsmask(ocean_color='skyblue', land_color=(0, 0, 0, 0), lakes=True, zorder = 3)

cbar = plt.colorbar()
plt.title('Temperature')

plt.show()

(This is modified code, used for something else. For detailed questions, other forums are more suitible. )

Answer 3

First you would have to read the data for example in matlab. Then you can get the whole field interpolated with the function griddata in matlab. And there also exist a landmask -function that allows you to further plot a map of your liking.

So yes this is possible and there are the necessary functions at least in matlab and I would guess that in other languages too.