Оси на минимальной протяженности, без отступов, на графиках растровых * объектов

Question

Есть ли способ убедиться, что рамка вокруг графика точно соответствует растровым экстентам? Далее в зависимости от пропорций устройства имеется зазор сверху и снизу или слева и справа от растра:

require(raster)
r = raster()
r[]= 1
plot(r, xlim=c(xmin(r), xmax(r)), ylim=c(ymin(r), ymax(r)))

Одним из элементов проблемы с растровыми объектами является то, что asp=1 обеспечивает правильное отображение. Следующая базовая диаграмма рассеяния имеет ту же проблему, когда asp=1:

plot(c(1:10), c(1:10), asp=1)

Попробуйте vectorplot(r) из пакета rasterVis, чтобы увидеть, как я хочу, чтобы оси выглядели.

EDIT:

Решения должны хорошо играть с наложениями SpatialPoints, а не отображать точки за пределами указанного растра:

require(raster)
require(maptools)

# Raster
r = raster()
r[]= 1

# Spatial points
x = c(-100, 0, 100)
y = c(100, 0, 100)
points = SpatialPoints(data.frame(x,y))

plot(r, xlim=c(xmin(r), xmax(r)), ylim=c(ymin(r), ymax(r)))
plot(points, add=T)

Answer 1

Возможно, вам лучше всего использовать одну из lattice функций для построения пространственных растровых объектов, предоставляемых пакетами raster и rasterVis. Вы обнаружили один из них в vectorplot(), но spplot() или levelplot() лучше соответствуют вашим потребностям в этом случае.

(Метод plot() на основе base graphics для "RasterLayer" объектов просто не позволяет вам легко устанавливать оси с соответствующим соотношением сторон. Для всех, кто заинтересован, я более подробно расскажу, почему это так. так в разделе внизу поста.)

В качестве примера типа сюжета, который levelplot() производит:

require(raster)
require(rasterVis)

## Create a raster and a SpatialPoints object.
r <- raster()
r[] <- 1:ncell(r)
SP <- spsample(Spatial(bbox=bbox(r)), 10, type="random")

## Then plot them    
levelplot(r, col.regions = rev(terrain.colors(255)), cuts=254, margin=FALSE) +
layer(sp.points(SP, col = "red"))

## Or use this, which produces the same plot.
# spplot(r, scales = list(draw=TRUE), 
#        col.regions = rev(terrain.colors(255)), cuts=254) +
# layer(sp.points(SP, col = "red"))

Любой из этих методов может по-прежнему отображать некоторую часть символа, представляющего точки, которые выпадают непосредственно за пределы растрового графика. Если вы хотите избежать этой возможности , вы можете просто поднастроить свой объект SpatialPoints, чтобы удалить любые точки, выходящие за пределы растра. Вот простая функция, которая сделает это за вас:

## A function to test whether points fall within a raster's extent
inExtent <- function(SP_obj, r_obj) {
    crds <- SP_obj@coord
    ext  <- extent(r_obj)
    crds[,1] >= ext@xmin  &  crds[,1] <= ext@xmax &
    crds[,2] >= ext@ymin  &  crds[,2] <= ext@ymax
}
## Remove any points in SP that don't fall within the extent of the raster 'r'
SP <- SP[inExtent(SP, r), ]

---

Дополнительная слабая деталь о том, почему трудно сделать plot(r) производить плотно прилегающие оси

Когда plot вызывается для объекта типа raster, растровые данные (в конечном итоге) выводятся с использованием либо rasterImage(), либо image(). Какой путь будет выбран, зависит от: (a) типа устройства, на котором выполняется печать; и (b) значение аргумента useRaster в исходном вызове plot().

В любом случае область построения настраивается таким образом, чтобы создавать оси, которые заполняют область построения, а не таким образом, который дает им соответствующее соотношение сторон.

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

• useRaster=TRUE

  ## Chain of functions dispatched by `plot(r, useRaster=TRUE)`
  getMethod("plot", c("RasterLayer", "missing"))
  raster:::.plotraster2
  raster:::.rasterImagePlot
  
  ## Call within .rasterImagePlot() that sets up the plotting region
  plot(NA, NA, xlim = e[1:2], ylim = e[3:4], type = "n",
             , xaxs = "i", yaxs = "i", asp = asp, ...)
  
  ## Example showing why the above call produces the 'wrong' y-axis limits
  plot(c(-180,180), c(-90,90), 
       xlim = c(-180,180), ylim = c(-90,90), pch = 16,
       asp = 1,
       main = "plot(r, useRaster=TRUE) -> \nincorrect y-axis limits")
  

