Установка последней версии SagaGIS в QGIS

В профильном сообществе попросили написать гайд по установке актуальной версии SagaGIS в QGIS. Исполняю обещание. В принципе, в составе QGIS SagaGIS уже включён, но во-первых, это временное явление (вроде как в новых версиях решили исключить SagaGIS из ядра), а во-вторых, в составе установочного пакета QGIS поставляется какая-то уж очень древняя версия Saga — кажется 7-я, тогда как актуальная версия на момент написания заметки — 9.2. 

По пунктам, кратенько. Потом постараюсь добавить подробностей, когда будет больше свободного времени. 

1. Скачиваем последнюю версию SagaGIS. 
2. Распаковываем скачанный архив в какую-нибудь папку. Я предпочитаю OSGeo4W64\apps\, но можно любую другую.
3. В QGIS отключаем плагин ядра SAGA GIS provider и скачиваем Processing Saga NextGen Provider. В Настройках последнего указываем путь к SagaGIS — куда мы его распаковали на предыдущем шаге. Закрываем QGIS. 
4. В папке с SagaGIS лежит архив saga4qgis.zip, внутри которого есть инструкция, что надо делать с содержимым этого архива. Не буду повторять её, но суть такова: содержимое архива надо распаковать (удалив предварительно уже существующее) в папку     C:\Users\...\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\processing_saga_nextgen 
   Только не делайте это «в лоб», сперва найдите где конкретно лежит файл SagaNameDecorator.py и папка \description\ — удалите их и скопируйте на их место содержимое архива.
5. В принципе всё готово, можно запускать QGIS и радоваться возможности работать с последней версией SagaGIS. Но есть еще одна маленькая хитрость, сильно облегчающая работу. Надо открыть файл SagaAlgorithm.py и заменить

elif isinstance(param, (QgsProcessingParameterString, QgsProcessingParameterField)):
                command += ' -{} "{}"'.format(param.name(), self.parameterAsString(parameters, param.name(), context))

на

elif isinstance(param, (QgsProcessingParameterString, QgsProcessingParameterField)):
                command += ' -{} "{:.10}"'.format(param.name(), self.parameterAsString(parameters, param.name(), context))

Для чего нужен последний пункт: дело в том, что инструменты SagaGIS, работающие с векторными исходниками, на входе принимают только shape файлы SHP. А те, в свою очередь, согласно стандарту не поддерживают имена полей длиннее 10 символов. Поэтому если вы попробуете подсунуть в SagaGIS на вход какой-нибудь Geopackage с именем поля 11 символов и больше, то в исходном варианте SagaAlgorithm.py вы получите ошибку. Мой же код автоматически обрезает имя поля по 10 символу и в результате формируется правильный shape-файл. 
К слову сказать, похожие улучшения можно вставлять и в другие Processing Provider плагины — они все сделаны более или менее похожим образом.

Как выудить в точки рельеф из SRTM или AW3D30

Иногда данных о рельефе территории исследования так мало, что приходится прибегать к помощи космических технологий. Точность у этой информации, прямо скажем, не ахти, особенно на  застроенных или лесистых территориях, но уж что есть —  дарёному коню, как говорится. Знаю, уже есть алгоритмы, которые с помощью нейросетей умеют удалять с этих массивов дома и небольшие лесочки, но пока не встречал их в открытом доступе, пригодном к тому же для использования неквалифицированному пользователю.

Конечной целью предполагается получение CSV файла вида X, Y, Z.  Потребуется следующий инструментарий:

1. QGIS
2. Плагины к нему:
• Какой-нибудь плагин для подгрузки подложек типа QuickMapServices. Ну или загрузить их вручную через XYZ Connections.
• SRTM-Downloader
• Point sampling tool
3. Аккаунты как минимум на одном из сайтов с рельефом со спутников:
• https://search.earthdata.nasa.gov/ - регистрация по ссылке https://urs.earthdata.nasa.gov//users/new
  
• https://www.eorc.jaxa.jp/ALOS/a/en/dataset/aw3d30/aw3d30_e.htm - регистрация по ссылке https://www.eorc.jaxa.jp/ALOS/a/en/dataset/aw3d30/registration.htm

Сразу предупрежу, что выудить данные с японского сайта — тот еще квест, но в детали вдаваться не буду. В сети есть инструкции. С NASA возни гораздо меньше, особенно если вам повезло и SRTM-Downloader работает как надо.

