180-й меридиан

georef-lrpt

Меня всегда поражало, насколько довольно простые с виду вещи, в реальности могут оказаться достаточно запутанными.

В далёком 2020 году я описывал процесс наложения данных, полученных со спутника Метеор-М2 на интерактивную карту. Эта задача достаточно сложная и запутанная, поэтому я решил оформить весь алгоритм в виде небольшой консольной утилиты - georef-lrpt. Она создаёт .png изображение и .vrt файл на основе бинарных данных полученных со спутника.

Запускается программа следующим образом:

java -jar georef-lrpt.jar --output-dir . --tle-file src/test/resources/2022-12-20.txt --vcdu-files "src/test/resources/*.vcdu"

На вход подаётся TLE и список бинарных файлов .vcdu. Внутри бинарных файлов должны быть VCDU-фреймы в том порядке, в котором были получены со спутника. На выходе получится два файла: 2022_12_20_20_25_30.png и 2022_12_20_20_25_30.vrt.

Второй файл - это специальный файл для GDAL в формате .vrt в котором хранится геопривязка пикселей. Этот файл можно использовать для того, чтобы создать GeoTIFF:

gdalwarp -tps -overwrite -of GTIFF 2022_12_20_20_25_30.vrt 2022_12_20_20_25_30.tiff

Который потом можно наложить на карту:

Выглядит идеально, но я решил скачать ещё данных, чтобы потестировать на граничных значениях и обнаружил интересное:

А что это такое в центре Атлантического океана? А это - несколько дней дебага!

180-ый меридиан

Если приглядется, то можно увидеть, что небольшой кусочек изображения отображается слева. Видимо он перенёсся с другой стороны из-за того, что Земля круглая. А вот с артефактом посередине немного сложнее. Если выключить слой прозрачности, то будет лучше видно:

Добро пожаловать в проблему 180-го меридиана или, как её называют во всём остальном мире - antimeridian issue.

Некоторые изображения могут начинаться от 170 меридиана и заканчиваться в -170. А некоторые начинаться в -170 и заканчиваться в 170-ом. Видимо GDAL не знает какой случай где и отображает и так, и эдак. Я попробовал искать опцию “wrapOn180=true”, которая бы говорила, что разрез приходится на 180-ый меридиан, но почему-то не нашёл ничего похожего в документации. На форумах советуют разрезать изображение на две части вручную.

Алгоритм

Итак, первым делом нужно определить границы изображения. Это делается достаточно легко:

if (y > top) {
	top = y;
}
if (y < bottom) {
	bottom = y;
}
if (!wrapDetected && previousX != UNDEFINED) {
	if (x < 0 && previousX > 170.0) {
		// going north
		wrapDetected = true;
		right = -180.0;
	}
	if (x > 170.0 && previousX < 0) {
		wrapDetected = true;
		left = 180;
	}
}
if (wrapDetected) {
	if (x < 0) {
		right = Math.max(x, right);
	}
	if (x > 0) {
		left = Math.min(x, left);
	}
} else {
	left = Math.min(x, left);
	right = Math.max(x, right);
}
previousX = x;

Где x и y - это longitude и latitude текущей GCP точки.

Далее нужно каким-то образом передать полученные координаты в .vrt, чтобы внешние скрипты смогли понять нужно ли делать два изображения или достаточно сделать одно. Если left > right , то произошёл перенос и нужно делать два.

И тут меня постигла очередная неудача. Я не смог найти ни одного способа, чтобы задать такие параметры в .vrt. Можно было попробовать добавить мета-поля, но gdalwarp их будет игнорировать. Поэтому я придумал ничего лучше, чем подставить эти значения в название файла. Поэтому вместо 2022_12_20_18_44_50.vrt у меня создаётся файл 2022_12_20_18_44_50_142_-55_-173_-11.vrt. Внешние скрипты должны распарсить это имя и получить настоящие границы изображения. После этого создать два GeoTIFF файла следующим образом:

gdalwarp -tps -overwrite -te 142 -55 180 -11 -of GTiff 2022_12_20_18_44_50_142_-55_-173_-11.vrt 2022_12_20_18_44_50-left.tiff
gdalwarp -tps -overwrite -te -180 -55 -173 -11 -of GTiff 2022_12_20_18_44_50_142_-55_-173_-11.vrt 2022_12_20_18_44_50-right.tiff

Если добавить левое и правое изображение на карту, то можно получить нужный результат:

А если показать в проекции EPSG:3832 ( PDC Mercator ) и отключить слой прозрачности, то можно увидеть как изображения склеиваются: