# глобальные опции отображения графиков
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 точка Москвы')