http://ift.tt/2iBq0Dp

I am using the focal function from raster package v2.5-8 to get max value in a 3x3 window. I expect the edges of both rows/columns to return as NA, instead the output returned is 9,9,9. Is this correct ?

Example :

library(raster); require(rasterVis)
r <- raster(nrows=3, ncols=3)
r[] <- 1:ncell(r)
plot(r);text(r);
r.class <- focal(r, w=matrix(1,nrow=3,ncol=3), fun=max) 
plot(r.class); text(r.class);

Output:

     [,1] [,2] [,3]
[1,]   NA   NA   NA
[2,]    9    9    9
[3,]   NA   NA   NA

Expected Output:

       [,1] [,2] [,3]
[1,]   NA   NA   NA
[2,]   NA    9   NA
[3,]   NA   NA   NA

This code works for adding a legend title to an active plot:

library(rasterVis)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
levelplot(r, margin=F)
trellis.focus("legend", side="right", clipp.off=TRUE, highlight=FALSE)
grid::grid.text(expression(m^3/m^3), 0.2, 0, hjust=0.3, vjust=1.5)
trellis.unfocus()

But it doesn't work for a plot saved to a variable, e.g.: p1 <- levelplot(r, margin=F).

This feature would be really useful, for example, for adding data units to the plot's legend.

One possible solution, achieved by changing some fundamental levelplot code, is suggested here: http://stackoverflow.com/questions/36874897/r-how-to-add-legend-title-to-levelplot-saved-to-a-variable

The streamplot plot help is not clear about the pc parameter of the droplet configuration : the example (0.5) is given in percentage, but according to the code, it should be given in integer between 0 and 100 (line 99: nDroplets <- ncell(cropField) * droplet$pc/100).

http://ift.tt/2lyfDxD

I have created a terrain map of my region of interest using:

library(dismo)
library(rasterVis)
gm = gmap(extent(c(-155.0181, -39.61941, 41.68144, 85.1355)),zoom=NULL,type="terrain",lonlat=TRUE,scale=2)
gm <- trim(gm)

e6 <- extent(-141, -96, 46.99998, 70.0)
e6pol <- as(e6, 'SpatialPolygons')
centroid <- coordinates(e6pol)

library(Cairo)
Cairo(file="whatever.png", 
      type="png",
      units="in", 
      width=4, 
      height=4, 
      pointsize=12*96/72, 
      dpi=1000)

levelplot(gm,maxpixel=ncell(gm),panel=panel.levelplot.raster,col.regions=gm@legend@colortable,
          interpolate=TRUE,colorkey=F,margin=FALSE,at=0:255,xlab=list(label="Longitude",cex=1),yscale.components = yscale.raster.subticks,
          xscale.components = xscale.raster.subticks,
          ylab=list(label="Latitude",cex=1),scales=list(x=list(cex=1),y=list(cex=1)),
          par.settings=list(axis.line=list(lwd=1.2), strip.border=list(lwd=1.2)),
          xlim=c(-155.0997, -76.0997),ylim=c(39.09998, 80.09997 ))
dev.off()

This gives me the plot below:

How can one chnage the colors such that water bodies are darkblue and land areas are green to gray, for example? Also add a colorkey.

I like gplot function because this allows to plot big rasters without having to show all tiles.
However, I sometimes need to plot layers below my raster object or plot 2 rasters one over the other. As far as I see it, gplot does not allow for this. Hence, I created a modified version of gplot to retrieve data in the tidy way instead of the ggplot layer. This also allows me to clean it a little while removing NA values for instance and to be more flexible on my ggplots (changing names, levels, aesthetics, ...).
I think that this function would have its place in {rasterVis} instead of my own SDMSelect dev package where it is currently.

#' Transform raster as data.frame to be later used with ggplot
#' Modified from rasterVis::gplot
#'
#' @param x A Raster* object
#' @param maxpixels Maximum number of pixels to use
#'
#' @details rasterVis::gplot is nice to plot a raster in a ggplot but
#' if you want to plot different rasters on the same plot, you are stuck.
#' If you want to add other information or transform your raster as a
#' category raster, you can not do it. With `SDMSelect::gplot_data`, you retrieve your
#' raster as a data.frame that can be modified as wanted using `dplyr` and
#' then plot in `ggplot` using `geom_tile`.
#' If Raster has levels, they will be joined to the final tibble.
#'
#' @export

gplot_data <- function(x, maxpixels = 50000)  {
  x <- raster::sampleRegular(x, maxpixels, asRaster = TRUE)
  coords <- raster::xyFromCell(x, seq_len(raster::ncell(x)))
  ## Extract values
  dat <- utils::stack(as.data.frame(raster::getValues(x)))
  names(dat) <- c('value', 'variable')

  dat <- dplyr::as.tbl(data.frame(coords, dat))

  if (!is.null(levels(x))) {
    dat <- dplyr::left_join(dat, levels(x)[[1]],
                            by = c("value" = "ID"))
  }
  dat
}

