Tag Info

Hot answers tagged

6

Excellent answer by @Martin but it seems you have 2 attributes you want to have a colour gradient (Plantation Cover and Forest Harvesting Intensity). If I'm mistaken then I will remove this post. You could set up a Rule-based style, create a new rule, edit the symbol layer from a Simple Fill to a Gradient Fill. Create a filter for your attribute column and ...


3

The native plot window of R does not have zoom capability, but here are some options: resize your plot window only plot a section of the map by specifying xlim and ylim arguments to plot() write the map to a very large pdf, e.g. precede the plot command by pdf("file.pdf", width=20, height=20) and send this to a plotter write to the plot to pdf, load in a ...


3

Here is a simple example to create a SpatialLinesDataFrame, which can be saved as a shapefile with rgdal::writeOGR(): # create example data set set.seed(1) dat <- matrix(stats::rnorm(2000), ncol = 2) ch <- chull(dat) coords <- dat[ch, ] plot(dat, pch=19) lines(coords, col="red") library("sp") library("rgdal") sp_line <- ...


3

I can reproduce your error. A current workaround is to use gpclib (free for non-commercial use) instead of rgeos: # install.packages("gpclib") library("gpclib") gpclibPermit() # [1] TRUE # Warning message: # In gpclibPermit() : # support for gpclib will be withdrawn from maptools at the next major release LA_new <- unionSpatialPolygons(LA, IDs, ...


3

Starting from Robin's solution, here is an alternative that uses less packages and spDists's ability to compute pairwise distances between two point sets: library(sp) set.seed(2014) x <- SpatialPoints(coords = matrix(rnorm(10), ncol = 2)) y <- SpatialPoints(coords = matrix(rnorm(10), ncol = 2)) plot(x, col = "red") points(y, col = "green") snap = ...


3

In QGIS 2.6 you can create a two color color ramp this way: Double click on the layer (or right click > Properties) > Style > Change from "Single Symbol" to "Categroized" > Color Ramp > Random colors> Select Gradient > Choose the colors and safe > Choose the column with the values > Classify > Apply In the "Column" field ...


3

Thiessen polygons are Voronoi diagrams - there is a 'voronoi' package available in the CRAN archives (not the main repository), but the 'deldir' package does the same job. require(deldir) # Create some points x <- c(32.5, 32.1, 33.5, 32.2, 33.0) y <- c(-2.2, -3.3, -2.3, -2.9, -3.0) # Calculate the Delaunay triangulation, then the tiles. z <- ...


2

