Услуга Google Планета Земля предоставляет возможность бесплатно работать с огромными объемами пространственной информации.
Например, за считанные минуты из миллиона космических снимков можно получить составную мозаику (составное изображение).
Если предположить, что каждая сцена Landsat 8 (набор спектральных каналов) занимает 1 ГБ в сжатом виде, то такой запрос обрабатывает объем информации порядка 1 ПБ.
И все это доступно бесплатно, быстро и в любое время.
Но есть мнение (ошибочное), что GEE на бесплатных аккаунтах позволяет обрабатывать и экспортировать только небольшие наборы данных.
На самом деле такое впечатление создается лишь тем, что начать программировать в GEE можно, даже не читая сервисную документацию, но извлечь много данных, не читая документацию, не получится.
Далее мы рассмотрим три различных решения проблемы векторизации растра и напишем серверную функцию GEE для вычисления геохеша двумя разными способами.
Введение
Нравится ли вам работать на суперкомпьютерах (например, университетских кластерах)? Мне не! Несмотря на всю их мощь, на практике полезность таких монстров зачастую совсем не очевидна.Можете упрекнуть меня в предвзятости (и это правда — для меня Linux Debian, конечно, лучший, да и Ubuntu, как его незадачливый потомок, тоже на что-то годится), но работа над Linux Centos, а это около десяти лет старый, совсем перебор (он еще новый) всегда второй по свежести, если посмотреть версии пакетов в нем).
Однако это все равно цветы.
Чтобы не тратить драгоценные ресурсы суперкомпьютеров, на них обычно установлен коммерческий компилятор Intel, причем очень старый (может быть старше самого Centos, который, судя по всему, уже когда-то обновлялся).
Но это не проблема, в конце концов.
Планировщики ресурсов в таких системах часто вызывают путаницу, и не зря — при попытке запустить задание, скажем, на 128 ядрах ничего не происходит, ни задание не запускается, ни утилита запуска не выдает никаких ошибок.
В процессе оказывается, что количество доступных нам ядер равно двум (всегда удивляюсь, почему по умолчанию дают так мало ресурсов), увеличение квоты требует административной работы (часто на этом этапе проще переключиться на какое-то Амазон AWS).
Но самое худшее происходит позже.
С завидной последовательностью бинарные модули вышеупомянутого проприетарного компилятора в конечном итоге устанавливаются в домашний каталог пользователя, который уже давно исчез вместе с его домашним каталогом, и без них компилятор, естественно, ничего не компилирует. Другого компилятора в системе, конечно, нет, потому что его код не был бы таким оптимальным.
Таким образом, ситуация проясняется - нам нужно найти подходящий Centos, установить его где-нибудь локально или в облаке, портировать необходимое программного обеспечения к нему, делаем статические сборки (поскольку на целевой системе некоторые системные пакеты могут быть заменены на что-то неизвестно и откуда) и тогда, возможно, мы сможем запустить все это на двух доступных нам ядрах кластера (невозможно предсказать, будут ли когда-либо доступны новые из них по разным причинам).
Я мог бы еще многое добавить, но вы уже понимаете, почему лично я не люблю суперкомпьютеры.
Теперь поговорим о другом суперкомпьютере — облаке и госсервисе.
Все, что нам нужно для начала, это открыть ссылку в браузере и можно начинать писать код на javascript и загружать данные (кстати, туда же можно загружать свои данные, а также можно использовать API сервиса для доступа из других разработок).
среду и языки).
В то же время сразу доступен огромный объем данных и вычислительных ресурсов.
На мой взгляд, все это стоит того, чтобы прочитать документацию и научиться достаточно эффективно работать в сервисе.
Однако обработка и получение растров вполне очевидна, тем более, что все наборы данных, представленные в сервисе, также имеют примеры их визуализации, и чтобы загрузить их на Google Диск (например) в растровом формате (обычно GeoTIFF), есть функция Экспорта.
достаточно .
image.toDrive. Но работа с векторными данными, особенно преобразование растров в векторы, далеко не так очевидна, и в Интернете часто можно встретить жалобы на то, что невозможно загрузить несколько сотен тысяч записей.
Специально сейчас проверил — за 21 минуту мне удалось загрузить на Google Диск 10 миллионов записей в формате GeoJSON, получив файл размером 2ГБ.
И это далеко не предел для бесплатного аккаунта (не говоря уже о том, что можно загружать данные частями).
Преобразование растровых данных сервиса в точечные данные
Рассмотрим три подхода:- Самый простой способ — создать геометрию, включающую необходимые точки, и получить для них атрибуты растра:
Метод подходит в основном при работе с векторными данными, которые мы хотим дополнить информацией из растров.var points_with_attributes = Image.reduceRegion({ reducer: ee.Reducer.toList(Image.bandNames().size()), geometry: points })
Здесь мы априори знаем координаты всех точек и получаем для них атрибуты растра.
- Более интересный метод предполагает преобразование всех растровых пикселей, попадающих в заданную область, в список.
Давайте сначала извлечем только атрибуты растровых пикселей:
var attributes = Image.reduceRegion({
reducer: ee.Reducer.toList(Image.bandNames().
size()),
geometry: some_area
})
Метод работает, но имеет подвох и ресурсозатратен.
В чем подвох? Дело в том, что в некоторых точках растра атрибутивные данные могут быть пропущены и такие точки в данном методе исключаются.
То есть метод возвращает от нуля до полного набора точек, попадающих в указанную область.
И как только мы захотим получить координаты возвращаемых точек, метод может перестать работать (при наличии пробелов в данных, как было отмечено выше).
Давайте посмотрим на пример кода: var latlon = ee.Image.pixelLonLat()
var points_with_attributes = image
.
addBands(latlon) .
reduceRegion({
reducer: ee.Reducer.toList(),
geometry: some_area
});
Функция ee.Image.pixelLonLat() создает растр с координатами каждой точки, который мы добавляем в рабочий растр.
Этот код прекрасно работает, но только до тех пор, пока в растровых данных нет пробелов и совпадают длины списков координат и атрибутов (по возможности все операции сервисом GEE выполняются параллельно и, как следствие, список для каждого канала формируется независимо).
Кроме того, необходимо добавить перепроецирование данных, но эта операция ресурсоемкая и сильно ограничивает наши возможности по обработке данных на бесплатном аккаунте: var latlon = ee.Image.pixelLonLat().
reproject(image.projection()) var points_with_attributes = image .
addBands(latlon) .
reduceRegion({
reducer: ee.Reducer.toList(),
geometry: some_area
});
К счастью, можно решить обе проблемы сразу: var latlon = ee.Image.pixelLonLat()
var points_with_attributes = image
.
addBands(latlon) .
add(image.select(0).
add(0)) .
reduceRegion({
reducer: ee.Reducer.toList(),
geometry: some_area
});
Здесь мы проводим бесполезные вычисления с данными (при условии, что первый канал растра содержит числовые данные), чтобы унифицировать длину списков атрибутов и списков координат, и в то же время выполняем перепроекцию координатного растра.
В результате все координаты, для которых не определены атрибуты, будут удалены (это может уменьшить количество полученных точек во много раз), координатный растр будет спроецирован на атрибутный растр, и данные будут обработаны корректно.
Это трюк, основанный на парадигме ленивых (ленивых) вычислений сервиса.
- И, наконец, самый простой и быстрый способ, который позволит обрабатывать большое количество данных:
var latlon = ee.Image.pixelLonLat()
var points_with_attributes = image
.
sample({
region: some_area,
geometries: true
});
В приведенном выше коде флаг «геометрия» указывает, нужны ли нам координаты точек или только их атрибуты.
Создание собственных функций сервера
Обычно помимо получения данных «как есть», мы еще хотим их как-то обработать.Это можно сделать локально, загрузив данные, или непосредственно в сервисе GEE. Во втором случае мы получаем огромные вычислительные ресурсы, но нам необходимо писать специальный код для серверной обработки.
Следует помнить, что в окне облачного редактора мы можем писать смесь клиентского и серверного кода (что часто приводит к различным ошибкам).
Когда вы впервые начинаете писать серверные функции GEE, нужны рабочие примеры, но их трудно найти.
Например, напишем функцию расчета геохеш или хэш z-кривой.
К слову, мне не удалось найти готовую реализацию такой функции, что еще раз подтверждает, насколько сложно новичкам найти примеры серверных функций GEE. Чем эта функция полезна? Помимо прочего, геохэш позволяет очень просто и без операций с геометрическими объектами производить объединение пространственных данных — нужно лишь уменьшить количество символов геохеша.
Кодирование и декодирование геохешей поддерживается в PostgreSQL/PostGIS, Google BigQuery и многих других системах.
Итак, давайте посмотрим на простую реализацию функции в стиле кода C: var geohash_accumulate = function(obj, list) {
// define common variables
var base32 = ee.String('0123456789bcdefghjkmnpqrstuvwxyz'); // (geohash-specific) Base32 map
// get previous state variables
var prev = ee.Dictionary(ee.List(list).
get(-1)); var lat = ee.Number(prev.get('lat',0)); var lon = ee.Number(prev.get('lon',0)); var idx = ee.Number(prev.get('idx',0)); var bit = ee.Number(prev.get('bit',0)); var evenBit = ee.Number(prev.get('evenBit',1)); var geohash = ee.String(prev.get('geohash','')); var latMin = ee.Number(prev.get('latMin',-90)); var latMax = ee.Number(prev.get('latMax',90)); var lonMin = ee.Number(prev.get('lonMin',-180)); var lonMax = ee.Number(prev.get('lonMax',180)); // calculate substep bit step idx // bisect E-W longitude var lonMid = ee.Number(ee.Algorithms.If(evenBit.gt(0), lonMin.add(lonMax).
divide(2), 0) ); idx = ee.Number(ee.Algorithms.If(evenBit.gt(0).
and(lon.gte(lonMid)), idx.multiply(2).
add(1), idx) ); lonMin = ee.Number(ee.Algorithms.If(evenBit.gt(0).
and(lon.gte(lonMid)), lonMid, lonMin) ); idx = ee.Number(ee.Algorithms.If(evenBit.gt(0).
and(lon.lt(lonMid)), idx.multiply(2), idx) ); lonMax = ee.Number(ee.Algorithms.If(evenBit.gt(0).
and(lon.lt(lonMid)), lonMid, lonMax) ); // bisect N-S latitude var latMid= ee.Number(ee.Algorithms.If(evenBit.eq(0), latMin.add(latMax).
divide(2), 0) ); idx = ee.Number(ee.Algorithms.If(evenBit.eq(0).
and(lat.gte(latMid)), idx.multiply(2).
add(1), idx) ); latMin = ee.Number(ee.Algorithms.If(evenBit.eq(0).
and(lat.gte(latMid)), latMid, latMin) ); idx = ee.Number(ee.Algorithms.If(evenBit.eq(0).
and(lat.lt(latMid)), idx.multiply(2), idx) ); latMax = ee.Number(ee.Algorithms.If(evenBit.eq(0).
and(lat.lt(latMid)), latMid, latMax) ); // check position evenBit = evenBit.not(); bit = bit.add(1); geohash = ee.String(ee.Algorithms.If(bit.eq(5), geohash.cat(base32.slice(idx,ee.Number(idx).
add(1))), geohash)); idx = ee.Number(ee.Algorithms.If(bit.eq(5), ee.Number(0), idx)); bit = ee.Number(ee.Algorithms.If(bit.eq(5), ee.Number(0), bit)); // return state var curr = prev .
set('idx', idx ) .
set('bit', bit ) .
set('evenBit', evenBit) .
set('geohash', geohash) .
set('latMin', latMin ) .
set('latMax', latMax ) .
set('lonMin', lonMin ) .
set('lonMax', lonMax ); return ee.List([curr]); }; function geohash_encode(lat, lon, precision) { var init = ee.Dictionary({lat: lat, lon: lon}) var state = ee.List.sequence(1,precision*5).
iterate(geohash_accumulate, ee.List([init])) return ee.String(ee.Dictionary(ee.List(state).
get(-1)).
get('geohash'))
}
Давайте также проведем простую проверку полученных геохешей // PostGIS check from https://postgis.net/docs/ST_GeoHash.html
print (ee.Algorithms.If(geohash_encode(48, -126, 20).
compareTo('c0w3hf1s70w3hf1s70w3'),'Error','OK')); // Google BigQuery check from https://cloud.google.com/bigquery/docs/reference/standard-sql/geography_functions#st_geohash print (ee.Algorithms.If(geohash_encode(47.62, -122.35, 10).
compareTo('c22yzugqw7'),'Error','OK'));
Вместо использования вложенных циклов в функции geohash_encode() мы создаем функцию geohash_accumulate() и вызываем ее 5 раз для расчета каждого символа геохеша.
Учитывая, что GEE ограничивает количество вызовов функций, это не очень масштабируемое решение.
Как можно улучшить функцию? Самый простой способ — продублировать участок кода под комментарием «вычислить подшаг бит шага idx» пять раз и вызвать полученную функцию geohash_accumulate() только один раз для каждого вычисляемого символа (при желании можно также уменьшить количество переменных, передаваемых между звонки).
Или можно переписать код внутренней функции в функциональном стиле.
В общем, я не собирался этого делать, но пока писал текст статьи, мне стало интересно попробовать и вот что получилось: var geohash_accumulate = function(obj, list) {
// define common variables
var range4 = ee.List.sequence(0,4).
map(function(val){return ee.Number(val).
multiply(1/4)}); var range8 = ee.List.sequence(0,8).
map(function(val){return ee.Number(val).
multiply(1/8)}); // get previous state var prev = ee.Dictionary(ee.List(list).
get(-1)) var lat = ee.Number(prev.get('lat',0)) var lon = ee.Number(prev.get('lon',0)) var n = ee.Number(prev.get('n',0)).
add(1) var geohash = ee.String(prev.get('geohash','')) var latMin = ee.Number(prev.get('latMin',-90)) var latMax = ee.Number(prev.get('latMax',90)) var lonMin = ee.Number(prev.get('lonMin',-180)) var lonMax = ee.Number(prev.get('lonMax',180)) // calculate step n var base32 = ee.String(ee.Algorithms.If(n.mod(2).
eq(1), '028b139c46df57eghksujmtvnqwyprxz', '0145hjnp2367kmqr89destwxbcfguvyz')); var latRange = ee.List(ee.Number(ee.Algorithms.If(n.mod(2).
eq(1), range4, range8))); latRange = latRange.map(function(item){return ee.Number(item).
multiply(latMax.subtract(latMin)).
add(latMin)}); var lonRange = ee.List(ee.Number(ee.Algorithms.If(n.mod(2).
eq(1), range8, range4))); lonRange = lonRange.map(function(item){return ee.Number(item).
multiply(lonMax.subtract(lonMin)).
add(lonMin)}); var latIdx = latRange.indexOf(latRange.filter(ee.Filter.lte('item', lat)).
get(-1)); latIdx = ee.Number(ee.Algorithms.If(latIdx.gte(latRange.size().
add(-1)), latIdx.add(-1), latIdx)); var lonIdx = lonRange.indexOf(lonRange.filter(ee.Filter.lte('item', lon)).
get(-1)); lonIdx = ee.Number(ee.Algorithms.If(lonIdx.gte(lonRange.size().
add(-1)), lonIdx.add(-1), lonIdx)); var idx = lonIdx.multiply(latRange.size().
add(-1)).
add(latIdx); // reset bounds latMin = latRange.get(latIdx) latMax = latRange.get(latIdx.add(1)) lonMin = lonRange.get(lonIdx) lonMax= lonRange.get(lonIdx.add(1)) // define geohash symbol var geohash = geohash.cat(base32.slice(idx,ee.Number(idx).
add(1))); // return state var curr = prev .
set('n', n ) .
set('geohash', geohash ) .
set('latMin', latMin ) .
set('latMax', latMax ) .
set('lonMin', lonMin ) .
set('lonMax', lonMax ); return ee.List([curr]); }; function geohash_encode(lat, lon, precision) { var init = ee.Dictionary({lat: lat, lon: lon}) var state = ee.List.sequence(1,precision).
iterate(geohash_accumulate, ee.List([init])) return ee.String(ee.Dictionary(ee.List(state).
get(-1)).
get('geohash'))
}
Да, это та же самая функция, только в стиле Javascript и требует в пять раз меньше вызовов внутренних функций, что позволит нам на бесплатном аккаунте обрабатывать примерно в пять раз больше данных.
Заключение
Сегодня мы рассмотрели еще один кусочек головоломки из мира обработки пространственных данных.Как обычно, я рассказал о том, что есть в свободном доступе для каждого и легко тиражируется, поскольку весь код перед вами.
Разумеется, этот код также будет доступен на моем GitHub (см.
репозиторий).
Фрагменты ГИС в подкаталоге GEE).
Сразу не выложил, потому что сначала надо навести порядок — оказывается, многого я не накопил за два десятилетия работы с пространственными данными и программированием в целом.
Теги: #Популярная наука #программирование #открытый код #JavaScript #Визуализация данных #Геоинформационные сервисы #земной двигатель #суперкомпьютеры
-
Межамериканский Банк Развития (Idb)
19 Oct, 24 -
Расширенный Ip
19 Oct, 24 -
Сквозное Тестирование: Что, Зачем, Почему
19 Oct, 24
