10
$\begingroup$

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?

$\endgroup$
  • 4
    $\begingroup$ this sounds like a question for GIS stack exchange $\endgroup$ – farrenthorpe Aug 11 '17 at 18:21
15
$\begingroup$

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:

enter image description here

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.

$\endgroup$
  • $\begingroup$ Thank you for that. Is there an easy way to use matplotlib to plot only over the land? Apparently I cannot use (mask = (xi > 0.5) & (xi < 0.6) & (yi > 0.5) & (yi < 0.6)) to define the landmask of the UK, $\endgroup$ – Stavros Keppas Aug 12 '17 at 13:15
  • $\begingroup$ You need mask data, e.g. a 2-d field with zeros over water and ones over land. If you don't have one, you can make one based on topography data, for example ETOPO 01 which I use often for this purpose. $\endgroup$ – milancurcic Aug 12 '17 at 15:21
  • $\begingroup$ As I am quite new in python could you please give me a more detailed description of what to do? $\endgroup$ – Stavros Keppas Aug 13 '17 at 12:18
  • 4
    $\begingroup$ @StavrosKeppas at some point you have to try to do it yourself. If you post here, or in GIS stackexchange, or in stackoverflow, you should try to do it as best you can. "Please write the program I need" questions often go unanswered. The best way to learn to use python is to try for a while. If it doesn't work, then you post your code that does not work and ask for assistance. This shows effort which is a big plus, and is likely encourage more helpful answers. Also, you can ask as many (good) questions as you like, don't keep extending the question in comments. $\endgroup$ – uhoh Aug 13 '17 at 13:41
  • $\begingroup$ I am trying to adapt this to resmple CFD result into a coarser mesh. This works excellent for single value of z (like temperature). What if the z is a velocity, which will have three values for the vector components. I am not looking into interpolation, just want to take the nearest value. Thoughts? $\endgroup$ – ery May 23 at 6:14
11
$\begingroup$

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. )

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for that. However, could you be a bit comprehensive giving some commands for python please? $\endgroup$ – Stavros Keppas Aug 11 '17 at 15:17
  • $\begingroup$ That basemap interpolation tutorial does not do what Stavros is looking for. $\endgroup$ – milancurcic Aug 11 '17 at 18:42
  • $\begingroup$ @milancurcic no, but it can hopefully put him on the right track. numpy and scipy are good packages for interpolation and all array processes, as shown in your example. The tricky thing is often to get basemap to do what one intend. $\endgroup$ – Tactopoda Aug 11 '17 at 23:08
  • $\begingroup$ @StavrosKeppas I added some code, but I'd go for GIS first. $\endgroup$ – Tactopoda Aug 13 '17 at 15:43
  • $\begingroup$ Thank you.That's great. I would like to plot the coastline as well. I added this two lines: m.drawcoastlines() m.drawcountries() but I get a coastline which is not a continuous line. $\endgroup$ – Stavros Keppas Aug 13 '17 at 16:45
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.