Итак, алгоритм следующий:

1. В QGIS любым известным способом открываем картографическую подложку. Это не обязательно, но как правило сильно облегчает работу.
2. Создаем временный слой типа Polygon. Рисуем в нем полигон по контуру территории, для которой нам нужен рельеф. Зуммируемся так, чтоб нарисованный полигон занимал большую часть экрана (Zoom to Layer).
3. С помощью строенного в QGIS инструмента Vector/Random points in polygons генерируем достаточное количество точек в пределах нарисованного в предыдущем пункте полигона. При этом создается новый временный слой с точками.
4. С помощью SRTM-Downloader скачиваем растры с рельефом для области, показанной на экране. Для чего жмём последовательно Set Canvas Extend и Download. Для данных ALOS придется скачивать данные через сайт и подгружать в QGIS вручную.
5. Запускаем Point sampling tool, указываем слой с точками из пункта 3 в качестве Layer containing sampling points, а в качестве Layers with fields/bands to get values from — растры (можно указать сразу несколько) из пункта 4. Тут важно отметить, что точки и растры должны быть в одной системе координат (WGS 84 -  EPSG 4326 в случае если мы работаем с SRTM). 
6. На выходе получаем векторный слой с точками с атрибутами в виде значений рельефа. Если область интереса попадает сразу на несколько растров, то столбцов с атрибутами будет несколько. Надеюсь, не надо рассказывать, как из нескольких столбцов получить один — это можно сделать как через Field Calculator в самом QGIS, так и в Excel или любом другом табличном редакторе.
7. Пересохраняем полученный точечный слой в виде CSV, не забывая выбрать нужную нам систему координат (если не выбирать, то по умолчанию  координаты будут представлены в виде градусов) и указать, что GEOMETRY сохраняется в виде AS_XY.

Отчеты в редактируемом формате

Довольно долгое время я передвал заказчикам отчеты исключительно в формате “doc”. Мне казалось, что это удобно — в случае чего там самостоятельно смогут добавить титульный лист или подправить что-нибудь в тексте по-мелочи. Все было отлично, пока сразу две организации независимо друг от друга за весьма короткий промежуток времени не подставили меня перед экспертами, переиначив втихаря мои выводы на фактичекски противоположные. Покумекав, я решил, что теперь только формат “pdf” — мой единственный друг. Пришлось поступиться некоторыми удобствами — оформлять титульные листы приходится самостоятельно и с мелкими правками воевать тоже.
Потом я и вовсе отказался от MSO в пользу LibreOffice и если заказчик вдруг просит текст отчета в т.н. «редактируемом формате» (хотя и “pdf” вполне себе редактируемый), то отдаю ему “docx” без формления (рамочек, штампа — которые из LibreOffice в этот формат нормально все-равно не сохраняются).
Для копирования в свои отчеты клиентам этого будет достаточно, а вот для переделки выводов придется покорячиться. В копилку к известной мудрости не передавать заказчикам предварительные результаты — не отдавать им отчеты в «легко фальсифицируемом формате».

Как бороться с осушенными ячейками

Я уже писал когда-то на эту тему. Пришло время немного обновить советы.
Модельные ячейки осушаются в том случае, когда уровень воды в безнапорной ячейке на каком-то из этапов расчета оказывается ниже её подошвы. В ряде случаев, это абсолютно нормальное явление, хотя по-возможности лучше этого избегать, но в большинстве ситуаций — такое поведение модели является результатом ошибки. Если площадные размеры модельной ячейки сильно превышают ее мощность, а градиент фильтрационного потока очень большой — с большой степени вероятности ячейка осушится. Когда подошва слоя активно «ныряет» вниз и «взлетает» вверх — тоже жди неприятностей.
Несколько способов борьбы с этим явлением:

1. Смириться. В ряде случаев, как я уже говорил, осушение — это вполне ожидаемая модельная картина, но сразу предупрежу: делать потом на такой модели прогнозные расчеты будет очень сложно, если вообще возможно. 
2. Дополнительно дискретизировать модельную сетку на проблемных участках. Способ очевидный, однако он не всегда помогает, а если и помогает, то увеличивает время счета, что критично при многовариантных расчетах или при использовании модулей автоматического подбора параметров типа PEST. 
3. Использовать опцию повторного обводнения (rewetting). Такая опция стала доступна в MODFLOW-2000/2005, однако будьте готовы к проблемам со сходимостью. Кроме того, скорость расчета после включения этого режима падает чуть ли не на порядок. 
4. Воспользоваться MODFLOW-NWT. Этот алгоритм был специально разработан для борьбы с нежелательным осушением ячеек. Работает в целом хорошо, стабильно и сравнительно быстро, но иногда излишне «сопротивляется» этому самому осушению, как-бы «размазывая» УГВ тонким слоем по подошве слоя.

