[scautoloc] Add location-dependent depth lookup (Slab2 backend)

apps/processing/scautoloc/CMakeLists.txt

@@ -16,6 +16,7 @@ SET(
 		stationconfig.cpp
 		stationlocationfile.cpp
 		picklog.cpp
+		depthlookup.cpp
 )
 
 SET(
@@ -33,6 +34,7 @@ SET(
 		stationconfig.h
 		stationlocationfile.h
 		picklog.h
+		depthlookup.h
 )
 
 SET(
@@ -50,5 +52,11 @@ SC_LINK_LIBRARIES_INTERNAL(${LOC_TARGET} client)
 SC_INSTALL_DATA(LOC ${LOC_TARGET})
 SC_INSTALL_INIT(${LOC_TARGET} ${INIT_TEMPLATE})
 
+# Install Slab2 depth-contour BNA files
+INSTALL(
+    DIRECTORY share/slabs/
+    DESTINATION share/${LOC_TARGET}/slabs
+)
+
 FILE(GLOB descs "${CMAKE_CURRENT_SOURCE_DIR}/descriptions/*.xml")
 INSTALL(FILES ${descs} DESTINATION ${SC3_PACKAGE_APP_DESC_DIR})

apps/processing/scautoloc/app.cpp

@@ -289,6 +289,18 @@ bool AutolocApp::initConfiguration() {
 	}
 	catch ( ... ) {}
 
+	try {
+		_config.depthLookupType = configGetString("autoloc.depthLookup");
+	}
+	catch ( ... ) {}
+
+	try {
+		_config.slabDir = Environment::Instance()->absolutePath(configGetString("autoloc.slab2.directory"));
+	}
+	catch ( ... ) {
+		_config.slabDir = Environment::Instance()->shareDir() + "/scautoloc/slabs";
+	}
+
 	try {
 		_config.minimumDepth = configGetDouble("locator.minimumDepth");
 	}
@@ -819,7 +831,7 @@ void AutolocApp::done() {
 		ar.setFormattedOutput(_formatted);
 		ar << _outputEP;
 		ar.close();
-		std::cerr << "Output to XML: " << objectCount << " objects(s)" << std::endl;
+		std::cerr << "Output to XML: " << objectCount << " object(s)" << std::endl;
 		_outputEP = nullptr;
 	}
 

apps/processing/scautoloc/autoloc.cpp

@@ -23,6 +23,7 @@
 #include "util.h"
 #include "scutil.h"
 #include "autoloc.h"
+#include "depthlookup.h"
 #include "stationlocationfile.h"
 
 
@@ -38,6 +39,8 @@ using namespace AutolocInternal;
 Autoloc::Autoloc() {
 	_associator.setOrigins(&_origins);
 	_relocator.setMinimumDepth(_config.minimumDepth);
+	// Safe default until init() is called with the actual config
+	_depthLookup = makeDepthLookup("Constant", _config.defaultDepth, _config.maxDepth);
 }
 // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
@@ -81,6 +84,14 @@ bool Autoloc::init() {
 	SEISCOMP_DEBUG("Setting configured locator profile: %s", _config.locatorProfile);
 	setLocatorProfile(_config.locatorProfile);
 
+	_depthLookup = makeDepthLookup(
+	    _config.depthLookupType,
+	    _config.defaultDepth,
+	    _config.maxDepth,
+	    _config.slabDir);
+	SEISCOMP_INFO("DepthLookup: using backend '%s'",
+	              _config.depthLookupType.c_str());
+
 	return true; // ready to start processing
 }
 // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -1320,7 +1331,7 @@ OriginPtr Autoloc::_xxlPreliminaryOrigin(const Pick *newPick)
 	// loop over several trial depths, which are multiples of the default depth
 	std::vector<double> trialDepths;
 	for ( int i = 0; dep <= _config.xxlMaxDepth; i++ ) {
-		dep = _config.defaultDepth*(1+i);
+		dep = _depthLookup->fetch(lat, lon)*(1+i);
 		trialDepths.push_back(dep);
 
 		// in case of "sticky" default depth, we don't need any more trial depths
@@ -1846,17 +1857,31 @@ bool Autoloc::_setDefaultDepth(Origin *origin)
 // May be set in an origin far outside the network where depth resolution is expected to be poor,
 // or in testing that depth resolution.
 {
-	OriginPtr test = new Origin(*origin);
+	auto candidates = _depthLookup->fetchCandidateDepths(
+	    origin->hypocenter.lat, origin->hypocenter.lon);
 
-	_relocator.setFixedDepth(_config.defaultDepth);
-	_relocator.useFixedDepth(true);
-	OriginPtr relo = _relocator.relocate(test.get());
-	if ( !relo ) {
+	OriginPtr bestRelo;
+	double bestScore = -1.0;
+
+	for ( double dep : candidates ) {
+		OriginPtr test = new Origin(*origin);
+		_relocator.setFixedDepth(dep);
+		_relocator.useFixedDepth(true);
+		OriginPtr relo = _relocator.relocate(test.get());
+		if ( !relo ) continue;
+		double score = _score(relo.get());
+		if ( score > bestScore ) {
+			bestScore = score;
+			bestRelo = relo;
+		}
+	}
+
+	if ( !bestRelo ) {
 		SEISCOMP_WARNING("_setDefaultDepth: failed relocation");
 		return false;
 	}
 
-	origin->updateFrom(relo.get());
+	origin->updateFrom(bestRelo.get());
 	origin->depthType = Origin::DepthDefault;
 
 	return true;
@@ -1889,7 +1914,8 @@ bool Autoloc::_setTheRightDepth(Origin *origin) {
 			return false;
 		}
 
-		double radius = 5*(relo->hypocenter.dep >= _config.defaultDepth ? relo->hypocenter.dep : _config.defaultDepth)/111.2;
+		double dd = _depthLookup->fetch(origin->hypocenter.lat, origin->hypocenter.lon);
+		double radius = 5*(relo->hypocenter.dep >= dd ? relo->hypocenter.dep : dd)/111.2;
 
 		// XXX This is a hack, but better than nothing:
 		// if there are at least 2 stations within 5 times the source depth, we assume sufficient depth resolution.
@@ -2079,7 +2105,7 @@ bool Autoloc::_rework(Origin *origin) {
 	}
 
 	if ( enforceDefaultDepth ) {
-		_relocator.setFixedDepth(_config.defaultDepth);
+		_relocator.setFixedDepth(_depthLookup->fetch(origin->hypocenter.lat, origin->hypocenter.lon));
 	}
 
 	bool keepDepth = adoptManualDepth || enforceDefaultDepth;
@@ -2128,7 +2154,7 @@ bool Autoloc::_rework(Origin *origin) {
 	_excludeDistantStations(origin);
 	_excludePKP(origin);
 
-	if ( origin->hypocenter.dep != _config.defaultDepth && origin->depthType == Origin::DepthDefault )
+	if ( origin->hypocenter.dep != _depthLookup->fetch(origin->hypocenter.lat, origin->hypocenter.lon) && origin->depthType == Origin::DepthDefault )
 		origin->depthType = Origin::DepthFree;
 
 	// once more (see also above)
@@ -2360,9 +2386,10 @@ bool Autoloc::_publishable(const Origin *origin) const
 	}
 
 
-	if ( origin->hypocenter.dep > _config.maxDepth ) {
+	double maxDep = _depthLookup->fetchMaxDepth(origin->hypocenter.lat, origin->hypocenter.lon);
+	if ( origin->hypocenter.dep > maxDep ) {
 		SEISCOMP_INFO("Origin %ld too deep: %.1f km > %.1f km (maxDepth)",
-			      origin->id, origin->hypocenter.dep, _config.maxDepth);
+			      origin->id, origin->hypocenter.dep, maxDep);
 		return false;
 	}
 
@@ -2412,7 +2439,7 @@ bool Autoloc::_store(Origin *origin)
 	}
 
 	if ( origin->depthType == Origin::DepthDefault &&
-	    origin->hypocenter.dep != _config.defaultDepth ) {
+	    origin->hypocenter.dep != _depthLookup->fetch(origin->hypocenter.lat, origin->hypocenter.lon) ) {
 		origin->depthType = Origin::DepthFree;
 	}
 
@@ -2506,7 +2533,7 @@ bool Autoloc::_associate(Origin *origin, const Pick *pick, const std::string &ph
 			bool fixed = false;
 			if ( _config.defaultDepthStickiness > 0.9 ) {
 				fixed = true;
-				_relocator.setFixedDepth(_config.defaultDepth);
+				_relocator.setFixedDepth(_depthLookup->fetch(origin->hypocenter.lat, origin->hypocenter.lon));
 			}
 
 //			else if ( origin->depthType == Origin::DepthManuallyFixed || origin->depthType == Origin::DepthPhases ) {
@@ -3497,7 +3524,7 @@ bool Autoloc::_depthIsResolvable(Origin *origin) {
 //		return true;
 //	}
 
-	if ( origin->depthType == Origin::DepthDefault && origin->hypocenter.dep != _config.defaultDepth )
+	if ( origin->depthType == Origin::DepthDefault && origin->hypocenter.dep != _depthLookup->fetch(origin->hypocenter.lat, origin->hypocenter.lon) )
 		origin->depthType = Origin::DepthFree;
 
 	OriginPtr test = new Origin(*origin);
@@ -3516,7 +3543,7 @@ bool Autoloc::_depthIsResolvable(Origin *origin) {
 	}
 
 	test = new Origin(*origin);
-	test->hypocenter.dep = _config.defaultDepth;
+	test->hypocenter.dep = _depthLookup->fetch(origin->hypocenter.lat, origin->hypocenter.lon);
 	_relocator.useFixedDepth(true);
 	relo = _relocator.relocate(test.get());
 	if ( !relo ) {

apps/processing/scautoloc/autoloc.h

@@ -27,6 +27,7 @@
 #include "config.h"
 #include "stationconfig.h"
 #include "picklog.h"
+#include "depthlookup.h"
 
 
 namespace Seiscomp {
@@ -364,6 +365,7 @@ class Autoloc {
 		AutolocInternal::OriginVector _origins;
 		AutolocConfig _config;
 		AutolocInternal::StationConfig _stationConfig;
+		DepthLookupPtr _depthLookup;
 
 		const Seiscomp::Config::Config *scconfig {nullptr};
 		const Seiscomp::DataModel::Inventory *scinventory {nullptr};

apps/processing/scautoloc/config.cpp

@@ -44,7 +44,11 @@ void AutolocConfig::dump() const {
 	SEISCOMP_INFO("    maxResidual                      %.1f s", maxResidualUse);
 	SEISCOMP_INFO("    maxResidual for keeping picks    %.1f s", maxResidualKeep);
 	SEISCOMP_INFO("    minPhaseCount                    %d",     minPhaseCount);
-	SEISCOMP_INFO("    maxDepth                         %.1f km", maxDepth);
+	SEISCOMP_INFO("    depthLookup                      %s",     depthLookupType.c_str());
+	if ( depthLookupType == "Slab2" )
+		SEISCOMP_INFO("    slab2.directory                  %s",     slabDir.c_str());
+	else
+		SEISCOMP_INFO("    maxDepth                         %.1f km", maxDepth);
 	SEISCOMP_INFO("    minStaCountIgnorePKP             %d",     minStaCountIgnorePKP);
 	SEISCOMP_INFO("    defaultDepthStickiness           %g",     defaultDepthStickiness);
 	SEISCOMP_INFO("    tryDefaultDepth                  %s",     tryDefaultDepth ? "true" : "false");

apps/processing/scautoloc/config.h

@@ -91,7 +91,18 @@ struct AutolocConfig {
 		//    maxResidualUse  = 2*maxRMS
 
 
-		// Use this depth if there is no depth resolution
+		// DepthLookup backend: "Constant" (default) or "Slab2"
+		// "Constant" uses defaultDepth/maxDepth below — same as previous behaviour
+		std::string depthLookupType{"Constant"};
+
+		// Directory containing Slab2 BNA depth-contour files.
+		// Only used when depthLookupType = "Slab2".
+		// Default: @DATADIR@/spatial/vector/slabs
+		std::string slabDir;
+
+		// Use this depth if there is no depth resolution.
+		// Used directly by the Constant backend; overridden geographically
+		// by Polygon and Slab2 backends.
 		double defaultDepth{10.0};           // unit: km
 		double defaultDepthStickiness{0.5};  // 0...1
 
@@ -105,7 +116,9 @@ struct AutolocConfig {
 		// Minimum depth in case there is depth resolution
 		double minimumDepth{5.0};          // uni: 5 km
 
-		// maximum depth of origin, checked before sending
+		// maximum depth of origin, checked before sending.
+		// Used directly by the Constant backend; overridden geographically
+		// by Polygon and Slab2 backends.
 		double maxDepth{1000.0};
 
 		// Max. secondary azimuthal gap