NOAA提取磁偏角数据

最近开发定位跟踪系统，计算需要用到磁偏角数据来校正磁北和真北方向，借鉴了px4源码中的利用计算航向调用磁偏角的方法，从NOAA Geomagnetic Field Calculator中去提取磁偏角。现简单记录下方法，由于地球磁场长期变化的，需要定时维护。

打开该网站后，选择“Magnetic Field Component Grid”计算和导出磁偏角。

Southern most lat&&Nourthern  most lat：纬度坐标设置在-60~60°之间，其他区域都在南北极了，一般人应该也用不到，可以看个人需求。

Western most long:&&Eastern most long:经度-180~180°之间，这没什么好说的，或者看个人需求。

Lat Step Size&&Lon Step Size:步进精度选择10°，看个人需求。

Elevation:：海报高度选择GPS  0km,也可以看个人需求。

“Magnetic component”磁性元件，选择磁偏角“Declination”，如果有其他磁倾角、磁场等需求也可以选择其他项。

Model:选择IGRF模型或者WMM都可以。

Start Date&&End Date&&Step size 默认选择当天即可。最后计算我是选择csv导出格式，通过python导出 Declination_sv in Degree列数据索引数组，代码如下：

import csv


with open('C:\\Users\\yanfabu\\Desktop\\igrfgridData.csv', mode='r') as file:
    reader = csv.reader(file)
    data = list(reader)


filtered_data = []
for row in data:
    if row[0].startswith('#') or row[0].startswith('Year'):
        continue
    filtered_data.append(row)


declination_table = []


for row in filtered_data:
    try:
        lat = float(row[1])
        lon = float(row[2])
        declination = float(row[4])


        print(f"Lat: {lat}, Lon: {lon}, Declination: {declination}")


        if len(declination_table) == 0 or len(declination_table[-1]) == 37:
            declination_table.append([declination])
        else:
            declination_table[-1].append(declination)
    except ValueError:
        continue


print("static const float declination_table[13][37] = {")
for row in declination_table:
    print("    {" + ", ".join(map(str, row)) + "},")
print("};")

借鉴px4 “geo_mag_declination”查找表方法，获取到该坐标的磁偏角校准真北方向。借鉴代码如下：

get_lookup_table_index(float *val, float min, float max)
{
	/* for the rare case of hitting the bounds exactly
	 * the rounding logic wouldn't fit, so enforce it.
	 */

	/* limit to table bounds - required for maxima even when table spans full globe range */
	/* limit to (table bounds - 1) because bilinear interpolation requires checking (index + 1) */
	*val = constrain(*val, min, max - SAMPLING_RES);

	return static_cast<unsigned>((-(min) + *val) / SAMPLING_RES);
}

static float
get_table_data(float lat, float lon, const int8_t table[13][37])
{
	/*
	 * If the values exceed valid ranges, return zero as default
	 * as we have no way of knowing what the closest real value
	 * would be.
	 */
	if (lat < -90.0f || lat > 90.0f ||
	    lon < -180.0f || lon > 180.0f) {
		return 0.0f;
	}

	/* round down to nearest sampling resolution */
	float min_lat = floorf(lat / SAMPLING_RES) * SAMPLING_RES;
	float min_lon = floorf(lon / SAMPLING_RES) * SAMPLING_RES;

	/* find index of nearest low sampling point */
	unsigned min_lat_index = get_lookup_table_index(&min_lat, SAMPLING_MIN_LAT, SAMPLING_MAX_LAT);
	unsigned min_lon_index = get_lookup_table_index(&min_lon, SAMPLING_MIN_LON, SAMPLING_MAX_LON);

	const float data_sw = table[min_lat_index][min_lon_index];
	const float data_se = table[min_lat_index][min_lon_index + 1];
	const float data_ne = table[min_lat_index + 1][min_lon_index + 1];
	const float data_nw = table[min_lat_index + 1][min_lon_index];

	/* perform bilinear interpolation on the four grid corners */
	const float lat_scale = constrain((lat - min_lat) / SAMPLING_RES, 0.0f, 1.0f);
	const float lon_scale = constrain((lon - min_lon) / SAMPLING_RES, 0.0f, 1.0f);

	const float data_min = lon_scale * (data_se - data_sw) + data_sw;
	const float data_max = lon_scale * (data_ne - data_nw) + data_nw;

	return lat_scale * (data_max - data_min) + data_min;
}

float get_mag_declination(float lat, float lon)
{
	return get_table_data(lat, lon, declination_table);
}