Анализ данных с использованием GRASS GIS и R
по адресу http://gis-lab.info/qa/grass-r.html
Примеры совместного использования ГИС GRASS и cтатпакета и языка R
ГИС GRASS является мощной геоинформационной системой с открытым исходным кодом, предназначенной для управления пространственными данными, обработки изображений (в том числе данных ДЗЗ), пространственного моделирования, визуализации данных и т.д. ГИС GRASS предоставляет пользователю множество модулей и функций, облегчающих анализ и обработку данных. R представляет собой платформу с открытым исходным кодом, предназначенную для статистических вычислений, моделирования и анализа. В данной статье рассматриваются примеры совместного использования R и ГИС GRASS.
Цель этой статьи - дать читателю отправную точку для дальнейшего изучения и использования связки GRASS GIS и R. Предполагается, что читатель имеет некоторое знакомство с данными системами. В частности, читатель должен иметь представление о синтаксисе языка R и иметь хотя бы небольшой опыт работы с ГИС GRASS.
Хотя для экспериментов мы будем использовать набор данных geosample, жесткой привязки к данному набору данных нет, и читатель с легкостью может адаптировать примеры под свои нужды и данные.
В статье использованы материалы книги "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())
Результат выполнения команды приводится ниже:
> 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, например:
> 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 <- readVECT6("admin", ignore.stderr=TRUE)
Просмотрим информацию:
> 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
Аналогично, загрузим данные по почвам:
> soils <- readVECT6("soils", ignore.stderr=TRUE)
> 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:
> modis <- readRAST6("modis", ignore.stderr=TRUE)
> 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):
> 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, подключаем библиотеки, импортируем данные:
> library(spgrass6)
> modis <- readRAST6("modis", ignore.stderr=TRUE)
> relief <- readRAST6("relief", ignore.stderr=TRUE)
Приведем данные к числовому виду, понятному функции для линейных регрессий lm:
> relief_mat <- as.matrix(relief) > relief_vec <- as.vector(relief_mat) > modis_mat <- as.matrix(modis) > modis_vec <- as.vector(modis_mat)
Построим линейную регрессию:
> model <- lm(modis_vec ~ relief_vec)
> 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(>|t|)
(Intercept) 4.625e+03 7.569e+00 611.1 <2e-16 ***
relief_vec 4.916e+00 9.682e-03 507.7 <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: < 2.2e-16
Интересно сравнить результаты, которые мы получили в R с регрессией, которую мы построим в самой ГИС GRASS. Для этого вызовем команду GRASS r.regression.line:
> 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 нет соответствующего модуля:
> model2 <- lm(relief_vec ~ modis_vec + I(modis_vec^2))
> 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(>|t|)
(Intercept) 3.340e+02 2.308e+00 144.72 <2e-16 ***
modis_vec -1.310e-02 4.810e-04 -27.23 <2e-16 ***
I(modis_vec^2) 3.351e-06 1.730e-08 193.70 <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: < 2.2e-16
Таким образом, на нескольких простых примерах было показан основной цикл работы в двух средах GRASS GIS и R.