Каков наилучший способ построить более одного растра с большими различиями между растрами в R?

Я пытался разработать лучшее решение для визуализации моих данных с помощью R. Однако мои усилия пока не увенчались успехом.

Предположим, у вас есть 6 растров (данные относятся к одной переменной, измеренной в разное время), где внутрирастровая вариация не так велика, а между растрами больше.

Вам нужна одна цветовая полоса для 6 растров. Что мне нужно, так это отобразить 6 растров, по-прежнему иметь возможность видеть вариации растров внутри и между ними и создать единый цветовой ключ с метками разрыва.

Одним из предложений является изменение масштаба данных и создание общей шкалы. Что-то вроде:

require(raster)
require(rasterVis)
s = stack(r1,...,r6)
s = (s-cellStats(s,"min"))/(cellStats(s,"max")-cellStats(s,"min"))

После масштабирования данных вариации между растрами уменьшаются и выглядят лучше.

Кроме того, как я могу создать легенду цветовой полосы, которая использует zlim=(min(s),max(s)), s как немасштабированные данные?

Большой воспроизводимый пример:

r1

r2

r3

r4

r5

r6


r maps raster r-raster
person code123    schedule 31.03.2015    source источник
comment
Пожалуйста, не могли бы вы опубликовать воспроизводимый пример?   -  person Oscar Perpiñán    schedule 31.03.2015
comment
Если вам нужно определить собственный цветовой ключ, используйте аргумент colorkey, как описано в этом разделе часто задаваемых вопросов. .   -  person Oscar Perpiñán    schedule 31.03.2015
comment
Спасибо за правки, Оскар. Конечно, благодаря вашему вкладу я могу определить разрывы и проделать большую часть работы с использованием растра и rasterVis. Как визуализировать данные, с которыми мы часто сталкиваемся в исследованиях климата, остается дилеммой. Если проблема, которую я здесь поднимаю, решена, тогда можно лучше использовать инструменты сопоставления в R. Я отредактировал вопрос, включив в него большой воспроизводимый пример через Dropbox.   -  person code123    schedule 31.03.2015

Ответы (1)


arrow_upward
0
arrow_downward

Это то, что у меня есть до сих пор:

library(hydroTSM)

x.tsZ<-as.numeric(precip[,1])# gives r1 points for interpolation etc

names(x.tsZ)<-EbroPPgisZ$ID

x.idw <- interpo(x.ts= x.tsZ, x.gis=EbroPPgisZ, 
                    X="EAST_ED50", Y="NORTH_ED50", sname="ID", bname="CHE_BASIN_NAME",
                    type= "cells",p4s= projjj,
                    subcatchments= EbroCatchmentsCHE,
                    cell.size= 18000, col.nintv=100,col.at=at,
                    ColorRamp= "PCPAnomaly2",  
                    main= "IDW Precipitation on the Ebro")


## Removing unncesary field
x.idw@data["var1.var"]  <- NULL

# Setting the projection of the interpolated object (ED50 Z30N)
proj4string(x.idw) <- projjj

#library(raster)
x.idw.r <- raster(x.idw)

# Defining the LatLong spatial reference
wgs84.p4s <- CRS("+proj=longlat +datum=NAD83 +ellps=GRS80 +no_defs")

# Re-projecting the from ED Z30N to LatLong WGS84.
x.idw.r.ll <- projectRaster(from=x.idw.r, crs=wgs84.p4s@projargs, method="ngb")

# Saving the raster file into a GeoTiff
#####                                                                     #####
##### ATENTION : The following line will WRITE A FILE into your hard disk #####
#####                                                                     #####
writeRaster(x.idw.r.ll, file="myfile1.tif", format="GTiff", overwrite=TRUE)

# display raster alongside boundaries
ras = raster("myfile1.tif")

cr <- crop(ras, extent(Prairie.Boundaries), snap="out") 
fr <- rasterize(Prairie.Boundaries, cr) 
r1  <- mask(x=cr, mask=fr)
# repeat the above for r2...26


    myat =unique(seq(10.87775,  16.74664,length.out=col.nintv))#col.nintv=45
    myat2=round(myat,digits = 1)# 
    a=unique(myat2)
    a[seq(0, length(a), 5)]
    myat1=c(11.3, 11.8, 12.3, 12.8, 13.3, 13.8, 14.3, 14.8, 15.3, 15.8, 16.3)# output of "a" defines where to place labels

    themes2 <- themes<- colorRampPalette(c("sienna4", "peachpuff","orange", "yellow", "lightskyblue1", "royalblue3", "darkblue"))#(length(myat)-1)

    myColorkey <- list(at=myat,space = "bottom",labels=list(axis.line = list(col = NA),cex=1,font=1,at=myat1,rot=0),height=1,width=1)#also control lengend with
    #coordinates(SITES...TTEST) <- ~x+y
    #detach ggplot2 before running this

    if (dev.cur() == 1) x11(width=10,height=6)
    s <- stack(r1,r2,r3,r4,r5,r6)
    levelplot(s, layout=c(3, 2), index.cond=list(c(1, 3, 5, 2, 4, 6)),col.regions=themes2,margin=FALSE,xlab=NULL,ylab=NULL,contour=F,
              par.strip.text=list(cex=0),scales=list(x=list(draw=FALSE),y=list(draw=FALSE)),colorkey=myColorkey,at=myat)+# add points to map

      layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_JJA1 < 0.05, 3, 1), cex=.4, col=138), rows=1,columns=1) +layer(sp.polygons(watersheds,lwd=0.5))+# add polygons to map
      layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_JJA2 < 0.05, 3, 1), cex=.4, col=138), rows=2,columns = 1)+
      layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_JJA3 < 0.05, 3, 1), cex=.4, col=138), rows=1,columns = 2)+
      layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_JJA4 < 0.05, 3, 1), cex=.4, col=138), rows=2,columns = 2)+
      layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_JJA5 < 0.05, 3, 1), cex=.4, col=138), rows=1,columns = 3)+
      layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_JJA6 < 0.05, 3, 1), cex=.4, col=138), rows=2,columns = 3)

введите здесь описание изображения

person code123    schedule 02.04.2015