Расчет площадей поверхностей с учетом рельефа средствами GRASS
Обсудить в форуме Комментариев — 2Редактировать в вики
по адресу http://gis-lab.info/qa/3darea-grass.html
Время от времени возникает задача рассчитать площадь поверхности с учетом ее наклона, т.е. того, что эта поверхность находится в трехмерном пространстве. Эта статья рассматривает пример того, как эту задачу можно решить в GRASS GIS.
Для примера возьмем набор geosample (http://gis-lab.info/qa/geosample.html). Данный набор содержит растр высот (слой relief) и векторную карту охраняемых природных территорий (слой oopt). Требуется посчитать площади каждой из охраняемых территорий как с учетом рельефа, так и без этого.
[править] Перепроецирование
Поскольку набор geosample содержит данные в географической системе координат (широта/долгота), а измерения требуется произвести в метрической системе координат, то их предварительно нужно перепроецировать, например в равноплощадную проекцию Альберса. Для этого создадим еще один набор geosample_sib_aea в новой проекции со следующими параметрами:
north: 6184472.19979012 south: 5008416.36981374 west: 16373043.44324575 east: 17599228.40562719 nsres: 659.6 ewres: 659.6
Заходим в этот набор и перепроецируем растровые данные (здесь и ниже галочка ">" означает приглашение командной строки GRASS, после нее вводится сама команда, в данном случае команда "r.proj loc=geosample in=relief out=relief"):
> r.proj loc=geosample in=relief out=relief
Аналогично перепроецируем векторные данные:
> v.proj loc=geosample in=oopt out=oopt
Посмотрим поближе на исходные данные (слой relief):
> r.info relief +----------------------------------------------------------------------------+ | Layer: relief Date: Wed Aug 17 18:08:10 2011 | | Mapset: PERMANENT Login of Creator: dima | | Location: geosample_sib_aea | | DataBase: /home/dima/laboro/GRASSDATA | | Title: ( relief ) | | Timestamp: none | |----------------------------------------------------------------------------| | | | Type of Map: raster Number of Categories: 4197 | | Data Type: CELL | | Rows: 1782 | | Columns: 1859 | | Total Cells: 3312738 | | Projection: Albers Equal Area | | N: 6184472.19979012 S: 5009075.96366118 Res: 659.59384743 | | E: 17599228.40562719 W: 16373043.44324575 Res: 659.5938474 | | Range of data: min = 20 max = 4197 | | | | Data Description: | | generated by r.proj | | | | Comments: | | r.proj input="relief" location="geosample" output="relief" method="n\ | | earest" | | | +----------------------------------------------------------------------------+
Аналогично по слою oopt:
> v.info oopt +----------------------------------------------------------------------------+ | Layer: oopt | | Mapset: PERMANENT | | Location: geosample_sib_aea | | Database: /home/dima/laboro/GRASSDATA | | Title: | | Map scale: 1:1 | | Map format: native | | Name of creator: dima | | Organization: | | Source date: Wed Aug 17 16:36:58 2011 | |----------------------------------------------------------------------------| | Type of Map: vector (level: 2) | | | | Number of points: 0 Number of areas: 8 | | Number of lines: 0 Number of islands: 8 | | Number of boundaries: 8 Number of faces: 0 | | Number of centroids: 8 Number of kernels: 0 | | | | Map is 3D: No | | Number of dblinks: 1 | | | | Projection: Albers Equal Area | | N: 5929488.79257984 S: 5114805.91166458 | | E: 17470617.17073559 W: 16788858.38019306 | | | | Digitization threshold: 0 | | Comments: | | | +----------------------------------------------------------------------------+
[править] Расчет площадей
GRASS содержит модуль r.surf.area (http://grass.osgeo.org/gdp/html_grass64/r.surf.area.html), который производит требуемые расчеты. Рассчитаем общую площадь всей территории охвата слоя relief с использованием данного модуля:
> r.surf.area relief Null value area ignored in calculation 5.439092e+11 Plan area used in calculation: 1.440478e+12 Surface Area Calculation(low, high, avg): 9.009841e+11 9.024818e+11 9.017329e+11 Current Region plan area: 1.442062e+12 Estimated Region Surface Area: 9.027246e+11 Done.
Остановимся подробнее на результатах расчета:
- Null value area ignored in calculation: площадь ячеек растра, которые не участвовали в расчетах (были замаскированы).
- Plan area used in calculation: площадь всего растра как площадь проекции растра на горизонтальную плоскость.
- Surface Area Calculation(low, high, avg): площадь растра с учетом рельефа. Рассчитывается путем построения треугольников, вершины которых находятся в центре ячеек растра. Поскольку для соседних ячеек можно построить несколько вариантов триангуляции (см. пример на рисунке), то и результаты расчета площадей могут различаться. Поэтому алгоритм сохраняет минимальное, максимальное и среднее ((минимальная + максимальная)/2) значение расчетной площади для каждой пары треугольников.
- Current Region plan area: площадь текущей области (без учета рельефа).
- Estimated Region Surface Area: итоговая расчетная площадь территории.
Таким образом, получаем, что величина 9.027246e+11 -- искомая площадь территории с учетом рельефа, 1.440478e+12 -- общая площадь без учета рельефа (если необходимо произвести сравнение площадей с учетом и без учета рельефа, то из общей площади 1.440478e+12 сначала нужно вычесть число 5.439092e+11 -- площадь пустых ячеек, которые в расчетах не участвовали).
Нужно заметить, что результаты расчетов модуля зависят от текущего разрешения (что вполне естественно, т.к. результат и должен зависеть от точности исходных данных).
Посмотрим на текущее разрешение:
> g.region -p projection: 99 (Albers Equal Area) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: krassovsky north: 6184472.19979012 south: 5008416.36981374 west: 16373043.44324575 east: 17599228.40562719 nsres: 659.59384743 ewres: 659.59384743 rows: 1783 cols: 1859 cells: 3314597
Огрубим разрешение:
> g.region res=6596 -p projection: 99 (Albers Equal Area) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: krassovsky north: 6184472.19979012 south: 5008416.36981374 west: 16373043.44324575 east: 17599228.40562719 nsres: 6607.05522459 ewres: 6592.39227087 rows: 178 cols: 186 cells: 33108
И рассчитаем площади с новым разрешением:
> r.surf.area relief Null value area ignored in calculation 5.440617e+11 Plan area used in calculation: 1.426251e+12 Surface Area Calculation(low, high, avg): 8.825813e+11 8.827467e+11 8.826640e+11 Current Region plan area: 1.442062e+12 Estimated Region Surface Area: 8.924489e+11
Как видим, результаты разнятся, и это нужно учитывать при работе.
Вернем разрешение к исходному разрешению растра:
> g.region rast=relief -p projection: 99 (Albers Equal Area) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: krassovsky north: 6184472.19979012 south: 5009075.96366118 west: 16373043.44324575 east: 17599228.40562719 nsres: 659.59384743 ewres: 659.59384743 rows: 1782 cols: 1859 cells: 3312738
[править] Расчет площадей по отдельным участкам
Однако, наша задача несколько иная: нужно расчитать площади по каждому из заповедников в отдельности. Для того, чтобы это сделать, необходимо следующее:
- Извлечь из карты oopt интересующий нас полигон.
- Создать маску из этого полигона, чтобы команда "r.surf.area relief" производила расчеты только в замаскированной области
- Вызов команды r.surf.area relief и сохранение результата
Реализация: Посмотрим содержимое таблицы атрибутов карты oopt:
> v.db.select oopt cat|NAME_PRT_R|NAME_R|TYPE 1|Катунский|Катунский|Заповедник 2|Кузнецкий Алатау|Кузнецкий Алатау|Заповедник 3|Участок "Ханхаринский"|Тигирекский|Заповедник 4|Участок "Белорецкий"|Тигирекский|Заповедник 5|Участок "Тигирекский"|Тигирекский|Заповедник 6|Алтайский|Алтайский|Заповедник 7|Шорский|Шорский|Национальный 8|Кирзинский|Кирзинский|Заказник
Создадим из векторной карты oopt растровую:
> v.to.rast in=oopt out=oopt use=attr col=cat label=NAME_PRT_R
Посмотрим, что из этого получилось, заодно рассчитаем статистику по площадям:
> r.report oopt un=k 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: geosample_sib_aea Wed Aug 17 18:38:50 2011| |-----------------------------------------------------------------------------| | north: 6184472.19979012 east: 17599228.40562719 | |REGION south: 5009075.96366118 west: 16373043.44324575 | | res: 659.59384743 res: 659.59384743 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: Labels (oopt in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | square | |#|description |kilometers| |-----------------------------------------------------------------------------| |1|Катунский . . . . . . . . . . . . . . . . . . . . . . . . . . . | 1501| |2|Кузнецкий Алатау. . . . . . . . . . . . . . . . . . . . . . . . | 3993| |3|Участок "Ханхаринский". . . . . . . . . . . . . . . . . . . . . | 4| |4|Участок "Белорецкий". . . . . . . . . . . . . . . . . . . . . . | 388| |5|Участок "Тигирекский". . . . . . . . . . . . . . . . . . . . . .| 17| |6|Алтайский . . . . . . . . . . . . . . . . . . . . . . . . . . . | 8770| |7|Шорский . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 4201| |8|Кирзинский . . . . . . . . . . . . . . . . . . . . . . . . . . .| 1196| |*|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 1,421,182| |-----------------------------------------------------------------------------| |TOTAL | 1,441,253| +-----------------------------------------------------------------------------+
Создадим маску, соответствующую участу "Катунский" (категория 1):
> r.mask oopt maskcat=1
Расчитаем площадь для этого участка:
> r.surf.area relief Null value area ignored in calculation 1.439753e+12 Plan area used in calculation: 1.439670e+12 Surface Area Calculation(low, high, avg): 1.447004e+09 1.467648e+09 1.457326e+09 Current Region plan area: 1.441253e+12 Estimated Region Surface Area: 1.458929e+09 Done.
Таким же образом можно получить информацию и по остальным участкам.
Эту процедуру легко автоматизировать, сохранив данные расчетов в файл, например, area.txt:
> for id in `seq 8` do r.mask -r r.mask oopt maskcat=$id r.category oopt cat=$id >> area.txt r.surf.area relief >> area.txt done
Начальные строки файла area.txt представлены ниже:
1 Катунский Null value area ignored in calculation 1.439753e+12 Plan area used in calculation: 1.439670e+12 Surface Area Calculation(low, high, avg): 1.447004e+09 1.467648e+09 1.457326e+09 Current Region plan area: 1.441253e+12 Estimated Region Surface Area: 1.458929e+09 Done. 2 Кузнецкий Алатау Null value area ignored in calculation 1.437260e+12 Plan area used in calculation: 1.439670e+12 Surface Area Calculation(low, high, avg): 3.870987e+09 3.888315e+09 3.879651e+09 Current Region plan area: 1.441253e+12 Estimated Region Surface Area: 3.883918e+09 Done.
Обсудить в форуме Комментариев — 2Редактировать в вики
Последнее обновление: 2014-05-14 21:52
Дата создания: 19.01.2012
Автор(ы): Дмитрий Колесов
© GIS-Lab и авторы, 2002-2021. При использовании материалов сайта, ссылка на GIS-Lab и авторов обязательна. Содержание материалов - ответственность авторов. (подробнее).