title: Nearest Neighbor or KD Tree with 2D Lat and Lon Python
body: <h2>What I am trying to do:</h2>
  
  <p>What I want to do is have a function for nearest neighbor where the user provides the latitude and longitude arrays and a singular point. Within the function I do not mind applying a function from a different library (such as the scipy kdtree). Previously I applied the below function that works for finding the nearest neighbor for 1D lat and lon arrays: </p>
  
  <pre><code>    import numpy as np
      def find_index_of_nearest_xy(lat_array, lon_array, lat_point, lon_point, num_neighbors = 4):
          dist = np.sqrt(np.abs((lon_array-lon_point)**2 + (lat_array - lat_point)**2))
          sorted_dist = sorted(dist)
          values = np.array(sorted_dist[:num_neighbors])
          indices = np.in1d(dist,values)
  
          return indices
  </code></pre>
  
  <p>I know there is a more efficient way to go about finding the index of the nearest point but this method was quick for the application I was doing previously. </p>
  
  <h2>More details:</h2>
  
  <p>The issue I am having now is that my function only works if lat and lon are 1D. With my new dataset, my lat and lon arrays are both 2D each having a shape of (5500, 5500) so raveling it and applying my function would be very inefficient as I am performing this analysis every 5 minutes in a 3 month period. An example of this approach would be: </p>
  
  <pre><code>    #current_AERONET is a pandas dataframe with aeronet lat, lon, time, and AOT data
      aeronet_lat = current_AERONET['Lat'].tolist()
      aeronet_lon = current_AERONET['Lon'].tolist()
      for ahi in ahi_files:
          openfile = netcdf.netcdf_file(ahi, 'r')
          ahi_lat = np.copy(np.squeeze(openfile.variables['lat'][:])) #(5500, 5500)
          ahi_lon = np.copy(np.squeeze(openfile.variables['lon'][:])) #(5500,5500)
          new_lat = np.ravel(lat) #(30250000,)
          new_lon = np.ravel(lon) #(30250000,)            
  
          idx = find_index_of_nearest_xy(new_lat, new_lon, aeronet_lat, aeronet_lon)
          current_lat = new_lat[idx]
          current_lon = new_lon[idx]
  </code></pre>
  
  <p>I found <a href="https://stackoverflow.com/questions/24850133/python-kd-tree-nearest-neigbour-where-distance-is-greater-than-zero">this link</a> that explains how to perform a kd tree using lat and lon coordinates but again the lat and lon are 1D. </p>
  
  <p>I want to return the indices since I am applying it to my lat, lon, and data arrays at a later step. </p>
  
  <p>So my general question is can anyone explain/provide an example (or a link) of how to extract a specific number of neighbors (in this case 4) for a specified latitude and longitude coordinate (for example 36.538537 N, 126.330194 E) provided that latitude and longitude arrays are 2D (i.e. the dataset I need to extract neighbors from).  </p>
  
  <p>If it is useful, I am working with satellite visible data. </p>
  
  <h2>What I have attempted already:</h2>
  
  <p>Previously I have attempted to ravel the array (as shown above) and that is how I discovered that it is extremely inefficient. I have also had it return True/False values then later looked for where my array was True but that was even more inefficient than the first approach. I attempted using nsmallest from the heapq library with no luck since I did not understand how to use the library appropriately. </p>
  
  <p>I have performed an extensive search to attempt to find a solution for my problem with no prevail, so I am out options and am now asking the community.</p>