code: def property_on_depth_profile(lat, lon, depths, quantity_ID="DENSITY"):
      
  
      lon1 = np.array((lon,))
      lat1 = np.array((lat,))
      depths1 = np.array(depths)
  
      nlayers = len(l1_layer_decode)
      point_properties = np.empty_like(depths)
  
      layer_depths     = np.empty((nlayers))
      layer_properties = np.ones((nlayers+1)) * -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( np.radians(lon1), np.radians(lat1),
                                        _litho_data[i,l1_data_decode["DEPTH"]], order=1)
          layer_properties[i], err = _interpolator.interpolate( np.radians(lon1), np.radians(lat1),
                                        _litho_data[i,l1_data_decode[quantity_ID]], order=1)
  
  
      A = -layer_depths
      B = -depths1 * 1000.0
      C = np.searchsorted(A, B)
  
      point_properties = np.empty_like(depths)
      for i,layer in enumerate(C):
  
          point_properties[i] = layer_properties[layer]
  
      return C, point_properties
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