Автоматизированная классификация снимков ДЗЗ Landsat

Вот моя history, все пока несколько хаотично. Первый линк прочитайте обязательно, во втором есть ссылки, где взять программы.
К третьему линку у меня есть поправки, и еще будет скрипт по векторизации классифицированных тайлов
и заливке результатов.

В любом случае рисовать с помощью этих картинок гораздо удобнее, чем по кривым композитам
на общеизвестных сайтах. Сможем ли мы извлечь полезное из обработки IRS таким же способом,
увидим позже.


# http://toroid.org/ams/landsat-and-grass
# http://www.custom-scenery.org/Satellite-Image.304.0.html
# http://www.custom-scenery.org/Building-Scener.331.0.html

wget ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p164/r020/p164r020_7x20010619.ETM-EarthSat-Orthorectified/*_nn*.tif.gz

grass

# in GUI: create new location using georeferenced data: p164r020_7p20010619_z41_nn80.tif

r.in.gdal input=p164r020_7p20010619_z41_nn80.tif output=20010619_B80
r.in.gdal input=p164r020_7t20010619_z41_nn10.tif output=20010619_B10
r.in.gdal input=p164r020_7t20010619_z41_nn20.tif output=20010619_B20
r.in.gdal input=p164r020_7t20010619_z41_nn30.tif output=20010619_B30
r.in.gdal input=p164r020_7t20010619_z41_nn40.tif output=20010619_B40
r.in.gdal input=p164r020_7t20010619_z41_nn50.tif output=20010619_B50
r.in.gdal input=p164r020_7t20010619_z41_nn70.tif output=20010619_B70

r.mapcalc "20010619_ndvi=1.0*(20010619_B40-20010619_B30)/(20010619_B30+20010619_B40)"
r.colors map=20010619_ndvi color=ndvi
qgis

i.fusion.brovey -l ms1=20010619_B20 ms2=20010619_B40 ms3=20010619_B50 pan=20010619_B80 outputprefix=20010619_brovey_2_4_5
i.landsat.rgb -p  strength=97 red=20010619_brovey_2_4_5.red green=20010619_brovey_2_4_5.green blue=20010619_brovey_2_4_5.blue

d.mon x0
d.rgb  red=20010619_brovey_2_4_5.red green=20010619_brovey_2_4_5.green blue=20010619_brovey_2_4_5.blue

r.composite  red=20010619_brovey_2_4_5.red green=20010619_brovey_2_4_5.green blue=20010619_brovey_2_4_5.blue output=rgb_20010619_brovey_2_4_5
qgis

i.group group=20010619 subgroup=all input=20010619_B10,20010619_B20,20010619_B30,20010619_B40,20010619_B50,20010619_B70,20010619_B80

date=20010619; for band in 10 20 30 40 50 70 80; do (r.mapcalc ${date}_B${band}_mask="if(${date}_B${band},1)" &) ; done

date=20010619 ; r.mapcalc "${date}_mask=(${date}_B10_mask*${date}_B20_mask*${date}_B30_mask*${date}_B40_mask*${date}_B50_mask*${date}_B70_mask*${date}_B80_mask)"

g.copy rast=20010619_mask,MASK

g.region -a rast=MASK zoom=MASK

i.cluster group=20010619 subgroup=all sigfile=20010619_sign-20.txt classes=20 convergence=98.0 separation=1.5 reportfile=report_20010619-20.txt
i.maxlik group=20010619 subgroup=all sigfile=20010619_sign-20.txt class=20010619_maxlik-20 reject=20010619_maxlik_rej-20
qgis

i.cluster group=20010619 subgroup=all sigfile=20010619_sign-10.txt classes=10 convergence=98.0 separation=1.5 reportfile=report_20010619-10.txt
i.maxlik group=20010619 subgroup=all sigfile=20010619_sign-10.txt class=20010619_maxlik-10 reject=20010619_maxlik_rej-10
qgis

Круто! Просто мечта - Вкалывают роботы а не человек! :slight_smile:
Кто же спорит с тем что в рукопашную мы всю територию замумукаемся рисовать, только вот, что бы такое вкурить “обычному человеку”, чтобы в этом разобраться и попробовать сделать самому? :slight_smile:
И, отсюда, вопрос. А надо ли пробовать? В том смысле, какие мощности могут потребоваться, что бы обсчитать всю территорию тайги и т.п.? Какая группа людей может с этим справится? Если для этого нужно больше 5-10 человек, то как этот обсчет должен координироваться? Ведь если я правильно понял при различных фильтрах, настройках (может термины не те, но смысл я думаю понятен) мы можем получить различный результат?? Значит надо использовать в одном районе одни настройки всем, кто его обсчитывает? Как их согласовать?
Понимаю, что может быть ответы на вопросы есть в данных линках. И я тороплюсь с вопросами. Но они в первую очередь приходят в голову.
Конечно попробую вкурить. Начал читать линки. Но думаю вопросы обязательно возникнут.
Может ветку на форуме организовать с faq-ом и ответами на вопросы по этому поводу?? Я так понял, что работа в этом направлении движется?
Но если по теме обсуждения. Посмотрел примеры на http://forum.openstreetmap.org/viewtopic.php?pid=73310#p73310. Хотя на примерах карты обработанных подобным образом ясно видно, что это не человек рисовал, однако тег “происхождение информации” в этом случае, надо ставить, как мне кажется, обязательно.
А с другой стороны, даже такая информация это лучше чем белое пятно.

Надо придумать простой способ (без поднятия WMS сервера и тому подобного) сделать результат классификации фоном в josm.

По приведенной history пока осилил только до подсчета NDVI - GRASS крайне суровая ГИС :slight_smile:
Появилось несколько вопросов:

Весь свежий Landsat ETM+ идет с полосами (неисправный сенсор). Мне кажется, что он почти не пригоден для наших задач классификации. Или с полосами как-то научились бороться? Из свежих снимков скорее всего придется использовать Lansat TM, основное отличие которого от ETM+ это отсутствие восьмого панхроматического канала. Что нужно поменять далее в скриптах при отсутствии панхроматического канала?

не понял что делаем в этом месте

Где бы скачать файлы сигнатур. 20010619_sign-20.txt и 20010619_sign-10.txt - это измененные сигнатуры из CORINE? Интересно, подойдут ли они для широт Московской области.

WMS для самой классификации не нужен, и qgis тоже.

Дядюшкой Сэмом разработана. Меня несколько напрягала излишняя объектная ориентированность
(особенно в импорте/перепроектировании), но для пакетного перепроцессирования GRASS это то что
надо. Если конечно написать доступный tutorial для “реальных”, типа ОСМ, а не абстрактных задач.

Не будет только Brovey fusion, и разрешение будет 30м/pix. Классификация от разрешения не зависит.
Но для нашей задачи это неважно.
Если даже нам дать мультиспектральный worldview 2, мы утонем в совершенно других проблемах:
тени и т.д.

Маски у каналов смещены относительно друг друга, поэтому нужно пересечение масок, чтобы
избежать “цветных” краев.

Не путайте unsupervised и supervised classification.
Эти сигнатуры генерирутся i.cluster
http://grass.osgeo.org/wiki/Image_classification
самый простой вариант, они ориентированы на максимальную разделяемость классов.
Я использовал самый сложный (но и самый качественный) вариант radiometric & geometric
supervised, в нем сигнатуры рассчитывает i.gensigset используя тренировочную карту.
Обратите внимание, что сигнатуры i.gensigset и i.cluster несовместимы.
Сигнатуры i.gensigset вообще-то тоже плохо переносимы с одного места на другое,
у меня даже не получилось прилично использовать август 2002 года для августа 2009 года
на той же row/path.
В качестве тренировочной карты я использовал CORINE landcover (CLC2000) для той
части сцены, которая находится в Финляндии, поэтому и классы у меня соответственные,
ели и сосны в одном классе. Но такой фокус проходит только в пограничных с ЕС районах.
Чтобы разнести ели и сосны по разным классам, мне придется создавать на основе текущей
классификации новую тренировочную карту, где они разделены.
То есть ходить и смотреть, где чисто еловый лес, а где чисто сосновый.
При этом нужен некий компромисс: чем больше я нарисую вручную, тем лучше статистика,
лучше спектральное распознавание, но и больше ручной работы.
А вам в Московской области надо 1) решить какие классы использовать 2) нарисовать их
на карте 3) использовать эту карту как тренировочную для i.gensigset
Вот как-то так, если кратко.
Результаты unsupervised classification можно использовать для составления этой карты.
Например, вода уже хороща видна на NDVI (низкие значения <0.1, правда нужно учитывать и
низкое альбедо, но вы в любом случае знаете примерно где находится вода, по тому же ГШ например)

source= там стоит. Никакие 1000 китайцев или индусов лучше точно не нарисуют.
Мне интересно, как вы себе представляете сам процесс получения лучшей информации.

Так в том то и дело, что я не представляю! :slight_smile:
Дело в том, мне кажется, что мы немного забыли о теме дискуссии.
Тема, насколько я понял, была использовать ли спутниковую информацию для отрисовки карты или лучше оставить белое пятно? Вот как раз по к этому и относилась моя фраза “А с другой стороны, даже такая информация это лучше чем белое пятно”. Признаю, что немного неуклюжая фраза. (не надо было использовать слово “даже”. Не писатель я, однако.) Но, если вам показалось, что я чем то принизил ценность подобной информации, то приношу извинения и ей и вам. :slight_smile: На самом деле, я и руками и ногами за. Т.к. поездив по Сибири, знаю, что по большей части этой территории ни один ОСМер ножками при жизни ближайших трех-четырех-пяти и т.д поколений не пройдет. (Если только не профукаем эти территории и не отдадим их китайцам.) А карта (с известной степенью достоверности) в любом случае нужна.
И я готов этому содействовать. Понять бы еще как. :slight_smile: (понять, надеюсь поможет FAQ, и описание последовательности действий на вики, буде таковые появятся)
По поводу ссылок.
Я смотрел кусок карты по этой ссылке http://www.openstreetmap.org/?lat=60.95675&lon=29.12733&zoom=15
Теги смотрел в Потлаче. Сегодня еще раз посмотрел. Потыкав по полигонам нашел только два одиноких quarry и wetland с тегом source.
на остальных только inner и outer. Может я что не так смотрел?
А по поводу канала. Ничего не могу сказать хорошо отрисован. Но тег source есть только у канала и еще у пары дорог highway=road. И может так и надо, но у меня не отрисовались леса и несколько озер вокруг него. Зато есть куча каких то безымянных точек на западе территории. Это баг какой то, или что?

Необходимо создавать статью в вики, я думаю. А тему надо прилепить. Что скажите, товарищи модераторы?

Меня идея частичной механизации очень заинтересовала.

Разделил, прикрепил.

2 usm78-gis
Тему правильно назвал?

Ага, osm2mp уже на ЛО валится.

Processing ways...        gpc malloc failure: edge table creation

Что там за мегаполигоны?

По площади или числу вертексов ?
mapnik не без труда, но уже поднялся с колен:
http://www.openstreetmap.org/?lat=61.1543&lon=29.433&zoom=11
Заливаю я, кстати, josm’ом. bulk_upload_sax.py не смог осилить эти файлы.
До нарезки на полосы ogr2osm откушал 22GB RAM при процессировании :O)

Как будто бы ни то ни другое…
Мультиполигоны он все нормально пережёвывает, а валится уже на обычных домиках.
То ли память в этом gpc утекает, то ли ещё чего…

Надо поднимать сервер, а нам клиентские приложения, что бы они скачивали выданные сервером куски, обрабатывали их, и отсылали обратно. Я так понимаю, вычислительной работы много в этом деле. ДЗЗ@home, в общем))
Что бы это не было уделом гиков. Потому что, как я понял, дело определения идёт успешно.

UPD: mapsёrfer тоже одолел тайлы, красиво!

“Автоматизированная” звучит вроде лучше. “автоматическая” как-то не предусматривает вмешательство
человека, хотя я еще тот лингвист :slight_smile:

Пилите Шура, пилите.

Можно попробовать заюзать BOINC для этого дела:

http://boinc.berkeley.edu/

Наверянка найдётся N людей, которые компьютер на ночь не выключают и готовы предоставить его в качестве вычислительного ресурса для полезных дел. На мой взгляд use case для которого BOINC создавался как раз наш случай, ибо задача квантуется и большого информационного обмена между отдельными квантами не требуется.

ИМХО, скоро ждите финов в гости - взвоют от такой детализации. :slight_smile:

Там проходит такой highway=unclassified параллельно границе. С железным забором 3.5м и
безконтактными датчыками, финны любят фотографии публиковать на тему luovutettu Karjala.
Даже тот забор, который первыми снесли венгры на границе с Австрией во времена перестройки,
был всего 2.8м.
Кстати, вот здесь есть 1:20000 WMS 80 летней давности, с которым интересно сравнивать классификацию:
http://www.karjalankartat.fi
А если когда будет mapnik style с поддержкой всех типов леса, то картинки будут еще разнообразнее :slight_smile:

Я думаю, что СЗФО в принципе реально обсчитать используя CORINE в Финляндии, Эстонии и Латвии как
тренировочную карту, расширяя ее за счет уже классифицированных и перекрывающихся сцен.
Но для этого нужна серьезная подготовительная работа и это само по себе большая работа.
Налогоплательщику бы это обошлось как минимум в несколько миллионов евро, судя по тому, как работа
по созданию CORINE проводилась в ЕС.

Я бы мог залить и Ломоносовский район, и далее до границы с Эстонией, но там местные рисовальщики
уже много чего вручную нарисовали.

Курил последнюю страницу, но чего-то не понял: вы массово импортируете данные? А откуда эти данные?