In my SDM Vignette, I have an example on how to use it then:

# Get partial raster data to plot in ggplot
pred.r.gg <- gplot_data(pred.r)
# Plot
ggplot() +
  geom_tile(
    data = dplyr::filter(pred.r.gg, variable == "resp.fit"), 
            aes(x = x, y = y, fill = value)) +
  scale_fill_gradient("Probability", low = 'yellow', high = 'blue') +
  coord_equal()

Tell me if you think, this could fit into {rasterVis} and if I can propose a PR.

http://ift.tt/2jF42Qu

I want to overlay images made by levelplot(). Overlaying a list of 16 trellis objects on each panel of the raster is a nightmare. I can overlay the maps individually using grid resources but I prefer the beauty of levelplot(). Let me make a reproducible code of how the maps are realised.

lapply(c("raster", "rasterVis", "RColorBrewer"), require, character.only = TRUE)  
# create a rasterbrick
my.brick <- vector("list", 16)
for(i in 1:16){
  my.brick[[i]] <- raster(xmn = 30, xmx = 42, ymn = -6, ymx = 6)
  values(my.brick[[i]]) <- runif(ncell(my.brick[[i]]), 1, 250)
}
my.brick <- stack(my.brick) # the rasterlayers

cols <- c("white", brewer.pal(9, "Reds")) # my colours
my.at <- seq(0, 250, 25)

p1 <- levelplot(my.brick, col.regions = cols, at = my.at) # plot the raster
coods <- vector("list", length = 16) # object to hold my point data
my.plots <- coods  # An empty list to hold trellis objects

for(i in 1:16){ # A loop to create the spatial points
  coods[[i]] <- data.frame(lat = runif(10, 30, 50), lon = runif(10, 0, 20), 
                           val = runif(10, 1, 250))
  attach(coods[[i]])  ## let the headers become R objects
  coordinates(coods[[i]]) <- ~ lon + lat # convert to a sp object.
  my.plots[[i]] <- levelplot(val ~ lon + lat, col.regions = cols, at = my.at, 
                             panel = panel.levelplot.points, cex = 1.3) + 
    layer(sp.points(coods[[i]], pch = 21, bg = "white", col = "black", 
                    lwd = 2, cex = 1.5)) # Plots. The layer is meant enhance 
                                         # the width of the symbol outline
}

How do I overlay themy.plots which is a list to the raster plot p1? There seems to be a few answers about this on stack exchange but I didn't find them sufficient.

http://stackoverflow.com/a/33147520/964866

http://ift.tt/2mW0vyE

I'm an undergraduate student who was recently given an assignment within R which I have minimal capability of approaching. Its mapping population density on an empty schematic of the country. I have no clue how to properly link the documents to make the below program work. Its very simple for most of you I assume. Pls help. The Poppath, ken08 and Ken adm are the documents I am assuming but have no clue how to make them actually link properly.

rm(list=ls())

library(sp)
library(rgdal)
library(raster)
library(rasterVis)
library(latticeExtra)

PopPath <- "/Users/YUE/Documents/Teaching/Eco 329 17/R/example1"


kenya <- raster( paste(PopPath, "ken08povmpi-uncert.tif", sep="/") )
kenya
MyColors <- colorRampPalette( colors=c("red", "yellow", "green"), 
                              space="Lab")
MyTheme <- rasterTheme(region=MyColors(30))
p0 <- levelplot(kenya, par.settings=MyTheme, 
                contour=TRUE, margin=FALSE)
print(p0)


KEN_level2 <- readRDS( file=paste(PopPath, "KEN_adm2.rds", sep="/") )

p1 <- p0 + layer(sp.polygons(KEN_level2, lwd=0.4, col="red"))
print(p1)

http://ift.tt/2nTMbn3

I've been researching this for a while now but haven't come across any solution that fit my needs or that I can transform sufficiently to work in my case:

I have a large car sharing data set for multiple cities in which I have the charging demand per location (e.g. row = carID, 55.63405, 12.58818, charging demand). I now would like to split the area over the city (example above is Copenhagen) up into a hexagonal grid and tag every parking location with an ID (e.g. row = carID, 55.63405, 12.58818, charging demand, cell ABC) so I know which hexagonal cell it belongs to.

So my question is twofold: (1) how can I create such a honeycomb grid with a side length of 124 meters (about 40000 sqm which is equivalent to 200x200 meters but nicer in hexagonal) in this area:

my_area <- structure(list(longitude = c(12.09980, 12.09980, 12.67843, 12.67843), 
                       latitude = c(55.55886, 55.78540, 55.55886, 55.78540)), 
                  .Names = c("longitude", "latitude"), 
                  class = "data.frame", row.names = c(NA, -4L))

