0

我一直在尝试开发一种更好的解决方案来使用 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

4

1 回答 1

0

这是我到目前为止所拥有的：

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)

在此处输入图像描述

于 2015-04-02T01:51:57.327 回答