ользовал с помощьюggplot2 построить климатические данные с привязкой к сетке за годы. Обычно это проецируемые файлы NetCDF. Ячейки имеют квадратные координаты модели, но в зависимости от того, какую проекцию использует модель, в реальном мире это может быть не так.

Мой обычный подход - сначала переназначить данные на подходящую регулярную сетку, а затем построить график. Это вносит небольшую модификацию в данные, обычно это приемлемо.

Тем не менее, я решил, что это уже недостаточно хорошо: я хочу отобразить проецируемые данные напрямую, без переотображения, как другие программы (например,ncl) можно, если я не ошибаюсь, сделать, не касаясь выходных значений модели.

Однако я сталкиваюсь с некоторыми проблемами. Ниже я подробно опишу возможные решения, от самых простых до самых сложных, и их проблемы. Можем ли мы их преодолеть?

РЕДАКТИРОВАТЬ: благодаря ответу @ lbusett я получилэта хорошая функция это включает в себя решение. Если вам это нравится, пожалуйста, upvote@ lbusett's answer!Начальная настройка

#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 непригодными для публикаций.