So, meanwhile, I found a good solution thanks to a colleague who at some point implemented something similar (based on this blog post ).
Old, very slow approach using slim
Firstly, here my reference implementation uses a slender, which is just a slightly refined version of my first “naive” approach in opening a message. It is easy to understand and work, but very slow.
import numpy as np from shapely.geometry import Polygon, Point def spatial_contour_frequency_shapely(paths,lon,lat): frequency = np.zeros(lon.shape) contours = [Polygon(path) for path in paths] for (i,j),v in np.ndenumerate(frequency): pt = Point([lon[i,j],lat[i,j]]) for contour in contours: if contour.contains(pt): frequency[i,j] += 1 return frequency
New, very fast solution using PIL
My (almost) final solution is no longer used, but instead uses image manipulation methods from PIL (Python Imaging Library). This solution is much, much, much faster, although in this form only for regular rectangular grids (see. Comment below).
import numpy as np from PIL import Image, ImageDraw def _spatial_contour_frequency_pil(paths,lon,lat,regular_grid=False, method_ind=None): def get_indices(points,lon,lat,tree=None,regular=False): def get_indices_regular(points,lon,lat): lon,lat = lon.T,lat.T def _get_ij(lon,lat,x,y): lon0 = lon[0,0] lat0 = lat[0,0] lon1 = lon[-1,-1] lat1 = lat[-1,-1] nx,ny = lon.shape dx = (lon1-lon0)/nx dy = (lat1-lat0)/ny i = int((x-lon0)/dx) j = int((y-lat0)/dy) return i,j return [_get_ij(lon,lat,x,y) for x,y in points] def get_indices_irregular(points,tree,shape): dist,idx = tree.query(points,k=1) ind = np.column_stack(np.unravel_index(idx,lon.shape)) return [(i,j) for i,j in ind] if regular: return get_indices_regular(points,lon,lat) return get_indices_irregular(points,tree,lon.T.shape) tree = None if not regular_grid: lonlat = np.column_stack((lon.T.ravel(),lat.T.ravel())) tree = sp.spatial.cKDTree(lonlat) frequency = np.zeros(lon.shape) for path in paths: path_ij = get_indices(path,lon,lat,tree=tree,regular=regular_grid) raster_poly = Image.new("L",lon.shape,0) rasterize = ImageDraw.Draw(raster_poly) rasterize.polygon(path_ij,1) mask = sp.fromstring(raster_poly.tobytes(),'b') mask.shape = raster_poly.im.size[1],raster_poly.im.size[0] frequency += mask return frequency
It should be noted that the result of these two approaches is not the same. The functions identified by the PIL approach are slightly larger than those defined with the effective approach, but in reality one is not better than the other.
Delays
Below are some timings created with a reduced data set (but not semi-artificial example data from the initial message):
spatial_contour_frequency/shapely : 191.8843 spatial_contour_frequency/pil : 0.3287 spatial_contour_frequency/pil-tree-inside : 2.3629 spatial_contour_frequency/pil-regular_grid : 0.3276
The most time-consuming step is to search for indices on an irregular grid of lon / lat contour points. The biggest part of the time is the cKDTree construct, so I moved it from get_indices . To put this in perspective, pil-tree-inside is the version in which the tree is created inside get_indices . pil-regular-grid has regular_grid=True , which for my dataset gives incorrect results, but gives an idea of how quickly this will work on a regular grid.
In general, I managed to largely eliminate the effect of an irregular grid ( pil vs. pil-regular-grid ), which I could hope for at the beginning! :)