Расчет-недорасчет

Заметка уже не актуальна.
Бывает так, что нормативная база входит в противоречие со здравым смыслом. Скажем больше, в отдельных областях хозяйственной деятельности такое положение вещей является скорее нормой, чем исключением. В нашем случае такое тоже случается. К примеру, по московским нормативам, если подземная часть проектируемого здания будет заглублена более чем на 10 метров от поверхности земли, то для такого здания надо делать расчет влияния на грунтовые воды. Даже если вода залегает на 25 метрах и даже в страшном сне не поднимается до отметок, сопоставимых с дном котлована. А «гидропрогноз» делать надо, никуда не денешься, иначе экспертизу не пройдешь. Что делать?
Что делать, что делать. Выкручиваться. Естественный уровень грунтовых вод в течение года «гуляет» вверх-вниз с амплитудой примерно в 1-2 метра. Соответственно, если от дна котлована до УГВ меньше метра, то можно посчитать прогноз для повышенных (паводковых) уровней подземных вод, задав повышенную инфильтрацию, к примеру. Котлован естественно затопит, с чем мы доблестно будем бороться методами открытого водоотлива и даже будем знать примерную величину водопритока. Ерунда, конечно, но формально требования законодательства мы выполнили.
Однако бывает так, что никакими сезонными колебаниями не заставишь подземные воды затопить котлован. В этом случае я даю прогнозную картинку с «паводковыми» уровнями, показываю на ней, что даже в случае ужасного паводка котлован останется сухим. Далее для наукообразия делаю простой расчет: Δω=ω∙S, где: ω — величина инфильтрации, а S — площадь здания, ну а Δω, соответственно, — величина дефицита инфильтрации, возникающего в результате строительства. Халтура, но опять же, формально влияние здания на подземные воды учтено.

Перевод отметок рельефа из Autocad в табличный вид

При подготовке данных для моделирования весьма часто возникает необходимость в конвертации отметок рельефа из геоподосновы, сохраненной в формате Autocad DWG, в табличный вид (типа X, Y, Z) для того, чтоб скормить эти данные какому-нибудь интерполятору (Surfer или встроенный в PmWin “Field Interpolator”).
Я почти уверен точно знаю, что эта задача может быть легко и быстро решена с помощью самого Autocad — достаточно запустить соответствующую программу на LISP-е и радоваться жизни. К сожалению, я LISP-а не знаю, да и вообще не являюсь большим специалистом в автокаде.
В современных версиях автокада эта проблема решается еще проще: через инструмент, расположенный в пункте меню Tools\Data Extraction. Инструмент довольно мощный, но в нестандартных случаях возможно придется повозиться.
Но проблему как-то надо решать. Я предлагаю использовать для этого MapInfo (согласен, для кого-то это выглядит сменой шила на мыло). Далее по пунктам:

1. Конвертируем слой с отметками (желательно, чтоб слой содержал только отметки в текстовом виде, без самих точек) из формата Autocad DWG  в формат MapInfo TAB (с помощью встроенного в MapInfo мини-приложения Universal Translator).
2. Подчищаем полученную таблицу от нетекстовых элементов: это можно сделать несколькими способами, наиболее удобный и быстрый — мини-приложение MapCad, но можно и с помощью Query Select и функции ObjectInfo(obj, 1), но там придется сначала создать дополнительный столбец в таблицу, занести него результат выполнения функции ObjectInfo(obj, 1), а уж потом делать Query Select по этому столбцу, выбирая значения, отличные от 10 (а 10 — это как раз текстовые).
3. Еще разок запускаем Update Column (создайте новый столбец с типом float или смените тип существующего столбца) с той же функцией, но с другими параметрами: ObjectInfo(obj, 3). Если все сделано правильно, то в вашей таблице появится столбец типа float с отметками рельефа.
4. Запускаем мини-приложение Coordinate Extractor: в таблице теперь будут столбцы с координатами центра текстовой подписи отметки рельефа. Вот тут важно отметить явный недостаток рассматриваемого метода: наши точки будут немного смещены относительно реальных отметок — ровно на столько, на сколько отличаются координаты середины метки от координат точки замера. Если вы страдаете перфекционизмом, то эту проблему можно решить с помощью простейших математических операций со свежеполученными координатами.
5. Запускаем Table/Create Points, если хотим заменить тектовые метки на точки (а уж сами метки пусть MapInfo своими силами рисует, благо соответствующий столбец в таблице уже есть).
6. Экспортируем полученную таблицу в нужный текстовый формат (txt или csv).

