Просмотр кода страницы Анализ данных с использованием GRASS GIS и R

{{Статья|Опубликована|grass-r}}
{{Аннотация|Примеры совместного использования ГИС GRASS и cтатпакета и языка R}}


ГИС GRASS является мощной геоинформационной системой с открытым исходным кодом, предназначенной для управления пространственными данными, обработки изображений (в том числе данных ДЗЗ), пространственного моделирования, визуализации данных и т.д. ГИС GRASS предоставляет пользователю множество модулей и функций, облегчающих анализ и обработку данных. R представляет собой платформу с открытым исходным кодом, предназначенную для статистических вычислений, моделирования и анализа. В данной статье рассматриваются примеры совместного использования R и ГИС GRASS.

Цель этой статьи - дать читателю отправную точку для дальнейшего изучения и использования связки GRASS GIS и R. Предполагается, что читатель имеет некоторое знакомство с данными системами. В частности, читатель должен иметь представление о синтаксисе языка R и иметь хотя бы небольшой опыт работы с ГИС GRASS.

Хотя для экспериментов мы будем использовать [http://gis-lab.info/qa/geosample.html набор данных geosample], жесткой привязки к данному набору данных нет, и читатель с легкостью может адаптировать примеры под свои нужды и данные.

В статье использованы материалы книги [http://www.grassbook.org/ "Open Source GIS: A GRASS GIS Approach"] (авторы Markus Neteler, Helena Mitasova).

==Установка необходимых пакетов==

Предположим, что у нас есть база геоданных GRASS (location/mapset), которую мы хотим прочитать в R, определенным образом обработать/проанализировать и сохранить полученный результат обратно в базу геоданных. Для решения этой задачи существует несколько подходов, мы воспользуемся тем, что и GRASS GIS, и R предоставляют пользователю командную оболочку, в которой пользователь может вводить команды и немедленно получать результат. Таким образом, режим работы будет следующий:

* Заходим в нужную область (location/mapset) GRASS GIS;
* Из командной строки GRASS вызываем среду R;
* Производим требуемые манипуляции с геоданными GRASS, используя возможности анализа среды R;
* При необходимости сохраняем результаты;
* Выходим из R и GRASS.

Таким образом, большую часть времени мы будем находиться одновременно в двух системах: GRASS и R. Однако, чтобы такого рода "матрёшка" сработала, необходимо, чтобы в среде R были установлены пакеты, позволяющие читать данные, хранящиеся в формате GRASS.

Для работы нам понадобится пакет spgrass6 и его зависимости (однако для того, чтобы производить анализ пространственных данных, рекомендуется посмотреть также пакеты akima, fields, geoR, grid, gstat, lattice, MASS, scatterplot3d, spatial и stepfun).

Установка пакетов в среде R не должна вызвать каких-либо трудностей, например для установки пакетов spgrass6 и gstat достаточно выполнить следующие команды (требуется подключение к интернет):

 install.packages("spgrass6", "gstat", dependencies = TRUE)

Загрузить необходимую библиотеку (например, spgrass6) в средy R можно командой

 library(spgrass6)

==Чтение и запись данных GRASS в среде R==

===Первые шаги===

Запустим GRASS и из командной строки GRASS загружаем R. В конечном итоге перед нами окажется командная строка системы R.

Загрузим в R необходимые библиотеки:

 library(spgrass6)

После загрузки данной библиотеки уже можно взаимодействовать с GRASS GIS. Например, введем следующую команду:

 str(gmeta6())

Результат выполнения команды приводится ниже:

 
 &gt; str(gmeta6())
 List of 24
  $ GISDBASE     : chr "/home/dima/GIS/grass"
  $ LOCATION_NAME: chr "grass"
  $ MAPSET       : chr "PERMANENT"
  $ MONITOR      : chr "x0"
  $ GRASS_GUI    : chr "text"
  $ n            : num 57.3
  $ s            : num 49
  $ w            : num 75
  $ e            : num 90
  $ t            : num 1
  $ b            : num 0
  $ nsres        : num 0.00417
  $ nsres3       : num 1.04
  $ ewres        : num 0.00417
  $ ewres3       : num 1
  $ tbres        : num 1
  $ rows         : int 1992
  $ rows3        : int 8
  $ cols         : int 3600
  $ cols3        : int 15
  $ depths       : int 1
  $ cells        : chr "7171200"
  $ 3dcells      : chr "120"
  $ proj4        : chr "+proj=longlat +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000"
  - attr(*, "class")= chr "gmeta6"

Как легко заменить, данная команда вывела информацию о текущей области (region) GRASS и другие параметры.

Находясь в R, вы имеете возможность запускать команды GRASS, не выходя из сессии R. Для этой цели можно использовать команду system, например:

 
 &gt; system("g.region -p")
 projection: 3 (Latitude-Longitude)
 zone:       0
 datum:      wgs84
 ellipsoid:  wgs84
 north:      57:18N
 south:      49N
 west:       75E
 east:       90E
 nsres:      0:00:15
 ewres:      0:00:15
 rows:       1992
 cols:       3600
 cells:      7171200

===Импорт векторных данных===

Импортируем несколько карт из GRASS GIS в среду R, для этого используется команда readVECT6.

 admin &lt;- readVECT6("admin", ignore.stderr=TRUE)

Просмотрим информацию:

 
 &gt; summary(admin)
 Object of class SpatialPolygonsDataFrame
 Coordinates:
        min      max
 x 75.08194 89.86998
 y 49.08308 57.25027
 Is projected: FALSE 
 proj4string :
 [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0]
 Data attributes:
       cat                          NAME  
  Min.   :1.00   Алтайский край       :1  
  1st Qu.:1.75   Новосибирская область:1  
  Median :2.50   Республика Алтай     :2  
  Mean   :2.50                            
  3rd Qu.:3.25                            
  Max.   :4.00

Аналогично, загрузим данные по почвам:

 
 &gt; soils &lt;- readVECT6("soils", ignore.stderr=TRUE) 
 &gt; summary(soils)
 Object of class SpatialPolygonsDataFrame
 Coordinates:
        min      max
 x 75.08194 89.86998
 y 49.08308 57.25027
 Is projected: FALSE 
 proj4string :
 [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0]
 Data attributes:
       cat         FAOSOIL         UNIT                   TYPE           NAME   
  Min.   :  1   Ch1-3a : 16   I      :39   LITHOSOLS        :39   Черноземы:52  
  1st Qu.: 44   C1-3a  : 14   Ch     :17   Haplic Chernozems:17   Литосоли :39  
  Median : 87   Sm13-3a: 14   C      :16   CHERNOZEMS       :16   Солонцы  :18  
  Mean   : 87   Oe1-a  : 10   Sm     :15   Mollic Solonetz  :15   Сероземы :11  
  3rd Qu.:130   Mo1-2ab:  9   Ck     :10   Calcic Chernozems:10   Гистосоли:10  
  Max.   :173   (Other):110   (Other):76   (Other)          :76   (Other)  :43  
  NA's   :  1   NA's   :  1   NA's   : 1   NA's             : 1   NA's     : 1  
    FAOSOIL_L    
  Min.   : 1.00  
  1st Qu.: 8.00  
  Median :21.00  
  Mean   :20.26  
  3rd Qu.:34.00  
  Max.   :40.00  
  NA's   : 1.00

При необходимости можно отобразить данные на карте:

 plot(admin, axes=TRUE)

Добавим еще один слой:

 plot(soils, add=TRUE, col="red")

===Импорт растровых данных===

Импорт растровых данных производится аналогично. Для этого предназначена команда readRAST6:

 readRAST6(vname, cat=NULL, ignore.stderr = FALSE, NODATA=NULL, plugin=NULL, mapset=NULL, useGDAL=TRUE, close_OK=TRUE)

В первую очередь, нас будет интересовать параметр vname - название растровой карты GRASS. Считаем, к примеру, растр под названием modis:

 
 &gt; modis &lt;- readRAST6("modis", ignore.stderr=TRUE)
 &gt; summary(modis)
 Object of class SpatialGridDataFrame
 Coordinates:
   min  max
 x  75 90.0
 y  49 57.3
 Is projected: FALSE 
 proj4string :
 [+proj=longlat +a=6378137 +rf=298.257223563 +no_defs
 +towgs84=0.000,0.000,0.000]
 Number of points: 2
 Grid attributes:
   cellcentre.offset    cellsize cells.dim
 x          75.00208 0.004166667      3600
 y          49.00208 0.004166667      1992
 Data attributes:
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
       0    4328    5418    7028    7499   32770  209650 
 

===Экспорт данных в GRASS GIS===

Для экспорта растровых данных существует команда:

 writeRAST6(x, vname, zcol = 1, NODATA=NULL, ignore.stderr = FALSE, useGDAL=TRUE, overwrite=FALSE, flags=NULL)

Здесь x - фрейм пространственных данных, который будет сохранен как растровая карта GRASS под именем vname (подробнее см. в справке R).

Аналогично, для работы с векторными данными используется команда:

 writeVECT6(SDF, vname,  v.in.ogr_flags=NULL, ignore.stderr = FALSE)

Здесь SDF - фрейм пространственных данных, который будет сохранен как векторная карта GRASS под именем vname (подробенее см. в справке R).

==Пример обработки данных==

В наборе geosample содержатся несколько растровых изображений: это данные MODIS (растр modis) и данные рельефа (растр relief). В качестве примера построим регрессию зависимости яркости modis от высоты relief.

Сначала загрубим разрешение, чтобы не обрабатывать излишнее количество данных, сделаем это в GRASS (хотя это же можно и через R):

  
 &gt; g.region res=0:01:0 -p
 projection: 3 (Latitude-Longitude)
 zone:       0
 datum:      wgs84
 ellipsoid:  wgs84
 north:      57:18N
 south:      49N
 west:       75E
 east:       90E
 nsres:      0:01
 ewres:      0:01
 rows:       498
 cols:       900
 cells:      448200

Заходим в R, подключаем библиотеки, импортируем данные:

 
 &gt; library(spgrass6)
 &gt; modis &lt;- readRAST6("modis", ignore.stderr=TRUE)
 &gt; relief &lt;- readRAST6("relief", ignore.stderr=TRUE)

Приведем данные к числовому виду, понятному функции для линейных регрессий lm:

 
 &gt; relief_mat &lt;- as.matrix(relief)
 &gt; relief_vec &lt;- as.vector(relief_mat)
 &gt; modis_mat &lt;- as.matrix(modis)
 &gt; modis_vec &lt;- as.vector(modis_mat)

Построим линейную регрессию:

  
 &gt; model &lt;- lm(modis_vec ~ relief_vec)
 &gt; summary(model)
 
 Call:
 lm(formula = modis_vec ~ relief_vec)
 
 Residuals:
      Min       1Q   Median       3Q      Max 
 -15950.9  -2035.4   -182.4   1276.3  20414.5 
 
 Coefficients:
              Estimate Std. Error t value Pr(&gt;|t|)    
 (Intercept) 4.625e+03  7.569e+00   611.1   &lt;2e-16 ***
 relief_vec  4.916e+00  9.682e-03   507.7   &lt;2e-16 ***
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 
 Residual standard error: 3893 on 433741 degrees of freedom
   (14457 observations deleted due to missingness)
 Multiple R-squared: 0.3728,	Adjusted R-squared: 0.3728 
 F-statistic: 2.578e+05 on 1 and 433741 DF,  p-value: &lt; 2.2e-16 
 

Интересно сравнить результаты, которые мы получили в R с регрессией, которую мы построим в самой ГИС GRASS. Для этого вызовем команду GRASS r.regression.line:

 
 &gt; system("r.regression.line map1=modis map2=relief")
  100%
 y = a + b*x
    a: offset
    b: gain
    R: sumXY - sumX*sumY/tot
    N: number of elements
    medX, medY: Means
    sdX, sdY: Standard deviations
 a  b  R  N  F medX  sdX  medY  sdY
 -44.4697 0.0758367 0.610568 433743 -0.372793 7025.41 4915.05 488.314 610.483

Построим теперь квадратичную регрессию, для вычисления которой в GRASS GIS нет соответствующего модуля:

 
 &gt; model2 &lt;- lm(relief_vec ~ modis_vec + I(modis_vec^2))
 &gt; summary(model2)
 
 Call:
 lm(formula = relief_vec ~ modis_vec + I(modis_vec^2))
 
 Residuals:
      Min       1Q   Median       3Q      Max 
 -1930.25  -245.41  -157.72    88.94  2803.82 
 
 Coefficients:
                  Estimate Std. Error t value Pr(&gt;|t|)    
 (Intercept)     3.340e+02  2.308e+00  144.72   &lt;2e-16 ***
 modis_vec      -1.310e-02  4.810e-04  -27.23   &lt;2e-16 ***
 I(modis_vec^2)  3.351e-06  1.730e-08  193.70   &lt;2e-16 ***
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 
 Residual standard error: 463.8 on 433740 degrees of freedom
   (14457 observations deleted due to missingness)
 Multiple R-squared: 0.4227,	Adjusted R-squared: 0.4227 
 F-statistic: 1.588e+05 on 2 and 433740 DF,  p-value: &lt; 2.2e-16 
 

Таким образом, на нескольких простых примерах было показан основной цикл работы в двух средах GRASS GIS и R.