Такие гифки видели все: город из космоса, по кадру в год. Серое пятно с каждым кадром больше. Google Earth Timelapse в прошлом году выложил такие таймлапсы для всей Земли. Я всегда верил тому, что они показывают. Ни разу не спросил, сколько в этом движении собственно города. Этой осенью я захотел такую же гифку для Минска, 10 лет по Sentinel-2, потому что его данные бесплатные.
10 лет не получилось. Sentinel-2A запустили в июне 2015 года. Первая сцена над городом в архиве датирована 29 июля 2015 года. В декабре 2022 года в архиве 7 полных летних сезонов, с 2016 по 2022. Для 10 лет понадобился бы Landsat с разрешением 30 метров.
Первая гифка
Я взял сцены Sentinel-2 L2A с Planetary Computer из тайла 35UNV. Окно около 32 на 28 км, в нём город вместе с МКАД и пригороды. Для каждого лета, с мая по сентябрь, я выбрал сцену с наименьшей облачностью в окне и собрал из 7 кадров гифку.
Это было похоже на растущий город: поля прыгали между бурым и зелёным, водохранилище на северо-западе меняло цвет. Потом я попробовал это посчитать. Пиксель “голый”, если его индекс растительности (NDVI) меньше 0,25 и это не вода. Между кадром 22 августа 2016 года и кадром 23 мая 2022 года голыми стали 64 км². За то же время 51 км² перестал быть голым, так что чистый прирост 13 км² к 108. Если смотреть только на прирост, как это делает глаз на гифке, голая площадь за неполные 6 лет выросла в полтора с лишним раза.
Я посмотрел, что лежит под этими 64 км² на карте ESA WorldCover за 2021 год. 60% из них пашня. В мае на многих полях ещё ничего не выросло, поэтому они выглядят как голая почва. Мой индекс не отличает их от новостройки. В более строгий счёт ниже из этих 64 км² попало только 5%.
Медианы вместо одиночных сцен
Обычно это лечат медианой. Для каждого лета я взял до 12 самых чистых сцен и замаскировал облака и тени по слою классификации сцены (SCL) из L2A. В каждом пикселе остаётся медиана его значений за лето. Потом я посчитал пиксели, которые меняют состояние, голый или не голый, хотя бы 2 раза за 7 кадров. По одиночным сценам это 157 км², по медианам 59 км², что тоже немало. Для одиночных сцен счёт идёт по меньшей площади, потому что у них больше пикселей под маской.
Голая площадь по медианам каждый год держится между 81 и 87 км² без всякого тренда, больше всего её в 2018 году.
В данных 2022 года была ещё одна ловушка. С 25 января 2022 года ESA обрабатывает сцены версией обработки (processing baseline) 04.00, которая прибавляет 1000 к каждому значению коэффициента отражения. Я эту добавку вычел. Без этого то же правило даёт 154 км² голой земли в 2022 году вместо 83 км².
Правило: смена должна продержаться
Тогда я стал требовать больше одного кадра. Пиксель считается новым городом, если он был зелёным (NDVI 0,5 и выше) в 2016 и 2017 годах, а голым в 2021 и 2022. Ещё он должен лежать в пятне хотя бы из 5 пикселей. Получается 3,46 км² в 854 пятнах. Обратная смена, из голого в зелёное при том же правиле пятен, даёт 0,51 км².
Этому числу я тоже не поверил. Поэтому я выбрал из этих пятен 120 случайных пикселей и посмотрел на каждый в медианах 2016, 2019 и 2022 годов. Я думал, что поля и лесные вырубки составят хотя бы половину. Я ошибся. На 77 пикселях к 2022 году стояло здание, дорога или асфальтированный двор. 33 оказались стройплощадками или землёй, расчищенной под стройку. 1 был полем, а 9 я не смог разобрать при разрешении 10 метров, в основном в частном секторе.
20 пикселей из 120 попали в одно место, на бывший аэропорт Минск-1. Его закрыли 23 декабря 2015 года, так что летом 2016 года, первым для Sentinel-2, полоса ещё на месте. В медианах она лежит в зелёном поле до 2018 года, в 2019 году появляется песок. К 2022 году на северо-западной половине стоят жилые дома. В радиусе 1 км от середины полосы лежит 0,54 км² изменившейся площади, 15% всех изменений в окне. Самое большое отдельное пятно тоже там.
Сверка с готовыми картами
Ещё до счёта я собирался сравнить результат с 2 готовыми картами. Первая это GHSL Еврокомиссии, выпуск R2022A, вышедший этим летом, застроенная поверхность с шагом 5 лет. На моей сетке она даёт 73,9 км² застройки в 2015 году и 77,2 км² в 2020 году, то есть на 3,3 км² больше. Мои 3,46 км² с 2016 по 2022 год почти столько же. Я этому обрадовался.
Потом я посмотрел, где находятся эти 3,3 км². 7% прироста GHSL приходится на полосу 100 метров вокруг моих пятен. Эта полоса занимает 3,4% окна, то есть совпадение примерно вдвое больше случайного. И всё равно остальные 93% в других местах. На клетках по 500 метров корреляция двух карт 0,06. На территории аэропорта GHSL почти не меняется: на моих изменившихся пикселях доля застройки у неё 4,6% в 2015 году и 4,9% в 2020. Тогда я наконец внимательно прочитал описание GHSL. Настоящих наблюдений в ней 2: Landsat 2014 года и композит Sentinel-2 за 2018 год. Эпоху 2015 года GHSL интерполирует, эпоху 2020 года экстраполирует. Для обеих работает модель, которая учитывает среди прочего рельеф, расстояние до воды и уже наблюдавшийся рост. Так что её число за 2020 год это не наблюдение.
С WorldCover для этого вопроса ещё хуже. Карта 2020 года сделана версией алгоритма v100, а карта 2021 года версией v200. В примечаниях к выпуску авторы предупреждают, что разница между двумя картами складывается из изменений на земле и из смены алгоритма. В моём окне застроенный класс за один год растёт со 153,9 до 186,8 км², а моё правило даёт меньше 1 км² в год.
Что я из этого вынес
В моём окне вокруг Минска нового города около 3,5 км² за 6 лет вместе со стройплощадками. Думаю, настоящее число скорее больше, чем меньше. Пиксель должен был быть зелёным 2 лета подряд перед сменой. Застройка мелкими кусками не проходит правило пятен. Но сколько правило пропустило, я не измерял.
Чего я не проверял
Пороги 0,25 и 0,5 я выбрал один раз и не подбирал. Земля, расчищенная до 2016 года, в счёт не попала, а в окне есть ещё и деревни за городом. Выборку я размечал один при разрешении 10 метров, другой человек мог бы отнести несколько точек в другую группу. Где старые дома с садами сменились новыми кварталами, моё правило посчитало бы это новым городом, хотя земля была застроена и раньше. Более высокое здание на месте низкого оно не видит совсем.
На 2,2% окна хотя бы за одно лето не нашлось ни одного чистого наблюдения. В счёт они не вошли. Сначала я думал, что маска принимает светлые крыши за облака. Классы говорят другое: около половины этих пикселей “unclassified” и треть “dark area”, 2 класса, которые выбросила моя собственная маска. В гифках я заполнил дыры значениями ближайшего пикселя.