Метода только выглядит громоздко, на самом же деле, у меня уходит на все эти действия не больше минуты — главное не сбиваться и соблюдать порядок действий.

Моделирование в понижениях

Коллега в телефонной беседе напомнил про один очень полезный и неочевидный трюк. Как известно, нестационарные модели очень чувствительны к начальному распределению напоров. Получить это распределение можно двумя путями: интерполяцией фактически наблюденных напоров в наблюдательных точках (скважинах); либо с помощью откалиброванной стационарной модели. Очевидны недостатки обоих подходов: для первого не всегда хватает исходных данных, а второй чреват ошибками (все упирается в качество решения обратной задачи).
Однако, существует массивный пласт задач, для которых совершенно не обязательно знание о распределении напоров. К примеру: определение водопритока к котловану. Нам совершенно неважно какой там напор, в качестве результата нас устроит достигнутое понижение и расход.
В этом случае допустим такой приём: заменяем в нашем расчете термин «напор» на термин «понижение». Таким образом, «начальный напор» становится «начальным понижением», которое, очевидно, равно нулю. Все модельные слои должны быть напорными (Confined), даже если они не являются таковыми. Важно помнить: изолинии и массивы результатов, получаемые при таком расчете, будут показывать не напоры, а понижения. Расходы будут такими же, как и при расчете в термине «напор».
Совершенно ясно, что и этот подход не лишен недостатков и ограничений. Так, без дополнительных ухищрений невозможно учесть эффект осушения пластов, а величины расходов могут оказаться существенно завышенными для безнапорных горизонтов, как это всегда бывает, если их задать напорными.

Задание геометрии слоев

Не знаю как для вас, но для меня самая, если можно так выразиться, нелюбимая часть моделирования — это процесс задания рельефа кровель и подошв водоносных и водоупорных слоев. Основные сложности возникают с выклиниванием (MODFLOW как известно не поддерживает разрывов слоев по простиранию и слоев нулевой мощности) и экстраполяцией (ну, т.е. когда имеющиеся скважины пробурены не на всю площадь моделирования — приходится иногда фантазировать).
Однако, с опытом набирается некий запас хитростей и know-how по решению таких задач. Выклинивающийся слой можно сделать околонулевой мощности — слой как-бы останется, но на решение уже влиять не будет. Другой вариант: в выклинившемся слое задать параметры аналогичные подстилающему или надстилающему слою — это сложнее, но с методологической точки зрения более грамотно. Важно помнить, что логичное казалось бы действие: сделать слой неактивным за границей выклинивания, приведет к образованию глухой непроницаемой границы для вертикального перетока. Само собой, это актуально для моделирования многослойной толщи и когда выклинивающийся слой находится посередине. Окажись он сверху или снизу — этот метод вполне применим. С экстраполяцией сложнее: приходится внимательно изучать геологические карты района, смотреть разрезы, выискивать архивные скважины. Творческая задача, одним словом.
У тех, кто работает с горно-складчатыми областями свои погремушки: bedrock folding, when not to use interpolation.

Импорт WMF в MapInfo

Возникла у меня не совсем типичная задачка: прислали мне векторную карту, сохраненную в формате WMF (Windows Media File). Мне для работы желательно загнать эту карту в MapInfo. Решение «в лоб»: привязать wmf, как растровую подложку. Не сработало. MapInfo в принципе это умеет делать, но с большими файлами (а у меня он весит почти 5 мегабайт) не дружит — показывает от силы половину карты, а то и того меньше. К тому же, как-то некомильфо получается — векторный исходник подкладывать растром в векторную же программу.
Пришлось идти по длинному пути: сначала открываем наш WMF в AutoCAD-е с помощью команды _wmfin, не меняя масштаба по осям, потом командой _explode разбиваем блок (на всякий случай, т.к. у MapInfo случаются затыки с импортом блоков), сохраняем файл в виде DXF. Транслируем этот DXF в MapInfo. С помощью MapBasic-утилиты Register Vector перемещаем и масштабируем карту. Увы, если у вас был на карте текст, то он скорее всего потеряется. Особенно, если он написан на кириллице. Для решения этой проблемы, по всей видимости, придется перемещать и масштабировать карту сразу в AutoCAD-е, еще до импорта ее в MapInfo.

