Просмотр кода страницы Расчет площадей поверхностей с учетом рельефа средствами GRASS

{{Статья|Черновик}}
Время от времени возникает задача рассчитать площадь поверхности с учетом ее наклона, т.е. того, что эта поверхность погружена в трехмерное пространство. Данная статья рассматривает пример того, как эту задачу можно реализовать в GRASS GIS.


Для определенности возьмем набор geosample (http://gis-lab.info/qa/geosample.html). Данный набор содержит растр высот relief и векторную карту охраняемых природных территорий oopt.
Требуется посчитать площади каждой из охраняемых территорий как с учетом рельефа, так и без этого.

== Перепроецирование ==
Поскольку набор geosample содержит данные в системе широта/долгота, то их предварительно нужно перепроецировать, например в равноплощадную проекцию Альберса. Для этого создадим еще один набор geosample_sib_aea в требуемой проекции со следующими параметрами:

<syntaxhighlight lang="bash">
north:      6184472.19979012
south:      5008416.36981374
west:       16373043.44324575
east:       17599228.40562719
nsres:      659.6
ewres:      659.6
</syntaxhighlight>

Заходим в этот набор и перепроецируем растровые данные (здесь и ниже галочка ">" означает приглашение командной строки GRASS, после нее вводится сама команда, в данном случае команда "r.proj loc=geosample in=relief out=relief"):
<syntaxhighlight lang="bash">
> r.proj loc=geosample in=relief out=relief
</syntaxhighlight>

Аналогично перепроецируем векторные данные:
<syntaxhighlight lang="bash">
> v.proj loc=geosample in=oopt out=oopt
</syntaxhighlight>

Посмотрим поближе на исходные данные:
<syntaxhighlight lang="bash">
> 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"                                                                 |
 |                                                                            |
 +----------------------------------------------------------------------------+
</syntaxhighlight>

Аналогично по oopt:
<syntaxhighlight lang="bash">
> 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:                                                                |
 |                                                                            |
 +----------------------------------------------------------------------------+
</syntaxhighlight>

== Расчет площадей ==
GRASS содержит модуль r.surf.area (http://grass.osgeo.org/gdp/html_grass64/r.surf.area.html), который производит требуемые расчеты. 
Нужно заметить, что результаты расчетов модуля зависят от текущего разрешения (что вполне естественно, т.к. результат и должен зависеть от точности исходных данных).
Посмотрим на текущее разрешение:
<syntaxhighlight lang="bash">
> 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
</syntaxhighlight>

Теперь расчитаем площади:
<syntaxhighlight lang="bash">
> 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.
</syntaxhighlight>
Величина 9.027246e+11 -- искомая площадь территории с учетом рельефа, 1.440478e+12 -- общая площадь без учета рельефа (если необходимо произвести сравнение площадей с учетом и без учета рельефа, то из общей площади 1.440478e+12 сначала нужно вычесть число 5.439092e+11 -- площадь пустых ячеек, которые в расчетах не участвовали).

Огрубим разрешение:
<syntaxhighlight lang="bash">
> 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
</syntaxhighlight>
И рассчитаем площади с новым разрешением:
<syntaxhighlight lang="bash">
> 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
</syntaxhighlight>
Как видим, результаты разнятся, и это нужно учитывать при работе.

Вернем разрешение к исходному разрешению растра:
<syntaxhighlight lang="bash">
> 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
</syntaxhighlight>

== Расчет площадей по отдельным участкам ==
Однако, наша задача несколько иная: нужно расчитать площади по каждому из заповедников в отдельности. Для того, чтобы это сделать, необходимо следующее:
#Извлечь из карты oopt интересующий нас полигон.
#Создать маску из этого полигона, чтобы команда "r.surf.area relief" производила расчеты только в замаскированной области
#Вызов команды r.surf.area relief и сохранение результата

Реализация:
Посмотрим содержимое таблицы атрибутов карты oopt:
<syntaxhighlight lang="bash">
> v.db.select oopt
cat|NAME_PRT_R|NAME_R|TYPE
1|Катунский|Катунский|Заповедник
2|Кузнецкий Алатау|Кузнецкий Алатау|Заповедник
3|Участок "Ханхаринский"|Тигирекский|Заповедник
4|Участок "Белорецкий"|Тигирекский|Заповедник
5|Участок "Тигирекский"|Тигирекский|Заповедник
6|Алтайский|Алтайский|Заповедник
7|Шорский|Шорский|Национальный
8|Кирзинский|Кирзинский|Заказник
</syntaxhighlight>

Создадим из векторной карты oopt растровую:
<syntaxhighlight lang="bash">
> v.to.rast in=oopt out=oopt use=attr col=cat label=NAME_PRT_R
</syntaxhighlight>

Посмотрим, что из этого получилось, заодно рассчитаем статистику по площадям:
<syntaxhighlight lang="bash">
> 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|
+-----------------------------------------------------------------------------+
</syntaxhighlight>

Создадим маску, соответствующую участу "Катунский" (категория 1):
<syntaxhighlight lang="bash">
> r.mask oopt maskcat=1
</syntaxhighlight>


Расчитаем площадь для этого участка:
<syntaxhighlight lang="bash">
> 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.
</syntaxhighlight>


Таким же образом можно получить информацию и по остальным участкам.

Эту процедуру легко автоматизировать, сохранив данные расчетов в файл, например, area.txt:
<syntaxhighlight lang="bash">
> 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
</syntaxhighlight>


Начальные строки файла area.txt представлены ниже:
<syntaxhighlight lang="bash">
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.
</syntaxhighlight>