<jupyter_start><jupyter_text>Using numpy and KD-trees with netCDF dataThere is now a Unidata Developer's [Blog entry](http://www.unidata.ucar.edu/blogs/developer/en/entry/accessing_netcdf_data_by_coordinates) accompanying this iPython notebook.The goal is to demonstrate how to quickly access netCDF data based on geospatial coordinates instead of array indices.- First we show a naive and slow way to do this, in which we also have to worry about longitude anomalies- Then we speed up access with numpy arrays- Next, we demonstrate how to eliminate longitude anomalies- Finally, we use a kd-tree data structure to significantly speed up access by coordinates for large problems  Getting data by coordinates from a netCDF FileLet's look at a netCDF file from the *Atlantic Real-Time Ocean Forecast System*.  If you have cloned the [Unidata 2013 TDS-Python Workshop](https://github.com/Unidata/tds-python-workshop), this data is already available in 'data/rtofs_glo_3dz_f006_6hrly_reg3.nc'.  Otherwise you can get it from [rtofs_glo_3dz_f006_6hrly_reg3.nc](https://github.com/Unidata/tds-python-workshop/blob/master/data/rtofs_glo_3dz_f006_6hrly_reg3.nc).  Looking at netCDF metadata from PythonIn iPython, we could invoke the **ncdump** utility like this:        filename = 'data/rtofs_glo_3dz_f006_6hrly_reg3.nc'    !ncdump -h $filename    *if* we know that a recent version of **ncdump** is installed thatcan read compressed data from netCDF-4 classic model files.Alternatively, we'll use the netCDF4python package to show information aboutthe file in a form that's somewhat less familiar, but contains the informationwe need for the subsequent examples.  This works for any netCDF file format:<jupyter_code>import netCDF4
filename = 'data/rtofs_glo_3dz_f006_6hrly_reg3.nc'
ncfile = netCDF4.Dataset(filename, 'r')
print ncfile                # shows global attributes, dimensions, and variables
ncvars = ncfile.variables   # a dictionary of variables
# print information about specific variables, including type, shape, and attributes
for varname in ['temperature', 'salinity', 'Latitude', 'Longitude']:
    print ncvars[varname]<jupyter_output><type 'netCDF4.Dataset'>
root group (NETCDF4_CLASSIC file format):
    Conventions: CF-1.0
    title: HYCOM ATLb2.00
    institution: National Centers for Environmental Prediction
    source: HYCOM archive file
    experiment: 90.9
    history: archv2ncdf3z
    dimensions = (u'MT', u'Y', u'X', u'Depth')
    variables = (u'MT', u'Date', u'Depth', u'Y', u'X', u'Latitude', u'Longitude', u'u', u'v', u'temperature', u'salinity')
    groups = ()

<type 'netCDF4.Variable'>
float32 temperature(u'MT', u'Depth', u'Y', u'X')
    coordinates: Longitude Latitude Date
    standard_name: sea_water_potential_temperature
    units: degC
    _FillValue: 1.26765e+30
    valid_range: [ -5.07860279  11.14989948]
    long_name:   temp [90.9H]
unlimited dimensions = (u'MT',)
current size = (1, 10, 850, 712)

<type 'netCDF4.Variable'>
float32 salinity(u'MT', u'Depth', u'Y', u'X')
    coordinates: Longitude Latitude Date
    standard_name: sea_water_salinity
    units: psu
    _FillValue: 1.26765e+30
    valid[...]<jupyter_text>Here's a sparse picture (every 25th point on each axis) of what the grid looks like on which Latitude, Longitude, Temperature, Salinity, and other variables are defined:![Example lat-lon grid](http://www.unidata.ucar.edu/software/netcdf/workshops/images/20130811_rew_blog_grid.png)  Example query: sea surface temperature and salinity at 50N, 140W?- So **Longitude** and **Latitude** are 2D netCDF variables of shape 850 x 712, indexed by **Y** and **X** dimensions- That's 605200 values for each- There's no _direct_ way in this file (and many netCDF files) to compute grid indexes from coordinates via a coordinate system and projection parameters.  Instead, we have to rely on the latitude and longitude auxiliary coordinate variables, as required by the CF conventions for data not on a simple lat,lon grid.- To get the temperature at 50N, 140W, we need to find **Y** and **X** indexes **iy** and **ix** such that (**Longitude[iy, ix]**, **Latitude[iy, ix]**) is "close" to (50.0, -140.0).  Naive, slow way using nested loops- Initially, for simplicity, we just use Euclidean distance squared, as if the Earth is flat, latitude and longitude are $x$- and $y$-coordinates, and the distance squared between points $(lat_1,lon_1)$ and $(lat_0,lon_0)$ is $( lat_1 - lat_0 )^2 + ( lon_1 - lon_0 )^2$.- Note: these assumptions are wrong near the poles and on opposite sides of longitude boundary discontinuity.- So, keeping things simple, we want to find **iy** and **ix** to minimize    ``(Latitude[iy, ix] - lat0)**2 + (Longitude[iy, ix] - lon0)**2``    ![Flat Earth](http://www.unidata.ucar.edu/software/netcdf/workshops/images/unproj.gif)  Reading netCDF data into numpy arraysTo access netCDF data, rather than just metadata, we will also need NumPy:- A Python library for scientific programming.- Supports n-dimensional array-based calculations similar to Fortran and IDL.- Includes fast mathematical functions to act on scalars and arrays.With the Python netCDF4 package, using "[ ... ]" to index a netCDF variable object reads or writes a numpy array from the associated netCDF file.The code below reads latitude and longitude values into 2D numpy arrays named **latvals** and **lonvals**:  First version: slow and spatially challengedHere's a function that uses simple nested loops to find indices that minimize the distance to the desired coordinates, written as if using Fortran or C rather than Python.  We'll call this function in the cell following this definition ...<jupyter_code>import numpy as np
import netCDF4

def naive_slow(latvar,lonvar,lat0,lon0):
    '''
    Find "closest" point in a set of (lat,lon) points to specified point
    latvar - 2D latitude variable from an open netCDF dataset
    lonvar - 2D longitude variable from an open netCDF dataset
    lat0,lon0 - query point
    Returns iy,ix such that 
     (lonval[iy,ix] - lon0)**2 + (latval[iy,ix] - lat0)**2
    is minimum.  This "closeness" measure works badly near poles and
    longitude boundaries.
    '''
    # Read from file into numpy arrays
    latvals = latvar[:]
    lonvals = lonvar[:]
    ny,nx = latvals.shape
    dist_sq_min = 1.0e30
    for iy in range(ny):
        for ix in range(nx):
            latval = latvals[iy, ix]
            lonval = lonvals[iy, ix]
            dist_sq = (latval - lat0)**2 + (lonval - lon0)**2
            if dist_sq < dist_sq_min:
                iy_min, ix_min, dist_sq_min = iy, ix, dist_sq
    return iy_min,ix_min<jupyter_output><empty_output><jupyter_text>When we call the function above it takes several seconds to run, because it calculates distances one point at a time, for each of the 605200 $(lat, lon)$ points.  Note that once indices for the point nearest to (50, -140) are found, they can be used to access temperature, salinity, and other netCDF variables that use the same dimensions.<jupyter_code>ncfile = netCDF4.Dataset(filename, 'r')
latvar = ncfile.variables['Latitude']
lonvar = ncfile.variables['Longitude']
iy,ix = naive_slow(latvar, lonvar, 50.0, -140.0)
print 'Closest lat lon:', latvar[iy,ix], lonvar[iy,ix]
tempvar = ncfile.variables['temperature']
salvar = ncfile.variables['salinity']
print 'temperature:', tempvar[0, 0, iy, ix], tempvar.units
print 'salinity:', salvar[0, 0, iy, ix], salvar.units
ncfile.close()<jupyter_output>Closest lat lon: 49.9867 -139.982
temperature: 6.46312 degC
salinity: 32.6572 psu<jupyter_text>NumPy arrays instead of loops: fast, but still assumes flat earthThe above function is slow, because it doesn't make good use of NumPy arrays.  It's much faster to use whole array operations to eliminate loops and element-at-a-time computation.  NumPy functions that help eliminate loops include:- The `argmin()` method that returns a 1D index of the minimum value of a NumPy array- The `unravel_index()` function that converts a 1D index back into a multidimensional index<jupyter_code>import numpy as np
import netCDF4

def naive_fast(latvar,lonvar,lat0,lon0):
    # Read latitude and longitude from file into numpy arrays
    latvals = latvar[:]
    lonvals = lonvar[:]
    ny,nx = latvals.shape
    dist_sq = (latvals-lat0)**2 + (lonvals-lon0)**2
    minindex_flattened = dist_sq.argmin()  # 1D index of min element
    iy_min,ix_min = np.unravel_index(minindex_flattened, latvals.shape)
    return iy_min,ix_min

ncfile = netCDF4.Dataset(filename, 'r')
latvar = ncfile.variables['Latitude']
lonvar = ncfile.variables['Longitude']
iy,ix = naive_fast(latvar, lonvar, 50.0, -140.0)
print 'Closest lat lon:', latvar[iy,ix], lonvar[iy,ix]
ncfile.close()<jupyter_output>Closest lat lon: 49.9867 -139.982<jupyter_text>Spherical Earth with tunnel distance: fast _and_ correctThough assuming a flat Earth may work OK for this example, we'd like to not worry about whether longitudes are from 0 to 360 or -180 to 180, or whether points are close to the poles.The code below fixes this by using the square of "tunnel distance" between (lat,lon) points. This version is both fast and correct (for a _spherical_ Earth).<jupyter_code>import numpy as np
import netCDF4
from math import pi
from numpy import cos, sin

def tunnel_fast(latvar,lonvar,lat0,lon0):
    '''
    Find closest point in a set of (lat,lon) points to specified point
    latvar - 2D latitude variable from an open netCDF dataset
    lonvar - 2D longitude variable from an open netCDF dataset
    lat0,lon0 - query point
    Returns iy,ix such that the square of the tunnel distance
    between (latval[it,ix],lonval[iy,ix]) and (lat0,lon0)
    is minimum.
    '''
    rad_factor = pi/180.0 # for trignometry, need angles in radians
    # Read latitude and longitude from file into numpy arrays
    latvals = latvar[:] * rad_factor
    lonvals = lonvar[:] * rad_factor
    ny,nx = latvals.shape
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    # Compute numpy arrays for all values, no loops
    clat,clon = cos(latvals),cos(lonvals)
    slat,slon = sin(latvals),sin(lonvals)
    delX = cos(lat0_rad)*cos(lon0_rad) - clat*clon
    delY = cos(lat0_rad)*sin(lon0_rad) - clat*slon
    delZ = sin(lat0_rad) - slat;
    dist_sq = delX**2 + delY**2 + delZ**2
    minindex_1d = dist_sq.argmin()  # 1D index of minimum element
    iy_min,ix_min = np.unravel_index(minindex_1d, latvals.shape)
    return iy_min,ix_min

ncfile = netCDF4.Dataset(filename, 'r')
latvar = ncfile.variables['Latitude']
lonvar = ncfile.variables['Longitude']
iy,ix = tunnel_fast(latvar, lonvar, 50.0, -140.0)
print 'Closest lat lon:', latvar[iy,ix], lonvar[iy,ix]
ncfile.close()<jupyter_output>Closest lat lon: 49.9867 -139.982<jupyter_text>KD-Trees: faster data structure for lots of queriesWe can still do better, by using a data structure designed to support efficient nearest-neighbor queries: the [KD-tree](http://en.wikipedia.org/wiki/K-d_tree).  It works like a multidimensional binary tree, so finding the point nearest to a query point is _much_ faster than computing all the distances to find the minimum.  It takes some setup time to load all the points into the data structure, but that only has to be done once for a given set of points.      For a single point query, it's still more than twice as fast as the naive slow version above, but building the KD-tree for 605,200 points takes more time than the fast numpy search through all the points, so in this case using the KD-tree for a _single_ point query is sort of pointless ...<jupyter_code>import numpy as np
import netCDF4
from math import pi
from numpy import cos, sin
from scipy.spatial import cKDTree

def kdtree_fast(latvar,lonvar,lat0,lon0):
    rad_factor = pi/180.0 # for trignometry, need angles in radians
    # Read latitude and longitude from file into numpy arrays
    latvals = latvar[:] * rad_factor
    lonvals = lonvar[:] * rad_factor
    ny,nx = latvals.shape
    clat,clon = cos(latvals),cos(lonvals)
    slat,slon = sin(latvals),sin(lonvals)
    # Build kd-tree from big arrays of 3D coordinates
    triples = zip(np.ravel(clat*clon), np.ravel(clat*slon), np.ravel(slat))
    kdt = cKDTree(triples)
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    clat0,clon0 = cos(lat0_rad),cos(lon0_rad)
    slat0,slon0 = sin(lat0_rad),sin(lon0_rad)
    dist_sq_min, minindex_1d = kdt.query([clat0*clon0, clat0*slon0, slat0])
    iy_min, ix_min = np.unravel_index(minindex_1d, latvals.shape)
    return iy_min,ix_min
                
ncfile = netCDF4.Dataset(filename, 'r')
latvar = ncfile.variables['Latitude']
lonvar = ncfile.variables['Longitude']
iy,ix = kdtree_fast(latvar, lonvar, 50.0, -140.0)
print 'Closest lat lon:', latvar[iy,ix], lonvar[iy,ix]
ncfile.close()<jupyter_output>Closest lat lon: 49.9867 -139.982<jupyter_text>Timing the functionsIf you're curious about actual times for the versions above, the iPython notebook "%%timeit" statement gets accurate timings of all of them.  Below, we time just a single query point, in this case (50.0, -140.0).  To get accurate timings, the "%%timeit" statement lets us do untimed setup first on the same line, before running the function call in a loop.<jupyter_code>%%timeit ncfile = netCDF4.Dataset(filename,'r');latvar = ncfile.variables['Latitude'];lonvar = ncfile.variables['Longitude']
naive_slow(latvar, lonvar, 50.0, -140.0)
%%timeit ncfile = netCDF4.Dataset(filename,'r');latvar = ncfile.variables['Latitude'];lonvar = ncfile.variables['Longitude']
naive_fast(latvar, lonvar, 50.0, -140.0)
%%timeit ncfile = netCDF4.Dataset(filename,'r');latvar = ncfile.variables['Latitude'];lonvar = ncfile.variables['Longitude']
tunnel_fast(latvar, lonvar, 50.0, -140.0)
%%timeit ncfile = netCDF4.Dataset(filename,'r');latvar = ncfile.variables['Latitude'];lonvar = ncfile.variables['Longitude']
kdtree_fast(latvar, lonvar, 50.0, -140.0)
ncfile.close()<jupyter_output><empty_output><jupyter_text>Separating setup from queryThe above use of the KD-tree data structure is not the way it's meant to be used.  Instead, it should be initialized _once_ with all the k-dimensional data for which nearest-neighbors are desired, then used repeatedly on each query, amortizing the work done to build the data structure over all the following queries.  By separately timing the setup and the time required per query, the threshold for number of queries beyond which the KD-tree is faster can be determined.That's exactly what we'll do now.  We split each algorithm into two functions, a setup function and a query function.  The times per query go from seconds (the naive version) to milliseconds (the array-oriented numpy version) to microseconds (the turbo-charged KD-tree, once it's built).                                                                                                                                                                                                                                                   Rather than just using functions, we define a Class for each algorithm, do the setup in the class constructor, and provide a query method.<jupyter_code># Split naive_slow into initialization and query, so we can time them separately
import numpy as np
import netCDF4

class Naive_slow(object):
    def __init__(self, ncfile, latvarname, lonvarname):
        self.ncfile = ncfile
        self.latvar = self.ncfile.variables[latvarname]
        self.lonvar = self.ncfile.variables[lonvarname]
        # Read latitude and longitude from file into numpy arrays
        self.latvals = self.latvar[:]
        self.lonvals = self.lonvar[:]
        self.shape = self.latvals.shape

    def query(self,lat0,lon0):
        ny,nx = self.shape
        dist_sq_min = 1.0e30
        for iy in range(ny):
            for ix in range(nx):
                latval = self.latvals[iy, ix]
                lonval = self.lonvals[iy, ix]
                dist_sq = (latval - lat0)**2 + (lonval - lon0)**2
                if dist_sq < dist_sq_min:
                    iy_min, ix_min, dist_sq_min = iy, ix, dist_sq
        return iy_min,ix_min

ncfile = netCDF4.Dataset(filename, 'r')
ns = Naive_slow(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)
print 'Closest lat lon:', ns.latvar[iy,ix], ns.lonvar[iy,ix]
ncfile.close()
# Split naive_fast into initialization and query, so we can time them separately
import numpy as np
import netCDF4

class Naive_fast(object):
    def __init__(self, ncfile, latvarname, lonvarname):
        self.ncfile = ncfile
        self.latvar = self.ncfile.variables[latvarname]
        self.lonvar = self.ncfile.variables[lonvarname]        
        # Read latitude and longitude from file into numpy arrays
        self.latvals = self.latvar[:]
        self.lonvals = self.lonvar[:]
        self.shape = self.latvals.shape
        

    def query(self,lat0,lon0):
        dist_sq = (self.latvals-lat0)**2 + (self.lonvals-lon0)**2
        minindex_flattened = dist_sq.argmin()                             # 1D index
        iy_min, ix_min = np.unravel_index(minindex_flattened, self.shape) # 2D indexes
        return iy_min,ix_min

ncfile = netCDF4.Dataset(filename, 'r')
ns = Naive_fast(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)
print 'Closest lat lon:', ns.latvar[iy,ix], ns.lonvar[iy,ix]
ncfile.close()
# Split tunnel_fast into initialization and query, so we can time them separately
import numpy as np
import netCDF4
from math import pi
from numpy import cos, sin

class Tunnel_fast(object):
    def __init__(self, ncfile, latvarname, lonvarname):
        self.ncfile = ncfile
        self.latvar = self.ncfile.variables[latvarname]
        self.lonvar = self.ncfile.variables[lonvarname]        
        # Read latitude and longitude from file into numpy arrays
        rad_factor = pi/180.0 # for trignometry, need angles in radians
        self.latvals = self.latvar[:] * rad_factor
        self.lonvals = self.lonvar[:] * rad_factor
        self.shape = self.latvals.shape
        clat,clon,slon = cos(self.latvals),cos(self.lonvals),sin(self.lonvals)
        self.clat_clon = clat*clon
        self.clat_slon = clat*slon
        self.slat = sin(self.latvals)
 
    def query(self,lat0,lon0):
        # for trignometry, need angles in radians
        rad_factor = pi/180.0 
        lat0_rad = lat0 * rad_factor
        lon0_rad = lon0 * rad_factor
        delX = cos(lat0_rad)*cos(lon0_rad) - self.clat_clon
        delY = cos(lat0_rad)*sin(lon0_rad) - self.clat_slon
        delZ = sin(lat0_rad) - self.slat;
        dist_sq = delX**2 + delY**2 + delZ**2
        minindex_1d = dist_sq.argmin()                              # 1D index 
        iy_min, ix_min = np.unravel_index(minindex_1d, self.shape)  # 2D indexes
        return iy_min,ix_min

ncfile = netCDF4.Dataset(filename, 'r')
ns = Tunnel_fast(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)
print 'Closest lat lon:', ns.latvar[iy,ix], ns.lonvar[iy,ix]
ncfile.close()
# Split kdtree_fast into initialization and query, so we can time them separately
import numpy as np
import netCDF4
from math import pi
from numpy import cos, sin
from scipy.spatial import cKDTree

class Kdtree_fast(object):
    def __init__(self, ncfile, latvarname, lonvarname):
        self.ncfile = ncfile
        self.latvar = self.ncfile.variables[latvarname]
        self.lonvar = self.ncfile.variables[lonvarname]        
        # Read latitude and longitude from file into numpy arrays
        rad_factor = pi/180.0 # for trignometry, need angles in radians
        self.latvals = self.latvar[:] * rad_factor
        self.lonvals = self.lonvar[:] * rad_factor
        self.shape = self.latvals.shape
        clat,clon = cos(self.latvals),cos(self.lonvals)
        slat,slon = sin(self.latvals),sin(self.lonvals)
        clat_clon = clat*clon
        clat_slon = clat*slon
        triples = zip(np.ravel(clat*clon), np.ravel(clat*slon), np.ravel(slat))
        self.kdt = cKDTree(triples)

    def query(self,lat0,lon0):
        rad_factor = pi/180.0 
        lat0_rad = lat0 * rad_factor
        lon0_rad = lon0 * rad_factor
        clat0,clon0 = cos(lat0_rad),cos(lon0_rad)
        slat0,slon0 = sin(lat0_rad),sin(lon0_rad)
        dist_sq_min, minindex_1d = self.kdt.query([clat0*clon0,clat0*slon0,slat0])
        iy_min, ix_min = np.unravel_index(minindex_1d, self.shape)
        return iy_min,ix_min

ncfile = netCDF4.Dataset(filename, 'r')
ns = Kdtree_fast(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)
print 'Closest lat lon:', ns.latvar[iy,ix], ns.lonvar[iy,ix]
ncfile.close()<jupyter_output>Closest lat lon: 49.9867 -139.982<jupyter_text>Setup times for the four algorithms<jupyter_code>%%timeit ncfile = netCDF4.Dataset(filename, 'r')
ns = Naive_slow(ncfile,'Latitude','Longitude')
%%timeit ncfile = netCDF4.Dataset(filename, 'r')
ns = Naive_fast(ncfile,'Latitude','Longitude')
%%timeit ncfile = netCDF4.Dataset(filename, 'r')
ns = Tunnel_fast(ncfile,'Latitude','Longitude')
%%timeit ncfile = netCDF4.Dataset(filename, 'r')
ns = Kdtree_fast(ncfile,'Latitude','Longitude')<jupyter_output>1 loops, best of 3: 2.52 s per loop<jupyter_text>Query times for the four algorithms<jupyter_code>%%timeit ncfile = netCDF4.Dataset(filename, 'r'); ns = Naive_slow(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)
%%timeit ncfile = netCDF4.Dataset(filename, 'r'); ns = Naive_fast(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)
%%timeit ncfile = netCDF4.Dataset(filename, 'r'); ns = Tunnel_fast(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)
%%timeit ncfile = netCDF4.Dataset(filename, 'r'); ns = Kdtree_fast(ncfile,'Latitude','Longitude')
iy,ix = ns.query(50.0, -140.0)<jupyter_output>10000 loops, best of 3: 73.8 µs per loop<jupyter_text>In the next cell, we copy the results of the %%timeit runs into Python variables.  _(Is there a way to capture %%timeit output, so we don't have to do this manually?)_<jupyter_code>ns0,nf0,tf0,kd0 = 3.76, 3.8, 27.4, 2520  # setup times in msec
nsq,nfq,tfq,kdq = 7790, 2.46, 5.14, .0738 # query times in msec<jupyter_output><empty_output><jupyter_text>Summary of timingsThe naive_slow method is always slower than all other methods.  The naive_fast method would only be worth considering if non-flatness of the Earth is irrelevant, for example in a relatively small region not close to the poles and not crossing a longitude discontinuity.Total time for running initialization followed by N queries is:  - naive_slow:    $ns0 + nsq * N$  - naive_fast:    $nf0 + nfq * N$  - tunnel_fast:   $nt0 + ntq * N$  - kdtree_fast:   $kd0 + kdq * N$<jupyter_code>N = 10000
print N, "queries using naive_slow:", round((ns0 + nsq*N)/1000,1), "seconds"
print N, "queries using naive_fast:", round((nf0 + nfq*N)/1000,1), "seconds"
print N, "queries using tunnel_fast:", round((tf0 + tfq*N)/1000,1), "seconds"
print N, "queries using kdtree_fast:", round((kd0 + kdq*N)/1000,1), "seconds"
print
print "kd_tree_fast outperforms naive_fast above:", int((kd0-nf0)/(nfq-kdq)), "queries"
print "kd_tree_fast outperforms tunnel_fast above:", int((kd0-tf0)/(tfq-kdq)), "queries"<jupyter_output>10000 queries using naive_slow: 77900.0 seconds
10000 queries using naive_fast: 24.6 seconds
10000 queries using tunnel_fast: 51.4 seconds
10000 queries using kdtree_fast: 3.3 seconds

kd_tree_fast outperforms naive_fast above: 1054 queries
kd_tree_fast outperforms tunnel_fast above: 492 queries
{'fraction_non_alphanumeric': 0.10218105766357932, 'fraction_numerical': 0.029194587903879805, 'mean_word_length': 3.904751936361733, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 43, 'https://': 2, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '302971', 'license': 'MIT', 'n_tokens_mistral': 8191, 'n_tokens_neox': 7663, 'n_words': 2448, 'path': 'tutorials_and_examples/netcdf-by-coordinates.ipynb', 'repo_name': 'shaunwbell/ipythonnb', 'size': '23429'}