# глобальные опции отображения графиков
knitr::opts_chunk$set(fig.show = 'hold')
options(digits = 2)
Цель - нанести на спутниковую карту Москвы уровни освещённости (данные 2009 года).
В работе использованы пакеты:
library('R.utils') # gunzip() для распаковки архивов
library('dismo') # gmap() для загрузки Google карты
library('raster') # функции для работы с растровыми картами в R
library('OpenStreetMap') # работа с геоданными от OpenStreetMap
library('knitr') # отчёты в html, pdf, doc
Исходные данные:
* Данные по ночной освещённости, собранные Национальным управлением океанических и атмосферных исследований США и размещённые на сайте ngdc.noaa.gov. Эти данные загружаются в виде архива (11,2 Мб в архиве, почти 1 Гб после распаковки).
* Спутниковая карта мира от Google Earth, нужный участок которой загружен с помощью функции gmap()
пакета «dismo
».
# директория для данных
if(!file.exists('./data')) dir.create('data')
# загрузка карты ночной освещённости мира
url_radianceCalibrated <- 'ftp://ftp.ngdc.noaa.gov/DMSP/web_data/x_rad_cal/rad_cal.tar'
calibratedLights <- './data/rad_cal.tar'
hiResTif <- './data/world_avg.tif'
# скачать и распаковать данные
if(!file.exists(calibratedLights)) download.file(url_radianceCalibrated,
calibratedLights, mode = 'wb')
if(!file.exists(paste(hiResTif, '.gz', sep=''))) untar(calibratedLights,
exdir = './data')
if(!file.exists(hiResTif)) gunzip(paste(hiResTif, '.gz', sep=''))
# создать растровый объект из графического файла
r <- raster(hiResTif)
# выбрать область на карте, ограниченную координатами
# Москва
long.min <- 36.95 # восточная долгота, левая граница карты
long.max <- 38.45 # восточная долгота, правая граница карты
lat.min <- 55.35 # северная широта, нижняя граница карты
lat.max <- 56.10 # северная широта, верхняя граница карты
# делаем 'рамку' для участка карты
e <- extent(long.min, long.max, lat.min, lat.max)
# обрезаем карту освещённости по рамке
rc <- crop(r, e)
# загружаем спутниковую карту и тоже обрезаем
g <- gmap(e, type='satellite')
Нанесём на спутниковую карту Москвы, ограниченную 36.95 в.д. и 38.45 в.д., 55.35 с.ш. и 56.1 с.ш., уровни освещённости.
# пересчитываем карту освещённости в проекцию спутниковой карты
mrc <- projectRaster(rc, g)
# строим график освещённости
# чтобы избежать появления пустого графика,
# сохраняем график во внешний файл
# папка для рисунков к отчёту
if(!file.exists('./lab2-2_pic')) dir.create('./lab2-2_pic')
png('./lab2-2_pic/RPlot1.png')
# 1. сама карта
plot(g, interpolate = TRUE)
# 2. контуры уровней освещённости
contour(mrc, add = TRUE, lwd = 2, levels = 1:5*20, col = rev(heat.colors(5)))
dev.off()
Максимальная освещённость равна 122. Найдём её координаты.
# координаты самой освещённой точки
xy <- xyFromCell(rc, which.max(rc))
xy
## x y
## [1,] 38 56
Нанесём точку на карту.
mxy <- Mercator(xy)
if(!file.exists('./lab2-2_pic')) dir.create('./lab2-2_pic')
png('./lab2-2_pic/RPlot2.png')
# 1. сама карта
plot(g, interpolate = TRUE)
# 2. контуры уровней освещённости
contour(mrc, add = TRUE, lwd = 2, levels = 1:5*20, col = rev(heat.colors(5)))
# 3. наносим точку на карту
points(mxy, col = 'red', bg = 'red', pch = 21, cex = 0.5)
dev.off()
Построим карту уровня улиц в окресностях наиболее освещённой точки.
# типы карт
type.map <- c('osm', 'maptoolkit-topo', 'mapquest', 'bing',
'stamen-toner', 'stamen-watercolor',
'osm-german', 'esri', 'esri-topo',
'skobbler', 'opencyclemap', 'osm-transport',
'osm-public-transport')
# вектор с координатами самой освещённой точки
xy <- as.vector(xy)
# загрузка карты
map <- openmap(upperLeft = c(xy[2] + 0.008, xy[1] - 0.0175),
lowerRight = c(xy[2] - 0.008, xy[1] + 0.0175),
type = type.map[2])
plot(map, interpolate = T)
# добавляем точку и подпись к ней
points(Mercator(xy), col = 'blue', pch = '+', cex = 2)
text(Mercator(xy + 0.0007), col = 'blue', font = 3, cex = 0.8,
'Самая освещѐнная \n точка Москвы')