Я использовал ggplot2 для построения ggplot2 климатических данных в течение многих лет. Обычно это проецируемые файлы NetCDF. Ячейки имеют квадратные координаты модели, но в зависимости от того, какую проекцию использует модель, в реальном мире это может быть не так.
Мой обычный подход - сначала переназначить данные на подходящую регулярную сетку, а затем построить график. Это вносит небольшую модификацию в данные, обычно это приемлемо.
Тем не менее, я решил, что это уже недостаточно хорошо: я хочу ncl спроецированные данные напрямую, без переотображения, как это могут сделать другие программы (например, ncl), если я не ошибаюсь, не касаясь выходных значений модели.
Однако я сталкиваюсь с некоторыми проблемами. Ниже я подробно опишу возможные решения, от самых простых до самых сложных, и их проблемы. Можем ли мы их преодолеть?
РЕДАКТИРОВАТЬ: благодаря ответу @lbusett я получил эту прекрасную функцию, которая включает в себя решение. Если вам это нравится, пожалуйста, ответьте upvote @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 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
[email protected]$id <- rownames([email protected])
tmp <- fortify(pr2poly, region = "id")
tmp2 <- merge(tmp, [email protected], 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 непригодными для публикаций.
