Расчет площадей поверхностей с учетом рельефа средствами GRASS: различия между версиями
Нет описания правки |
|||
| (не показана 1 промежуточная версия 1 участника) | |||
| Строка 115: | Строка 115: | ||
* Estimated Region Surface Area: итоговая расчетная площадь территории. | * Estimated Region Surface Area: итоговая расчетная площадь территории. | ||
[[Файл: | [[Файл:Ambiguous_triangilation_1.png|thumb|На рисунке показаны два различных варианта построения треугольников при расчете площади поверхности (синими точками обозначены узлы решетки растра -- 9 ячеек)]] | ||
[[Файл:Ambiguous triangilation3D.png|thumb|На рисунке показаны два различных варианта построения треугольников при расчете площади поверхности (ячейки растра представлены вертикальными столбиками -- 4 ячейки)]] | [[Файл:Ambiguous triangilation3D.png|thumb|На рисунке показаны два различных варианта построения треугольников при расчете площади поверхности (ячейки растра представлены вертикальными столбиками -- 4 ячейки)]] | ||
Текущая версия от 06:18, 19 марта 2012
по адресу 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.