<p>You have a lot of options.</p>

<p>The easiest solution for this simple task would be to use a GIS software, e.g. the free <a href="http://www.qgis.org/en/site/" rel="noreferrer">QGIS</a>. Add delimited text layer and try raster interpolation. Download a <a href="http://www.diva-gis.org/gdata" rel="noreferrer">free coastline vector</a> and clip your raster with the coastline. A few searches at <a href="https://gis.stackexchange.com/">GIS SE</a> 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. </p>

<p>Alternatively (according to your updated question), you can use Python. This will somehow give you more control of your workflow. <a href="https://matplotlib.org/basemap/" rel="noreferrer">Basemap</a> is a useful package, see e.g. <a href="http://earthpy.org/interpolation_between_grids_with_basemap.html" rel="noreferrer">this tutorial</a> for a start. Python is also free and there is a great community at <a href="https://stackoverflow.com/questions/tagged/python">SE</a> 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.) <a href="http://www.gdal.org/" rel="noreferrer">GDAL</a> is a great library. </p>

<p>You can also use <a href="https://www.r-project.org/about.html" rel="noreferrer">R</a>, that might be a smart solution if you intend to do some more demanding statistical analysis later. There are <a href="http://www.milanor.net/blog/maps-in-r-plotting-data-points-on-a-map/" rel="noreferrer">some tutorials</a> that can put you on the right track. </p>

<p><a href="http://gmt.soest.hawaii.edu/" rel="noreferrer">GMT</a> should also able to make what you need and there is a <a href="https://github.com/GenericMappingTools/gmt-python" rel="noreferrer">python interface</a>, at least under development. </p>

<p>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. </p>

<p>Enjoy your map-making! </p>

<hr>

<p>Update: </p>

<p>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. </p>

<pre><code>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()
</code></pre>

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

<p><a href="https://i.sstatic.net/RoDgb.png" rel="noreferrer"><img src="https://i.sstatic.net/RoDgb.png" alt="enter image description here"></a></p>

CC BY-SA 3.0
2017-08-11T15:17:25.243Z
Thank you for that. However, could you be a bit comprehensive giving some commands for python please?
Stavros Keppas
CC BY-SA 3.0
2017-08-11T18:42:45.237Z
That basemap interpolation tutorial does not do what Stavros is looking for.
milancurcic
CC BY-SA 3.0
2017-08-11T23:08:01.817Z
@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.
user2821
CC BY-SA 3.0
2017-08-13T15:43:15.153Z
@StavrosKeppas I added some code, but I'd go for GIS first.
user2821
CC BY-SA 3.0
2017-08-13T16:45:49.177Z
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.
Stavros Keppas
CC BY-SA 3.0
2017-08-14T01:40:28.020Z
drawlsmask works a bit differently. For more and better options I'd have look at at the gdal library.
user2821
CC BY-SA 4.0
2022-03-22T04:48:02.433Z
Other interpolating functions will produce a better result than cubic interpolation. https://doi.org/10.3390/atmos12101318
norman_h
CC BY-SA 3.0
2017-08-11T12:35:18.487Z
2017-08-13T15:52:39.450Z
2017-08-13T15:52:39.450Z
user2821
<p>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.</p>

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

CC BY-SA 3.0
2017-08-11T10:38:46.427Z
2017-08-11T10:38:46.427Z
Communisty
<p>It is straightforward to do so with <code>numpy</code>, <code>scipy.interpolate.griddata</code>, and <code>matplotlib</code>. Here is an example:</p>

<pre><code>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 &gt; 0.5) &amp; (xi &lt; 0.6) &amp; (yi &gt; 0.5) &amp; (yi &lt; 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)
</code></pre>

<p>Result:</p>

<p><a href="https://i.sstatic.net/AMWOC.png" rel="noreferrer"><img src="https://i.sstatic.net/AMWOC.png" alt="enter image description here"></a></p>

<p>How to use this:</p>

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

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

