1
$\begingroup$

I have another question, is there a package that interpolates precipitation data taking into account mountains and oceans? I have so far used Numpy and Basemap but as you can see in the code, the data from Europe affect the data from North America! This distorts the results and must be prevented at all costs.

(I am aware that my interpolation is just very simplified)

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patches import Path, PathPatch

# data coordinates and values
x = np.array([6.47, 4.8, 1.94, -4.57, 5.78, 3.95, 5.45, 0.38, 7.13, 12.43, 11.59, -1.083333333, -60.996, -83.144 ])
y = np.array([46.37, 43.9, 47.83, 48.36, 44.95, 43.57, 43.45, 43.12, 43.92, 51.35, 48.22, 52.88333333, 50.958, 37.02])
z = np.array([-29.5, -27.6, -32.5, -32.4, -16.0, -0.1, -12.3, -15.2, -36.2,  -20.6,-38.2, -38.2, -2.0, -4.0  ])

# target grid to interpolate to
xi =np.arange(-100, 100, 1, dtype='float64')
yi =np.arange(-100,100, 1, dtype='float64')
xi, yi = np.meshgrid(xi, yi)

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

fig, ax = plt.subplots(figsize=(10,10))
m = Basemap(llcrnrlon=x.min()-0.1,llcrnrlat=y.min()-0.1,urcrnrlon=x.max()+0.1,urcrnrlat=y.max()+0.1, projection='merc', resolution='h',area_thresh=1000.,ax=ax)

m.drawcoastlines() #draw coastlines on the map
x,y=m(xi, yi) # convert the coordinates into the map scales
cs=ax.contourf(x, y, zi, np.linspace(-40, 0),extend='both',cmap='jet') #plot the data on the map.
plt.scatter(x, y)
plt.show()

Oh and if someone knows how I can delete these meshgrid points from Numpy, I would be very happy about a short answer, because this I have somehow not found (which I myself find very strange)

Edit: Normally the ocean is covered with a mask, but this would make the code much too long

Result:enter image description here

$\endgroup$
4
  • 2
    $\begingroup$ Have you considered putting constraints on the data? Such as creating a zone for Europe & another for North America & specifying that only data from within each zone can only be used for each land mass. $\endgroup$
    – Fred
    May 14 at 0:17
  • 1
    $\begingroup$ Hi @Fred, thanks for your reply that sounds promising, do you have a website or link explaining this step? All the best $\endgroup$
    – Weiss
    May 14 at 7:26
  • $\begingroup$ Depending on the application & what software is available, some packages allow regions to be specified using a screen digitized polygon & then specifying whether everything inside or outside of the polygon is to be included or excluded. A more laborious way would be to write a script to split the data file. One criteria that could be used in this situation is the co-ordinate for longitude. Choose something west of Iceland & everything west of that reports to the US date file & everything east of that reports to the Europe data file. $\endgroup$
    – Fred
    May 14 at 8:12
  • $\begingroup$ May I suggest checking out verde? It's much more suitable than griddata for the sort of thing you're doing. And I suggest replacing basemap, which is deprecated, with cartopy. $\endgroup$
    – kwinkunks
    May 14 at 19:31

0

Your Answer

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

Browse other questions tagged or ask your own question.