code: def property_at_lat_lon_depth_points(lat, lon, depth, quantity_ID="DENSITY"):
      
  
      nlayers = len(l1_layer_decode)
      shape = np.array(lon).shape
  
      lon1   = np.array(lon).reshape(-1)
      lat1   = np.array(lat).reshape(-1)
      depth1 = np.array(depth).reshape(-1)
      point_properties = np.empty_like(depth1)
  
      layer_depths     = np.empty((nlayers, lat1.shape[0]))
      layer_properties = np.ones((nlayers+1, lat1.shape[0])) * -99999.0  # if the point is not found it will end up in the overshoot !
  
      # should assert here that the three arrays are equal size
  
      for i in range(0, nlayers, 1 ):
          layer_depths[i], err = _interpolator.interpolate( lon1 * np.pi / 180.0, lat1 * np.pi / 180.0,
                                        _litho_data[i,l1_data_decode["DEPTH"]], order=1)
          layer_properties[i], err = _interpolator.interpolate( lon1 * np.pi / 180.0, lat1 * np.pi / 180.0,
                                        _litho_data[i,l1_data_decode[quantity_ID]], order=1)
  
  
      A = -layer_depths
      B = -depth1 * 1000.0
      C = divmod(np.searchsorted(A.ravel(), B), A.shape[1])[0] # YEP - this seems to be the best way !!
  
      # point_properties = np.diag(layer_properties[C[:],:])
  
      point_properties = np.empty_like(depth1)
      for i,layer in enumerate(C):
  
          point_properties[i] = layer_properties[layer,i]
  
      return C, point_properties.reshape(shape)
docs: Lat / Lon are in degrees
      Depth in km
      quantity_ID needs to match those in the litho1 model
  
      Points that are not found are given the out-of-range value of -99999