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?
It is straightforward to do so with
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)
How to use this:
yare locations of points - these correspond to
latvalues of your stations;
zare the values of points - this corresponds to your temperature observations from stations;
yiare target grid axes - these will be your target longitude and latitude coordinates, which must match your landmask field;
ziis 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
linear, this can also be
nearest. I suggest you play with each to see what yields the best result for your dataset.
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.
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!
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. )
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.