(2) How can I then associate all of my points on the map with a certain grid cell?

I'm really lost at this point, I tried to use tons of packages like rgdal, hexbin, sp, raster, rgeos, rasterVis, dggridR, ... but none of them got me to where I want to go. Help is much appreciated!

Parking data example:

                  id latitude longitude           timestamp charging_demand
1: WBY1Z210X0V307780 55.68387  12.60167 2016-07-30 12:35:07              22
2: WBY1Z210X0V307780 55.63405  12.58818 2016-07-30 16:35:07              27
3: WBY1Z210X0V307780 55.68401  12.49015 2016-08-02 16:00:08              44
4: WBY1Z210X0V307780 55.68694  12.49146 2016-08-03 13:40:07               1
5: WBY1Z210X0V307780 55.68564  12.48824 2016-08-03 14:00:07              66
6: WBY1Z210X0V307780 55.66065  12.60569 2016-08-04 16:19:15              74

Hi Oscar,

Sorry to bother you again.

Please consider the following example. The key.arrow for vectorplot appears on the upper left corner when used alone. But when vectorplot and levelplot are used together, the key.arrow doesn't appear. Do yo think there is an easy solution to that problem?

Best
Pascal

library(rasterVis)
library(gridExtra)

proj <- CRS('+proj=longlat +datum=WGS84')

df <- expand.grid(x=seq(-2, 2, .01), y=seq(-2, 2, .01))
df$z <- with(df, (3*x^2 + y)*exp(-x^2-y^2))
r1 <- rasterFromXYZ(df, crs=proj)

p1 <- 
  vectorplot(r1, , region = FALSE, key.arrow = list(size = 5, label = 'm/s'))

p2 <- 
  levelplot(r1, margin = FALSE) +
  vectorplot(r1, , region = FALSE, key.arrow = list(size = 5, label = 'm/s'))

grid.arrange(p1, p2, nrow = 1)

r <- raster(system.file("external/test.grd", package="raster"))
s <- stack(r, r*2)
names(s) <- c('meuse', 'meuse x 2')

library(ggplot2)

theme_set(theme_bw())
gplot(s) + geom_tile(aes(fill = value)) +
facet_wrap(~ variable) +
scale_fill_gradient(low = 'white', high = 'blue') +
coord_equal() +
theme(
plot.background = element_blank()
,panel.grid.major = element_blank()
,panel.grid.minor = element_blank()
,panel.border = element_blank()
,panel.background = element_blank())

I am trying to create a hovmoller plot for temperature anomaly data. I have set my color scale as reverse RdBu running from -2.5 to 2.5 with 10 steps. It seems to work fine, but I am concerned about what happens if the anomaly is outside the (-2.5, 2.5) range. Let's say I wind up with a value of 3.0. Does it show up as white? Or deep red like 2.5 does? Here's the code:

hovmoller(tmpplt, dirXY=y,
at = do.breaks(c(-2.5,2.5),10),
contour = F, interpolate = F,
par.settings=RdBuTheme(region=rev(brewer.pal(9,'RdBu'))),
main="HadSST3 Anomalies 1850-2018 (1961-1990 base period)")

Any help is appreciated. Thanks.

http://developer.r-project.org/blosxom.cgi/R-devel/NEWS/2015/06/29#n2015-06-29

In previous version of raster, sampleRegular didn't provide xy values when maxpixels > ncell. This has already been fixed.
histogram and densityplot have an if-else loop that isn't needed anymore.

I'm building rasterVis 0.37 from source on OS X 10.10.4 for use with R 3.1.3. I'm getting an error during the source build validation step. I believe this worked correctly in previous versions of rasterVis.

From rasterVis_0.37.tar.gz built with an automated script:

* installing *source* package ‘rasterVis’ ...
** package ‘rasterVis’ successfully unpacked and MD5 sums checked
** R
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded
Error in unloadNamespace(package) : 
  namespace ‘lattice’ is imported by ‘sp’ so cannot be unloaded
Error in library(pkg, character.only = TRUE, logical.return = TRUE, lib.loc = lib.loc) : 
  “lattice” version 0.20.33 cannot be unloaded.
Error: loading failed

Some comments by smilemichael about axis.margin:

1. label allows only one decimal number for example 0,1 is fine but 0.01 will show 0 .
2. scales.margin works fine but label of axis.margin is incorrect. Color ramp shows most of values are less than 0.1 and mean value supposed to be less than 0.1 but label of axis.margin was forced by 0.1. That is fine if axis is extended to the position of 0.1.

Can you give a simple example how to plot a streamline for a given uv components of a wind field? I tried all of the examples on the net but I couldn't manage to plot a meaningful streamlines of uv wind field. I get only points on a dark background. For a simple start, assume, we have a 2 matrices of u and v fields. then what?

