int GetTimezoneOffset(double lon, double lat, const char* shapefile_path, std::string* out_tzid) {
// 1. 坐标标准化和转换 (核心修复)
lon = std::fmod(lon, 360.0);
if (lon > 180) lon -= 360;
else if (lon < -180) lon += 360;
// 转换为shapefile期望的坐标顺序 (纬度, 经度)
const double point_coords[2] = {lat, lon}; // [Y, X] 格式
// 2. GDAL初始化
GDALAllRegister();
CPLSetConfigOption("SHAPE_RESTORE_SHX", "YES");
// 方案B: 不使用exportToWkt的替代方法
GDALDataset* poDS = static_cast<GDALDataset*>(
GDALOpenEx(shapefile_path, GDAL_OF_VECTOR | GDAL_OF_READONLY,
nullptr, nullptr, nullptr));
if (!poDS) {
FLOGW << "Open failed: " << CPLGetLastErrorMsg();
return 1;
}
std::unique_ptr<GDALDataset, decltype(&GDALClose)> ds_guard(poDS, GDALClose);
OGRLayer* poLayer = poDS->GetLayer(0);
if (!poLayer) {
FLOGW << "Layer 0 not found";
return 1;
}
// 获取图层的实际SRS
const OGRSpatialReference* layerSRS = poLayer->GetSpatialRef();
// 创建查询点 (使用图层SRS)
OGRPoint point;
point.setX(lon); // 经度
point.setY(lat); // 纬度
if (layerSRS) {
point.assignSpatialReference(layerSRS);
// 调试输出SRS信息
char* srsWkt = nullptr;
layerSRS->exportToPrettyWkt(&srsWkt);
FLOGD << "Layer SRS: " << (srsWkt ? srsWkt : "NULL");
CPLFree(srsWkt);
}
// 7. 精确点查询 (完全模拟ogrinfo行为)
poLayer->SetSpatialFilter(&point);
poLayer->ResetReading();
// 8. 收集所有匹配的时区ID
std::vector<std::string> matched_tzids;
OGRFeature* poFeature;
while ((poFeature = poLayer->GetNextFeature()) != nullptr) {
std::unique_ptr<OGRFeature, decltype(&OGRFeature::DestroyFeature)> feature_guard(
poFeature, OGRFeature::DestroyFeature);
OGRGeometry* geom = poFeature->GetGeometryRef();
if (!geom) continue;
// 统一空间参考系统
if (geom->getSpatialReference() == nullptr) {
geom->assignSpatialReference(poLayer->GetSpatialRef());
}
// 精确包含测试
if (geom->Contains(&point) || geom->Intersects(&point)) {
const char* tzid = poFeature->GetFieldAsString("tzid");
if (tzid) {
matched_tzids.push_back(tzid);
FLOGD << "Found matching TZ: " << tzid;
}
}
}
// 9. 验证结果 (关键调试步骤)
if (matched_tzids.empty()) {
FLOGW << "NO MATCHING TIMEZONE FOUND for (" << lon << "," << lat << ")";
// 备用策略：使用ogrinfo方法
FLOGW << "Falling back to bbox query method";
OGREnvelope env;
point.getEnvelope(&env);
poLayer->SetSpatialFilterRect(env.MinX, env.MinY, env.MaxX, env.MaxY);
poLayer->ResetReading();
while ((poFeature = poLayer->GetNextFeature()) != nullptr) {
std::unique_ptr<OGRFeature, decltype(&OGRFeature::DestroyFeature)> feature_guard(
poFeature, OGRFeature::DestroyFeature);
OGRGeometry* geom = poFeature->GetGeometryRef();
if (!geom) continue;
if (geom->Contains(&point)) {
const char* tzid = poFeature->GetFieldAsString("tzid");
if (tzid) {
matched_tzids.push_back(tzid);
FLOGD << "Fallback match: " << tzid;
}
}
}
}
// 10. 结果处理和时区选择
if (matched_tzids.empty()) {
FLOGW << "No timezone found for point (" << lon << "," << lat << ")";
return 1;
}
for(auto& id:matched_tzids) {
FLOGI << "tony ids:" << id;
}
// 选择策略：优先大陆时区
std::string selected_tzid;
auto it = std::find_if(matched_tzids.begin(), matched_tzids.end(),
[](const std::string& tz) {
return tz.find("Asia") != std::string::npos ||
tz.find("Europe") != std::string::npos;
});
selected_tzid = (it != matched_tzids.end()) ? *it : matched_tzids[0];
FLOGI << "Matched timezone: " << selected_tzid << " at (" << lon << "," << lat << ")";
// 11. 时区偏移计算 (跨平台兼容)
const char* orig_tz = getenv("TZ");
std::unique_ptr<char, void(*)(void*)> old_tz(
orig_tz ? strdup(orig_tz) : nullptr, [](void* p) { if(p) free(p); });
setenv("TZ", selected_tzid.c_str(), 1);
tzset();
time_t now = time(nullptr);
struct tm local_tm;
#ifdef _WIN32
localtime_s(&local_tm, &now);
#else
localtime_r(&now, &local_tm);
#endif
int offset = 0;
#if defined(__linux__) || defined(__APPLE__)
offset = static_cast<int>(local_tm.tm_gmtoff / 60);
#elif defined(_WIN32)
TIME_ZONE_INFORMATION tz_info;
GetTimeZoneInformation(&tz_info);
offset = -static_cast<int>(tz_info.Bias);
#else
struct tm utc_tm;
gmtime_r(&now, &utc_tm);
time_t utc_time = mktime(&utc_tm);
time_t local_time = mktime(&local_tm);
offset = static_cast<int>(difftime(local_time, utc_time) / 60);
#endif
// 恢复原始TZ设置
if (old_tz) {
setenv("TZ", old_tz.get(), 1);
} else {
unsetenv("TZ");
}
tzset();
if (out_tzid) *out_tzid = selected_tzid;
return offset;
}
其中时区库是combined-shapefile-with-oceans-now.shp2023版本，gdal库是3.7.0版本，geom->getGeometryType()的值为0，也就是wkbMultiPolygon类型。现在我在不改变函数功能的前提下，在代码功能的基础上增加一个“面积最小优先"优先级逻辑”，请你帮我完成。
博客内容为 introduction，未提供更多关键信息。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
470
