Я использовал ggplot2 для построения климатических данных с координатной сеткой в течение многих лет. Обычно это проецируемые файлы NetCDF. Ячейки имеют квадратные координаты в координатах модели, но в зависимости от того, какую проекцию использует модель, в реальном мире может быть не так.
Мой обычный подход - сначала переназначить данные на подходящую регулярную сетку, а затем построить. Это приводит к небольшому изменению данных, обычно это приемлемо.
Однако я решил, что этого уже недостаточно: я хочу построить спроецированные данные напрямую, без переназначения, как могут, если я не ошибаюсь, делать другие программы (например, ncl), не касаясь выходных значений модели.
Однако я столкнулся с некоторыми проблемами. Ниже я подробно расскажу о возможных решениях, от самых простых до самых сложных, и их проблемах. Сможем ли мы их преодолеть?
РЕДАКТИРОВАТЬ: благодаря ответу @ lbusett я получил эту прекрасную функцию, которая включает решение. Если вам это нравится, проголосуйте за ответ @ lbusett!
Начальная настройка
#Load packages
library(raster)
library(ggplot2)
#This gives you the starting data, 's'
load(url('https://files.fm/down.php?i=kew5pxw7&n=loadme.Rdata'))
#If you cannot download the data, maybe you can try to manually download it from http://s000.tinyupload.com/index.php?file_id=04134338934836605121
#Check the data projection, it's Lambert Conformal Conic
projection(s)
#The data (precipitation) has a 'model' grid (125x125, units are integers from 1 to 125)
#for each point a lat-lon value is also assigned
pr <- s[[1]]
lon <- s[[2]]
lat <- s[[3]]
#Lets get the data into data.frames
#Gridded in model units:
pr_df_basic <- as.data.frame(pr, xy=TRUE)
colnames(pr_df_basic) <- c('lon', 'lat', 'pr')
#Projected points:
pr_df <- data.frame(lat=lat[], lon=lon[], pr=pr[])
Мы создали два фрейма данных, один с координатами модели, другой с реальными точками пересечения широты и долготы (центрами) для каждой ячейки модели.
Необязательно: использование меньшего домена
Если вы хотите более четко видеть формы ячеек, вы можете подгруппировать данные и извлечь только небольшое количество ячеек модели. Просто будьте осторожны, вам может потребоваться настроить размер точек, пределы сюжета и другие удобства. Вы можете подмножество, подобное этому, а затем повторить приведенную выше часть кода (без load()):
s <- crop(s, extent(c(100,120,30,50)))
Если вы хотите полностью разобраться в проблеме, возможно, вы захотите попробовать как большой, так и малый домен. код идентичен, меняются только размеры точек и лимиты карты. Приведенные ниже значения относятся к большому целому домену. Хорошо, а теперь приступим к сюжету!
Начните с плитки
Наиболее очевидное решение - использовать плитку. Давай попробуем.
my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())
my_cols <- scale_color_distiller(palette='Spectral')
my_fill <- scale_fill_distiller(palette='Spectral')
#Really unprojected square plot:
ggplot(pr_df_basic, aes(y=lat, x=lon, fill=pr)) + geom_tile() + my_theme + my_fill
И вот результат: 
Хорошо, теперь кое-что более продвинутое: мы используем настоящие LAT-LON, используя квадратные плитки.
ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.2, height=1.2) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat)) #the result is weird boxes...
Хорошо, но это не настоящие квадраты модели, это хитрость. Кроме того, блоки модели расходятся в верхней части области и все ориентированы одинаково. Не хорошо. Давайте спроецируем сами квадраты, хотя мы уже знаем, что это неправильно ... может быть, это выглядит хорошо.
#This takes a while, maybe you can trust me with the result
ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.5, height=1.5) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))
Во-первых, на это уходит много времени. Недопустимо. И еще раз: это неправильные модельные ячейки.
Попробуйте использовать очки, а не плитки
Может быть, мы можем использовать круглые или квадратные точки вместо плиток, а также проецировать их!
#Basic 'unprojected' point plot
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) + geom_point(size=2) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_cols + my_theme +
coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat))
Мы можем использовать квадратные точки ... и проецировать! Мы подходим ближе, хотя знаем, что это все равно неверно.
#In the following plot pointsize, xlim and ylim were manually set. Setting the wrong values leads to bad results.
#Also the lambert projection values were tired and guessed from the model CRS
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) +
geom_point(size=2, shape=15) +
borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_cols +
coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))
Достойные результаты, но не полностью автоматические, и построения точек недостаточно. Я хочу, чтобы клетки реальной модели с их формой были видоизменены проекцией!
Может быть, полигоны?
Итак, как вы можете видеть, я ищу способ правильно построить блоки модели, спроектированные в правильной форме и положении. Конечно, блоки модели, которые представляют собой квадраты в модели, после проецирования становятся формами, которые больше не являются правильными. Так, может быть, я смогу использовать полигоны и проецировать их? Я пытался использовать rasterToPolygons и fortify и следить за этим сообщением, но мне это не удалось. Я пробовал это:
pr2poly <- rasterToPolygons(pr)
#http://mazamascience.com/WorkingWithData/?p=1494
pr2poly@data$id <- rownames(pr2poly@data)
tmp <- fortify(pr2poly, region = "id")
tmp2 <- merge(tmp, pr2poly@data, by = "id")
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill
Ладно, попробуем заменить латинские на латинские ...
tmp2$long <- lon[]
tmp2$lat <- lat[]
#Mh, does not work! See below:
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill
(извините, что поменял цветовую гамму в сюжетах)
Мммм, даже не стоит пробовать с проекцией. Может, мне стоит попытаться вычислить широты и долготы углов ячеек модели, создать для них многоугольники и перепроецировать их?
Вывод
- Я хочу нанести данные проектируемой модели на ее собственную сетку, но у меня не получилось. Использование тайлов неверно, использование точек - хакерское занятие, а использование полигонов, похоже, не работает по неизвестным причинам.
- При проецировании через
coord_map()линии сетки и метки осей неправильные. Это делает спроектированные ggplots непригодными для публикаций.
ggmap? - person ulfelder schedule 25.04.2017readAll()перед сохранением не делал. Теперь он должен работать, вы можете его протестировать? - person AF7 schedule 27.04.2017sрезультатов как"+proj=lcc +lat_1=30. +lat_2=65. +lat_0=48. +lon_0=9.75 +x_0=-25000. +y_0=-25000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs", что является метрической проекцией. Однако ваши координаты - долгота-широта, а разрешение растра - 1. Где вы взяли набор данных? Делали ли вы над ним какие-то пространственные разработки? - person lbusett schedule 27.04.2017+units=m, и я не эксперт по прогнозам, возможно, для этого есть причина. Странный. Тем не менее, похоже, что это сработает, если я перепроектирую материал. - person AF7 schedule 28.04.2017mapview, вы получите совершенно неправильную систему координат. - person lbusett schedule 28.04.2017RasterStack. Это может привести к путанице, поскольку, например, стандартные процедуры обработкиraster, такие какprojectRaster, не будут работать должным образом. На мой взгляд, было бы лучше сохранить его какSpatialPixelsDataFrame. - person lbusett schedule 28.04.2017raster. Если проекция необходима внутри R, я следую этим отличным инструкциям: gis.stackexchange.com/questions/120900/ - person AF7 schedule 28.04.2017NCLможет делать более точные графики, это меня очень раздражает :) С другой стороны, другие программы для построения графиков NetCDF (такие как GrADS и, как мне кажется, Ferret) просто автоматически переназначаются на свой внутренний максимум. -res сетка LAT-LON, не сообщая об этом пользователю. Эти сюжеты считаются приемлемыми, но я хочу добиться большего. - person AF7 schedule 29.04.2017