• useRaster=FALSE

  ## Chain of functions dispatched by `plot(r, useRaster=FALSE)`
  getMethod("plot", c("RasterLayer", "missing"))
  raster:::.plotraster2
  raster:::.imageplot
  image.default
  
  ## Call within image.default() that sets up the plotting region
  plot(NA, NA, xlim = xlim, ylim = ylim, type = "n", xaxs = xaxs, 
       yaxs = yaxs, xlab = xlab, ylab = ylab, ...)
  
  ## Example showing that the above call produces the wrong aspect ratio
  plot(c(-180,180), c(-90,90), 
       xlim = c(-180,180), ylim = c(-90,90), pch = 16,
       main = "plot(r,useRaster=FALSE) -> \nincorrect aspect ratio")
  

Answer 2

Так я решил эту проблему

require(raster)
r = raster()

# default for raster is 180 row by 360 cols =  64800 cells
# fill with some values to make more interesting
r[]= runif(64800, 1, 1000)

# Set margin for text
par(mar=c(2, 6, 6, 2))

# Set some controls for the raster cell colours and legend

MyBrks<-c(0,1,4,16,64,256,1E20)                                  
MyLbls<-c("<1","<4","<16","<64","<256","<Max")       
MyClrs<-c("blue","cyan","yellow","pink","purple","red")  

# Plot raster without axes or box or legend
# Note xlim and ylim don't seem do much unless you want to trim x and y
plot(r, 
     col=MyClrs,
     axes=FALSE,
     box=FALSE,
     legend=FALSE
     )

# Set up the ranges and intervals for axes - you can get the min max
# using xmin(r) and ymax(r) and so on if you like
MyXFrm <- -180
MyXTo <- 180
MyXStp <- 60
MyYFrm <- -90
MyYTo <- 90
MyYStp <- 30

# Plot the axes
axis(1,tick=TRUE,pos=ymin(r),las=1,at=seq(MyXFrm,MyXTo ,MyXStp )) 
axis(2,tick=TRUE,pos=xmin(r),las=1,at=seq(MyYFrm ,MyYTo ,MyYStp ))

# Plot the legend use xpd to plot the legend outside the plot region
par(xpd=TRUE)
legend(MyXTo ,MyYTo , 
       legend=MyLbls[1:6],                    
       col= MyClrs,             
       fill=Clrs[1:6],                            
       bg=rgb(0,0,0,0.85),           
       cex=0.9,
       text.col="white",
       text.font=2,                      
       border=NA                 
      ) 

# Add some axis labels and a title
text(-220,0,"Y",font=2)
text(0,-130,"X",font=2) 
text(0,120,"My Raster",font=4,cex=1.5) 

Answer 3

Чувак, я был озадачен и закончил тем, что просто переключил цвет переднего плана на сюжет.Затем вы можете воспользоваться тем, что метод растрового графика вызывает fields:::image.plot, что позволяет просто построить легенду (во второй раз, на этот раз с чернилами!)Это не элегантно, но сработало в этом случае:

    par("fg" = NA)
    plot(r, xlim = c(xmin(r), xmax(r)), ylim = c(ymin(r), ymax(r)), axes = FALSE)
    par(new = TRUE,"fg" = "black")
    plot(r, xlim = c(xmin(r), xmax(r)), ylim = c(ymin(r), ymax(r)), axes = FALSE, legend.only = TRUE)
    axis(1, pos = -90, xpd = TRUE)
    rect(-180,-90,180,90,xpd = TRUE)
    ticks <- (ymin(r):ymax(r))[(ymin(r):ymax(r)) %% 20 == 0]
    segments(xmin(r),ticks,xmin(r)-5,ticks, xpd = TRUE)
    text(xmin(r),ticks,ticks,xpd=TRUE,pos=2)
        title("sorry, this could probably be done in some more elegant way")

Answer 4

Я думаю, что лучшим (или самым простым) решением является использование image():

library(raster)

# Raster
r = raster()
r[]= rnorm(ncell(r))

# Spatial points
x = c(-100, 0, 100)
y = c(100, 0, 100)
points = SpatialPoints(data.frame(x,y))

# plot
image(r)
plot(points, add=T, pch=16, cex=2)