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

def property_on_depth_profile(lat, lon, depths, quantity_ID="DENSITY"):
    """
    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

    """

    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