CC BY-SA 3.0
2017-08-12T13:15:56.643Z
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,
Stavros Keppas
CC BY-SA 3.0
2017-08-12T15:21:19.840Z
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](https://www.ngdc.noaa.gov/mgg/global/) which I use often for this purpose.
milancurcic
CC BY-SA 3.0
2017-08-13T12:18:56.600Z
As I am quite new in python could you please give me a more detailed description of what to do?
Stavros Keppas
CC BY-SA 3.0
2017-08-13T13:41:25.357Z
@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.
uhoh
CC BY-SA 4.0
2019-05-23T06:14:35.810Z
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?
ery
CC BY-SA 4.0
2019-05-24T18:52:41.050Z
@ery If you're taking nearest value, I assume you're not concerned with preserving energy or momentum in your field. In that case, just apply the same procedure (as for scalar) to each velocity component.
milancurcic
CC BY-SA 4.0
2020-01-06T19:30:33.163Z
this method does not take into account that the taken data are on a curved sphere, correct?
K.-Michael Aye
CC BY-SA 4.0
2020-01-16T14:58:30.810Z
@K.-MichaelAye Interpolation is independent from the projection, in other words, it interpolates in whatever space the input data is. Here, the input data is on a sphere, and so is the result. If the result needs to be cast to, say, Cartesian space, apply the transformation after the interpolation step.
milancurcic
CC BY-SA 4.0
2020-01-17T15:03:35.430Z
But latitude and longitude are not linear on a sphere, specifically near the poles, that should influence the weighting of interpolation, shouldn't it?
K.-Michael Aye
CC BY-SA 4.0
2020-01-19T19:55:53.163Z
I see. Correct, if you need conserved interpolated quantity in physical (Cartesian) space, project it from spherical (lon, lat) to Cartesian (x, y) first, then interpolate. If the data is in in (lat, lon) and this answer taken as is, the result will only be conservative in (lat, lon) space, not (x, y).
milancurcic
CC BY-SA 4.0
2021-06-02T23:23:12.680Z
`griddata(..., method='cubic')` might provide more realistic looking results..
naught101
CC BY-SA 3.0
2017-08-11T18:37:51.740Z
2017-08-11T18:37:51.740Z
milancurcic
<p>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?</p>

CC BY-SA 3.0
2017-08-11T18:21:43.080Z
this sounds like a question for GIS stack exchange
f.thorpe
<meteorology><temperature><python><interpolation>
How to interpolate scattered data to a regular grid in Python?
Stavros Keppas
# How to interpolate scattered data to a regular grid in Python?

(Asked by: username_0 on 2017-08-11)

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?

**Tags:** meteorology, temperature, python, interpolation

> **Comments:**
> **username_1**: this sounds like a question for GIS stack exchange  

---

## 3 answers

### Answer by username_3 on 2017-08-11. Score: 28 | Accepted answer

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](https://i.sstatic.net/AMWOC.png)](https://i.sstatic.net/AMWOC.png)

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.

> **Comments:**
> **username_0**: 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,  
> **username_3**: 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](https://www.ngdc.noaa.gov/mgg/global/) which I use often for this purpose.  
> **username_0**: As I am quite new in python could you please give me a more detailed description of what to do?  
> **username_6**: @username_0 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.  
> **username_7**: 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?  
> **username_3**: @username_7 If you're taking nearest value, I assume you're not concerned with preserving energy or momentum in your field. In that case, just apply the same procedure (as for scalar) to each velocity component.  
> **username_8**: this method does not take into account that the taken data are on a curved sphere, correct?  
> **username_3**: @username_8 Interpolation is independent from the projection, in other words, it interpolates in whatever space the input data is. Here, the input data is on a sphere, and so is the result. If the result needs to be cast to, say, Cartesian space, apply the transformation after the interpolation step.  
> **username_8**: But latitude and longitude are not linear on a sphere, specifically near the poles, that should influence the weighting of interpolation, shouldn't it?  
> **username_3**: I see. Correct, if you need conserved interpolated quantity in physical (Cartesian) space, project it from spherical (lon, lat) to Cartesian (x, y) first, then interpolate. If the data is in in (lat, lon) and this answer taken as is, the result will only be conservative in (lat, lon) space, not (x, y).  
> **username_9**: `griddata(..., method='cubic')` might provide more realistic looking results..  

---

### Answer by username_2 on 2017-08-11. Score: 13

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](http://www.qgis.org/en/site/). Add delimited text layer and try raster interpolation. Download a [free coastline vector](http://www.diva-gis.org/gdata) and clip your raster with the coastline. A few searches at [GIS SE](https://gis.stackexchange.com/) 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](https://matplotlib.org/basemap/) is a useful package, see e.g. [this tutorial](http://earthpy.org/interpolation_between_grids_with_basemap.html) for a start. Python is also free and there is a great community at [SE](https://stackoverflow.com/questions/tagged/python) 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](http://www.gdal.org/) is a great library.

You can also use [R](https://www.r-project.org/about.html), that might be a smart solution if you intend to do some more demanding statistical analysis later. There are [some tutorials](http://www.milanor.net/blog/maps-in-r-plotting-data-points-on-a-map/) that can put you on the right track.

[GMT](http://gmt.soest.hawaii.edu/) should also able to make what you need and there is a [python interface](https://github.com/GenericMappingTools/gmt-python), 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](https://i.sstatic.net/RoDgb.png)](https://i.sstatic.net/RoDgb.png)

> **Comments:**
> **username_0**: Thank you for that. However, could you be a bit comprehensive giving some commands for python please?  
> **username_3**: That basemap interpolation tutorial does not do what username_0 is looking for.  
> **username_2**: @username_3 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.  
> **username_2**: @username_0 I added some code, but I'd go for GIS first.  
> **username_0**: 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.  
> **username_2**: drawlsmask works a bit differently. For more and better options I'd have look at at the gdal library.  
> **username_4**: Other interpolating functions will produce a better result than cubic interpolation. https://doi.org/10.3390/atmos12101318  

---

### Answer by username_5 on 2017-08-11. Score: 0

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.


---

