code: def layer_depth( lat, lon, layerID="LID-BOTTOM"):
      
  
      ## Must wrap longitude from 0 to 360 ...
  
      lon1 = np.array(lon)%360.0
      lat1 = np.array(lat)
  
      # ## Must wrap longitude from -180 to 180 ...
      #
      # lon1[np.where(lon1 > 180.0)] = 360.0 - lon1[np.where(lon1 > 180.0)]
      #
      data, err = _interpolator.interpolate( np.radians(lon1), np.radians(lat1),
                                        _litho_data[l1_layer_decode[layerID], l1_data_decode["DEPTH"]], order=1 )
  
      return data
docs: Returns layer depth at lat / lon (degrees)
      where lat/lon may be arrays (of equal size).
      Depths are returned in metres.