A simple for loop will suffice. You can use readGDAL in the rgdal package but I would recommend raster in the raster package. You have to be a bit tricky and use strsplit in the assign function to strip off the ".tif" file extension. setwd("C:/rasters") rlist=list.files(getwd(), pattern="tif$", full.names=FALSE) for(i in rlist) { ...


2

Interesting finding, and thanks for the perfectly prepared reproducible example! Both variograms you mention are meant to prepare for linear kriging in the next step, and you may not want that. The robust variogram was (IIRC) deviced for normal data with some pollution (outliers), but not for count data. I would advice to look at model-based geostatistics, ...


2

I see a different output, namely: > m2 <- update( m1, corr = corExp(c(300, 0.7), form = ~ x + y, nugget = T) ) Error in lme.formula(fixed = log(zinc) ~ 1, data = meuse, random = ~1 | : nlminb problem, convergence error code = 1 message = false convergence (8) this is nlme_3.1-119 on R version 3.1.2 (2014-10-31) Platform: x86_64-pc-linux-gnu ...


2

sp provides a shorter form to select features based on spatial intersection, following the OP example: pts[ply,] as of: points(pts[ply,], col = 'red') Behind the scenes this is short for pts[!is.na(over(pts, geometry(ply))),] The thing to note is that there is a geometry method that drops attributes: over changes behaviour if its second argument has ...


2

In sp, SpatialPoints*, SpatialPixels* and SpatialGrid* (with * omitted or replaced by DataFrame) do support more than 2 spatial dimensions, as OP has done, but SpatialPolygons* and SpatialLines* do not. With gstat you can do 3-D block kriging with 3-D blocks (using block = c(10,10,10)), but you cannot do this for non-rectangular blocks, as OP wants. It is ...


2

The spatial data.frame is not correctly formed. This might work: library(rgeos) library(sp) library(rgdal) wa.map <- readOGR("ZillowNeighborhoods-WA.shp", layer="ZillowNeighborhoods-WA") sodo <- wa.map[wa.map$CITY == "Seattle" & wa.map$NAME == "Industrial District", ] # Don't use df as name, it is an R function # Better to set longitudes as ...


2

Well, you have two potential issues: 1) "HIshp" is in a "UTM Zone 4" projection and the raster you defined is in a geographic projection (lat/long). You can use "spTransform" to reproject "Hishp" so that it aligns with "new_ras". 2) You need to define the attribute in "HIshp" that you want to represent the raster values (ie., rasterize(HIshp, new_ras, ...


1

A NaN is different than NA. The NaN often results from a divide by zero error whereas NA is the R value for no data. These values behave in specific ways and it would be good for you to read some R background material to understand the behavior. Two useful operators to be aware of are: is.na() and is.nan(). y=c(0,1,2,3,4,NA) x=c(0,1,2,3,4,NA) (d=x/y) ...


1

I recent explored exactly this. Depending on what you want you could use the density funciton in the spatstat package. However, ppp.density returns an isotropic density, which is technically the intensity process, and there are no options for kernel type. Here is a toy example using a raster to define extent and cell size. You can define the weights from ...


1

Not sure it is what you are looking for, but here is a start with this workaround (based on this answer): library(osmar) osm_data <- get_osm(complete_file(), source = osmsource_file('extract.osm')) hw_ids <- find(osm_data, way(tags(k == "highway"))) hw_ids <- find_down(osm_data, way(hw_ids)) ways <- subset(osm_data, ids = hw_ids) way_ids <- ...


1

First, reproject the vector data to WGS84 (Lat/Long degrees): HI_WGS84 <- spTransform(HIshp, CRS("+proj=longlat +elips=WGS84")) Then rasterize a new result: new_ras <- raster(nrow=1520, ncol=2288) crs(new_ras) <- crs(HI_WGS84) extent(new_ras) <- c(-159.816, -154.668, 18.849, 22.269) HI_popdens <- rasterize(HI_WGS84, new_ras, ...


1

As pointed out by @mdsumner, You can read a SpatialGridDataFrame directly using readGDAL in the rgdal package. You can also easily coerce raster objects to a SpatialGridDataFrame or matrix. If your imagery is in separate files I would recommend reading them in as a stack and then coercing to a SpatialGridDataFrame. require(raster) require(rgdal) # Using ...


1

The ipdw package accounts for barriers during interpolation. It is an application of Inverse Distance Weighting using path distances. You can find the vignette here. You might also consider non-Euclidean kriging using the geoRcb package on github. I have not tried it but it looks promising.


1

Answer to question 1: it is possibly appropriate, and there are alternatives, and some of them have been implemented rather recently in gstat, see this vignette. Q 2: there is a function zerodist in sp; if the collocated data have different time stemps, it should not be a problem if you do proper ST kriging. Q 3: if gstat is informed that coordinates are ...


1

This can easily be solved using the overlay function from the raster package. Objects rst1 and rst2 are replicates of the initial 'volcano' layer, and a random sample of n = 1000 cells in rst2 is set to NA. Afterwards, overlay is applied and the associated function rejects all cells in rst1 that hold a valid value, i.e. different from NA, in rst2. ...


1

I would recommend raster::calc because according to the documentation: For large objects calc will compute values chunk by chunk. raster2<-reclassify(raster2,c(-Inf,Inf,0)) s->stack(raster1,raster2) rs1<-calc(s, sum) You might want to look into this vignette for instructions on processing large rasters.


1

Yes, it is appropriate. Prediction by kriging can theoretically only get better when you bring in more correlated information, and that is what you do when moving from kriging to co-kriging. In practice, the gain can be disappointing, considering the effort it takes. There can also be other reasons to favor co-kriging. An example is when you need the ...



Only top voted, non community-wiki answers of a minimum length are eligible