И еще раз об осушенных блоках модели

Возвращаясь к напечатанному. Довольно типичная ситуация: «стена в грунте», внутри стены задан скважинный дренаж. Откачивающие скважины заданы как граничное условие второго рода. Запускаем расчет и получаем вот такую картину:

Такое впечатление, что из шести скважин работают только верхние две. Самое забавное, что так оно и есть, но в чем же дело?
Вот как вся эта ситуация это выглядит при итеративном расчете:

1. Работает скважина, снижается уровень воды.
2. Вдруг уровень стал ниже подошвы слоя — отключили ячейку со скважиной.
3. О, скважина исчезла (ячейка то с ней отключена) — уровень воды начинается повышаться.

А проблема в том, что отключенная из-за осушения ячейка не включится, а следовательно и скважина не заработает.
Как бороться с такой напастью. Да очень просто, почти все рецепты я уже описал в своей предыдущей записи на эту тему, к которым можно добавить еще один, а именно: задать скважины не расходом, а напором (понижением) — кстати, вовсе не обязательно использовать для этого первый род, можно попытаться учесть сопротивление прискважинной зоны и задать скважины третьим родом (пакетами General Head или даже лучше Drain). А расход скважин можно получить уже из Water Budget.

Осушающиеся блоки модели

MODFLOW очень «не любит» маломощные безнапорные горизонты. Когда уровень воды в ячейке модели становится очень близок к отметке ее подошвы, MODFLOW ее делает неактивной и больше в расчете не использует. Даже если вокруг потоп и все окружающие эту ячейку блоки затоплены «под завязку», MODFLOW ее не включит в расчет и не обводнит. Следствием такой особенности являются проблемы с точностью и сходимостью расчета. Если отключенных ячеек оказывается относительно много, то модель почти гарантированно не сойдется при умолчальных настройках расчетных модулей (Solvers).
Для борьбы с этой напастью придумано масса способов:

• Можно поиграться с настройками солверов. Особенно помогает уменьшить Damping Parameter в настройках модуля PCG2.
• Иногда помогает просто сменить солвер. Неплохо сходится солвер GMG, но он поддерживает только MODFLOW-2000 и не все препроцессоры с ним совместимы.
• Полезно воспользоваться пакетом «обводнение» (Wetting Capability), который специально создан для борьбы с этим эффектом. К сожалению, само по себе использование этого пакета может привести к проблемам со сходимостью. За что боролись, как говорится. Хотя иногда он очень помогает, особенно при нестационарных расчетах процессов обводнения изначально сухих горизонтов.
• Самый радикальный способ — отказаться от расчета в безнапорной постановке и считать все слои напорными. В этом случае MODFLOW не будет отключать никакие сухие ячейки, в силу того, что проводимость в них не зависит от уровня воды. Такой подход приводит к заведомо большим расходам потока, но зато позволяет быстро получить хоть какой-нибудь результат. К тому же, при расчете тех же дренажей, небольшое увеличение расхода никакой беды не представляет, создавая некий запас прочности (главное помнить, что таких «запасов» часто оказывается не один и не два и в сумме они нарисуют такой  «запасище», что проектанты на стенку полезут).

Профильные модели

Все-таки есть небольшая польза от тупого просиживания штанов в офисе. Придумал способ, как избежать длительного и мучительного итеративного рисования свободной поверхности на профильных моделях в MODFLOW.
Нужно в качестве подошвы ячеек задать координату y их нижней границы. Следовательно, когда уровень в блоке становится ниже подошвы, то ячейка сама отключается, что нам и нужно. Не стоит забывать о сохранении единичной мощности пласта.
Остаются проблемы с инфильтрацией: не совсем понятно, куда ее задать (на всю «площадь» некорректно, да и приводит к тому же к необходимости итеративного же ее уточнения), а при неизвестной верхней границе — невозможно задать на ней инфильтрацию.
Ну и обязательно надо включить модуль wetting capability, что традиционно приводит к проблемам со сходимостью.