Dear Oscar,

I have found an issue with the legend and the colours of the raster dXY vector plot.
When I use vectorplot(isField = dXY), the function creates areas with negative magnitudes of the vectors. (The colours indicate the magnitude of the vector, not the direction, right?)
It would seem the most logical to me that the arrows indicate the direction and the colours the magnitude of the vector. As a vector cannot have a negative magnitude, the negative values seem wrong to me. A vector itsself cannot have a negative value, only its components
I think if you would plot the absolute values of the now calculated magnitudes in the plot, this issue will be fixed.

It all seems like a small error in the coding to me, so if I come across as trying to school you; that is not my intention. It would by the way be nice to keep the pastel shade for the low velocities.

On a sidenote:
Perhaps it would be a nice addition to the vector plot to include the physical quantity and the unit above the legend? I am thinking of the placement of short string like "v [m/s]" that can be defined by the user in the function, which is then placed above the legend by the function.

Kind regards,
Joeri

Already implemented in plot3d.

A simple solution was proposed here

library(rasterVis) 
f <- system.file("external/test.grd", package="raster") 
r <- round(raster(f)) 

## assign a color table 
## if the original raster (on disk) had a color table it should already be in this slot 
r@legend@colortable=rainbow(cellStats(r,max)) 
str(r) 

## use the color table instead. Could write a simple wrapper function for this... 
## this might not work for all color tables.... 
cols=r@legend@colortable 
levelplot(r,col.regions=cols,at=0:length(cols)) 

Would it be possible to add a superscript on the grid.text as the title of the scale?

For example:
grid::grid.text('Soil bulk density (kg m^3)'.......
Where the 3 should be over the m

Passing levels to predict randomForest with a raster http://ift.tt/1Gqc6Uo I try to predict forest stand types (n=43 classes coding landform, geology, water- and nutrient supply etc.) on a large 1027206000 cell raster by randomForests. Among many DEM derived parameters that I am using as covariates, I also have 2 rasters with ID numbers of a geological map and a soil map. Many categorical mapping units go with those IDs. I train the models with a dataframe and attach the categorical mapping units by "merge" to it. So far everything is ok. The model does what it should do and I can predict some test data held in a dataframe too. But now I intent to make some maps of the predictions. But when running the model with a rasterstack or brick gives only rasters with all NA's. My impression is, that I do something wrong in passing the factor levels to the rasterstack/rasterbrick. Her is some code that reproduces the problem. library(raster) library(rasterVis) library(randomForest) # make a raster set.seed(0) r <- raster(nrow=10, ncol=10) r[] <- runif(ncell(r)) * 10 is.factor(r) r <- round(r) # make faktor f <- as.factor(r) is.factor(f) # get some none-sense levels x <- levels(f)[[1]] x$code <- paste("A",letters[10:20]) x$code2 <- paste("B",letters[10:20]) x$code3 <- letters[10:20] levels(f) <- x f<-deratify(f) # make a brick levels(f) set.seed(2) # get some none-sense dataframe xx<-data.frame(code=sample(rep(paste("A",letters[10:20]),10)), code2=sample(rep(paste("B",letters[10:20]),10)), code3=rep( letters[10:20],10), y=as.factor(sample(rep(paste(rep(1:5)),22)))) # fit and predict a random forest with it ranfor<-randomForest(y~.,data=xx,ntree=100) predict(ranfor) # try to predict with a raster names(f)<-c("code","code2","code3") a<-predict(object=f,ranfor,na.omit=T,factors=list(code=levels(xx$code), code2=levels(xx$code2), code3=levels(xx$code3))) plot(a) # gives an empty raster # convert the raster to a dataframe and predict again x<-as.data.frame(f) names(x)<-c("code","code2","code3") aa<-predict(ranfor,x) plot(aa) # works just fine! Any suggestions? Thank you! R version 3.1.2 (2014-10-31) Platform: x86_64-w64-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] randomForest_4.6-10 rasterVis_0.35 latticeExtra_0.6-26 RColorBrewer_1.1-2 lattice_0.20-33 raster_2.4-15 sp_1.1-1 loaded via a namespace (and not attached): [1] grid_3.1.2 hexbin_1.27.0 Rcpp_0.11.6 rgdal_1.0-4 tools_3.1.2 zoo_1.7-12 via Active questions tagged r - Stack Overflow http://ift.tt/1p1c5QJ July 17, 2015 at 03:32PM

View on Trello

http://stackoverflow.com/questions/32252946/one-side-margin-in-levelplot

Dear Oscar,

Considering these data:

u <- v <- raster(xmn=0, xmx=2, ymn=0, ymx=2, ncol=1e3, nrow=1e3)
x <- init(u, v='x')
y <- init(u, v='y')
u <- y * cos(x)
v <- y * sin(x) 
field <- stack(u, v)
names(field) <- c('u', 'v')

Is it intentional for this to work:

vectorplot(field, isField='dXY', narrows=300, uLayer=1, vLayer=2, region=FALSE, col='blue')

And this one to fail? (region is turned from FALSE to TRUE):

vectorplot(field, isField='dXY', narrows=300, uLayer=1, vLayer=2, region=TRUE, col='blue')
Error: $ operator is invalid for atomic vectors

Pascal

R version 3.2.2
rasterVis version 0.37
lattice version 0.20-33

A warning appears while using contourplot function.

library(rasterVis)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
contourplot(r, labels=list(cex=0.4), cuts=12)
Warning in .local(x, data, ...) :
  FUN.margin is deprecated. Use margin as a list instead.
Warning in .local(x, data, ...) :
  axis.margin is deprecated. Use margin as a list instead.

Hi,

Any idea why the following doesn't work? It seems like the colors and key should follow the levels specified in rat$landcover.

r <- raster(nrow=10, ncol=10)
r[] = 1
r[51:100] = 3
r[3:6, 1:5] = 5
r <- ratify(r)

rat <- levels(r)[[1]]
rat$landcover <- factor(c('Pine', 'Oak', 'Meadow'), levels = c('Oak', 'Meadow', 'Pine'))
rat$class <- c('A1', 'B2', 'C3')
levels(r) <- rat

levelplot(r, col.regions=c('palegreen', 'midnightblue', 'indianred1'))

Thanks!

Hi Oscar-
I think the axes and labels should be col="black" or there should be an option to set these colors. Also, an option to change the label text size would be great.
I implemented this in a fork along with some changes to deal with the polygon grob baseline (I opened another issue for this).
Here is a plot made with the edited levelplot function. The call was:

levelplot(Shipping_lognorm, par.settings=myTheme,
                         main=list("Shipping", vjust=1),
                         margin=TRUE, FUN.margin=mean,
                        scales.margin=list(x=c(0,0.8), y=c(0,0.8)),
                        axis.margin=TRUE, cuts=50,
                        panel=function(...){
                           panel.levelplot(...)
                          sp.polygons(land, col=cont.border, fill=cont.fill)
                        })

Dear Oscar,

Is it possible to add/modify "names.attr" option for the "vectorplot" function? It seems to not work with multilayer RasterStack containing horizontal and vertical components:

Error in .local(x, data, ...) :
Length of names.attr should match number of layers.

Regards,
Pascal

http://ift.tt/2gHLUTD

I have a raster that I have reprojected from Cylindrical Equal-Area (Lambert) Central Meridian: -160. Datum: WGS 1984 (+proj=cea +lon_0=Central Meridian +lat_ts=Standard Parallel +x_0=False Easting +y_0=False Northing +ellps=WGS84) to latitude and longitude. The original raster looks like this:

library(rasterVis)
levelplot(r)

Original raster without missing points

r
#class       : RasterLayer 
#dimensions  : 64, 200, 12800  (nrow, ncol, ncell)
#resolution  : 2e+05, 2e+05  (x, y)
#extent      : -20037507, 19962493, -6363885, 6436115  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=cea +lon_0=-160 +lat_ts=Standard Parallel +x_0=False Easting +y_0=False Northing +ellps=WGS84 
#data source : N:\My documents\Data\Exposure\Reefs at risk\Global_Threats\Acidification\arag_380\w001001.adf
#names       : w001001
#values      : 1.025163, 4.11939  (min, max)

I have been able to use the projectRaster function from the Raster package by cutting the extent of the y-axis of the raster by 0.95 to the North and South. I was having a few problems reprojecting the raster without cutting the extent (see here: http://ift.tt/2gHAK0Y).

#Cut y-axis values because projectRaster failed using full extent
extent(r) <- c(xmin= -20037507, xmax= 19962493, ymin= 0.95*(-6363885), ymax= 0.95*(6436115))
Define the new Proj.4 spatial reference
sr <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
Project Raster
projected_r <- projectRaster(r, crs = sr, method = 'bilinear')

The problem I am having is that there are missing data in the new, reprojected raster.

levelplot(projected_r)

reprojected raster with missing points

The missing coordinates are not located in the 5% North and Souththat I cut from the extent so I am not sure why these data are missing? Any help would be appreciated!

http://r-sig-geo.2731867.n2.nabble.com/rasterVis-levelplot-Use-colortable-in-raster-with-attribute-table-tp7585716p7586573.html

gplot is great! Thanks for creating. I'm having one weird issue when aggregating a raster to a lower resolution prior to plotting it with gplot (rather than relying on the maxpixels arguement). This seems to be something deeper within ggplot, but it's related to the approach used by maxpixels, so there might be some way gplot can fix it.

I've been using gplot very successfully to map my data. Here is an example plot, where the supplied raster to gplot has a resolution of ~1e10. That's way too many pixels for my eye to see - and also crashes R if I try to plot it anyway.

midwest.rast <- raster("Midwest_PI_Raster.tif")

midwest.rast

class       : RasterLayer 
dimensions  : 75891, 127386, 9667450926  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 18855, 1292715, 1453065, 2211975  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : /Users/Midwest_PI_Raster.tif 
names       : Midwest_PI_Raster 
values      : 2, 18  (min, max)


counties <- readOGR("cb_2016_us_county_500k/cb_2016_us_county_500k.shp")
midwest_counties <- counties[counties$STATEFP %in% c(17, 39, 18, 29),]

gplot(midwest.rast, maxpixels = 1e+7) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis(option = "magma", na.value = "transparent") +
  theme_void() +
  labs(fill = NULL) +
  coord_equal() +
  guides(fill = guide_colourbar(title.hjust = 0.5, nbin = 100, label.position = "left")) +
  theme(legend.position = c(0.8, 0.25)) +
  geom_polygon(data = midwest_counties,
               aes(x = long, y = lat, group = group),
               color = "black",
               fill = "transparent",
               size = 0.5)

However, I realized that I could speed up my prior calculations by using aggregate to reduce the resolution of my raster first. I've spared the calculations here, but just simply aggregate the original raster using fact=31 to reduce the resolution to ~1e7. This is the same resolution that I was having gplot go to anyway.

midwest.rast.agg <- aggregate(midwest.rast, fact = 31)

midwest.rast.agg

class       : RasterLayer 
dimensions  : 2449, 4110, 10065390  (nrow, ncol, ncell)
resolution  : 310, 310  (x, y)
extent      : 18855, 1292955, 1452785, 2211975  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : /private/var/folders/4p/mpc8rykn1svfp0dd6639lkw40000gn/T/RtmpNE6q5T/raster/r_tmp_2017-07-15_100058_626_42210.grd 
names       : Midwest_PI_Raster 
values      : 2, 18  (min, max)

gplot(midwest.rast.agg, maxpixels = 1e+7) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis(option = "magma", na.value = "transparent") +
  theme_void() +
  labs(fill = NULL) +
  coord_equal() +
  guides(fill = guide_colourbar(title.hjust = 0.5, nbin = 100, label.position = "left")) +
  theme(legend.position = c(0.8, 0.25)) +
  geom_polygon(data = midwest_counties,
               aes(x = long, y = lat, group = group),
               color = "black",
               fill = "transparent",
               size = 0.5)

The resulting maps looks identical EXCEPT that the alpha value for the fill is clearly less than 1. This seems totally random to me. I've even tried to explicitly set the alpha in scale_fill_viridis to 1, but that's the default value anyway and doesn't do anything.

Furthermore, if you zoom in, you can see what appears to be a phantom grid from ggplot. Not sure why that's there either, since I used theme_void().

To be fair, even if I manually convert this aggregated raster to a SPDF and then supply it directly to ggplot and geom_tile (no gplot), the same issue happens. So I guess it because a question of (1) what is gplot doing with maxpixels that makes the saturation correct and (2) is there an update to gplot that could prevent this weird desaturation when supply an aggregated raster?

Presently, names.attr gets coerced to character before being passed to strip.custom's factor.levels. This means that expressions (e.g. to permit subscript/superscript/symbols) are broken. Allowing these would be nice. Related SO post here.

r <- raster(nrow=10, ncol=10)
r[] = 1
r[20:30] = 10
r[51:80] = 21
r[3:6, 1:5] = 25

r <- ratify(r)
rat <- levels(r)[[1]]

rat$labs <- letters[seq_len(nrow(rat))]
levels(r) <- rat

levelplot(r, col.regions=c('black', 'red', 'yellow', 'green'))

The "c" category (yellow) is not shown

Hello Oscar,

The vectorplot function seems to be unable to display rasterstack with at least one layer completely filled with NA.

Reproducible code:

u <- v <- raster(xmn=0, xmx=2, ymn=0, ymx=2, ncol=1e3, nrow=1e3)
x <- init(u, v='x')
y <- init(u, v='y')
u <- y * cos(x)
v <- y * sin(x)
u1 <- cos(y) * cos(x) * NA
v1 <- cos(y) * sin(x) * NA
u2 <- sin(y) * sin(x)
v2 <- sin(y) * cos(x)
field <- stack(u, u1, u2, v, v1, v2)
names(field) <- c('u', 'u1', 'u2', 'v', 'v1', 'v2')
vectorplot(field, isField='dXY', narrows=300, uLayer=1:3, vLayer=6:4)

Best,
Pascal

Using rasterVis version 0.36.
When comparing the original .Rd file and the compiled help page of hovmoller function, two differences lead to error when trying to run the examples:

1. listFich <- dir(pattern='\\.nc') became listFich <- dir(pattern='\.nc') in ?hovmoller
2. month <- function(x)format(x, '%m') became month <- function(x)format(x, ' in ?hovmoller

What could explain these differences?

Hi, Thanks for rasterVis! I'm having a small issue with rasterVis_0.4.

The example below produces vectors in the 2nd and 4th quadrant that aren't in the direction I expect (toward the origin).

## simple example
library(rasterVis)
r <- raster(nrows=2, ncols=2, xmn=0, xmx=5, ymn=0, ymx=5)
set.seed(111)
values(r) <- runif(4, 0, 5)
dX <- r
values(dX) <- c(1.5, -1.5, 1.25, -2.7)
dY <- r
values(dY) <- c(-1.5, -1.5, 1.25, 2.7)
comp <- stack(list(dx=dX, dy=dY))
mag <- r
values(mag) <- sqrt(values(dX)^2 + values(dY)^2)
ang <- r
values(ang) <- atan2(x=values(dX), y=values(dY))
vf <- stack(list(mag=mag, slope=ang))
# plot(comp)
# plot(vf)
vectorplot(vf, isField=TRUE, units="radians", lwd=5, col='orange',
           par.settings=viridisTheme(), title="vf of mean net movement")

Please advice how best to remove the title frame and its content as it's taking too much space for my particular small multiples. I'd like to do letter labeling in-frame.

Dear Oscar,

I have some trouble with the streamplot function. I want to plot the streamlines of a vector field I exported from Hydrus. I managed to get the vectors into a Rasterstack (by adding a nonsense projection, as I am dealing with a vertical crossection).

The thing is: When I try to use the vectorplot function with isField = 'dXY', there are no issues whatsoever: I get the same plot as Hydrus gave me; but when I try to use the streamplot function with 'dXY', I get the following error:
Error in if (isField) callNextMethod(object, layers, droplet, streamlet, :
argument is not interpretable as logical

I interpret this as that, in the streamplot function, the isField accepts either a TRUE or FALSE, but no
'dXY' . What am I doing wrong?

I hope that we can get this to work, as rasterVis currently is the only package for R that can generate streamlines from a vector field and the pictures that I got so far looked very good.

Kind regards,
Joeri

Plotting a logical Raster layer with levelplot works as expected (coercion to numeric), but this fails for a logical RasterStack. E.g.:

set.seed(1)
r <- raster(matrix(runif(9), 3))

levelplot(r > 0.5)
## works fine

plot(stack(r, r) > 0.5)
## works fine

levelplot(stack(r, r) > 0.5)
## Error in FUN(X[[i]], ...) : 
##   only defined on a data frame with all numeric variables 

levelplot((stack(r, r) > 0.5)[[1]])
## works fine

Hi Oscar,

I noticed a strange behaviour when using levelplot with mislabeled arguments. Accidentally I passed a theme to col.regions instead of par.settings, and if normal R script is used, you just get an incorrect plot and warnings() warning. However, if the same code is embedded in markdown chunk, this causes flooding of text warnings output which completely freezes R session for me. Here is some reproducible markdown example.

`

title: "test"
author: "Author"
date: "February 6, 2018"
output: html_document

library(raster)
library(rasterVis)

r <- raster(volcano)
levelplot(r, col.regions = RdBuTheme)

`

http://ift.tt/2nTMbn3

I've been researching this for a while now but haven't come across any solution that fit my needs or that I can transform sufficiently to work in my case:

I have a large car sharing data set for multiple cities in which I have the charging demand per location (e.g. row = carID, 55.63405, 12.58818, charging demand). I now would like to split the area over the city (example above is Copenhagen) up into a hexagonal grid and tag every parking location with an ID (e.g. row = carID, 55.63405, 12.58818, charging demand, cell ABC) so I know which hexagonal cell it belongs to.

So my question is twofold: (1) how can I create such a honeycomb grid with a side length of 124 meters (about 40000 sqm which is equivalent to 200x200 meters but nicer in hexagonal) in this area:

my_area <- structure(list(longitude = c(12.09980, 12.09980, 12.67843, 12.67843), 
                       latitude = c(55.55886, 55.78540, 55.55886, 55.78540)), 
                  .Names = c("longitude", "latitude"), 
                  class = "data.frame", row.names = c(NA, -4L))

(2) How can I then associate all of my points on the map with a certain grid cell?

I'm really lost at this point, I tried to use tons of packages like rgdal, hexbin, sp, raster, rgeos, rasterVis, dggridR, ... but none of them got me to where I want to go. Help is much appreciated!

Parking data example:

                  id latitude longitude           timestamp charging_demand
1: WBY1Z210X0V307780 55.68387  12.60167 2016-07-30 12:35:07              22
2: WBY1Z210X0V307780 55.63405  12.58818 2016-07-30 16:35:07              27
3: WBY1Z210X0V307780 55.68401  12.49015 2016-08-02 16:00:08              44
4: WBY1Z210X0V307780 55.68694  12.49146 2016-08-03 13:40:07               1
5: WBY1Z210X0V307780 55.68564  12.48824 2016-08-03 14:00:07              66
6: WBY1Z210X0V307780 55.66065  12.60569 2016-08-04 16:19:15              74

• viridis
• Change default theme?

https://github.com/hypertidy/graticule

https://github.com/rspatial/terra/blob/master/R/graticule.R

http://ift.tt/2i1majl

I am using the focal function from raster package v2.5-8 to get max value in a 3x3 window. I expect the edges of both rows/columns to return as NA, instead the output returned is 9,9,9. Is this expected ?

Example :

library(raster); require(rasterVis)
r <- raster(nrows=3, ncols=3)
r[] <- 1:ncell(r)
plot(r);text(r);
r.class <- focal(r, w=matrix(1,nrow=3,ncol=3), fun=max) 
plot(r.class); text(r.class);

Wrong Output:

     [,1] [,2] [,3]
[1,]   NA   NA   NA
[2,]    9    9    9
[3,]   NA   NA   NA

Correct Output:

       [,1] [,2] [,3]
[1,]   NA   NA   NA
[2,]   NA    9   NA
[3,]   NA   NA   NA

Hi Oscar,

I just noticed a difference using, for example, levelplot on NetCDF read with the raster package. With raster_2.3-12, the coordinates are correctly displayed (N-S, W-E), but not with raster_2.3-24, the most recent version. Do you have an idea to explain this difference?

Best
Pascal

raster_2.3-12:

raster_2.3-24:

Given a large RasterStack, plotting with the default raster plot method can be interrupted with Esc.

When I plot with levelplot, interrupting with Esc halts plotting but the console remains active indefinitely (okay, maybe the length of time can be defined... I lost patience). This happens when plotting to the RStudio gui device, but not when plotting to, e.g., an x11 device (which allows me to interrupt plotting, returning immediately to the console).

Here's an example that reproduces the problem, on my system at least...

library(rasterVis)
s <- stack(replicate(10, raster(matrix(runif(1e6), 1e3))))
levelplot(s)
# now hit Escape after plotting has commenced
# Eventually, kill the RStudio process and cry about your lost workspace

But no problem with the following:

x11()
levelplot(s)
# now hit Escape and get on with what you were doing, no worries

Similarly, no problem if the stack is small:

s <- stack(replicate(10, raster(matrix(runif(100), 10))))
levelplot(s)
# Hit escape really fast to beat the plotting

I'm using RStudio Version 0.99.837 under Windows 8 x64. My sessionInfo() is:

R version 3.2.2 (2015-08-14)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8 x64 (build 9200)

locale:
[1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252    LC_MONETARY=English_Australia.1252
[4] LC_NUMERIC=C                       LC_TIME=English_Australia.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rasterVis_0.38      latticeExtra_0.6-26 RColorBrewer_1.1-2  lattice_0.20-33     raster_2.4-20       sp_1.2-1           

loaded via a namespace (and not attached):
[1] zoo_1.7-12     parallel_3.2.2 tools_3.2.2    hexbin_1.27.1  Rcpp_0.12.2    grid_3.2.2

Dear Oscar,

Is it possible/feasible to derive streamlines, instead of droplets, from the "streamplot" function?

Regards,
Pascal

Hi Oscar,

Is it easy to add a (optional) legend to a vector plot? Typically, an arrow outside the frame giving a wind speed.

Best,
Pascal

Hi Oscar-
In working with the new axis.margin parameter, I came across another issue (I think). when you define scale.margin, it partially adjusts the margin plot. The plot is shifted according to the provided scale. However, because the baseline of the polygon grob is calculated as the minimum of the range of values (rVal[1]), the baseline of the polygon is not dropped to scale.x[1] or scale.y[1].
I forked your code and made some edits to fix this and to change the margin axes to black and make the axes text larger. If you want I can make a pull request, but you will probably come up with a more elegant way to do this.
I think the ideal would be to allow the specification of either both a minimum and maximum x and y in scale.margin or just a minimum, letting the code calculate the maximum. Of course, you could still leave scale.margin=NULL to allow the code to calculate the range and apply to the margin plot. The reason I say this is because in some cases, plotting just the range of values is not really representative, but instead there may be some potential range for the data on which the margins should be plotted. If this is the case, then plotting just min-max is not properly representative. I hope this makes sense.
Cheers,
Cotton

Hi Oscar

Is there a way to use continuous color plates in rastervis instead of the discrete default?