Tag Info

New answers tagged

1

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 ...


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

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 ...


0

Now that I have significantly more R experience under my belt, it's time to answer this question... First was to make sure the geocoded address was converted to a spatialpoint cleaner: AddressDetail <- geocode(match_address) coordinates(AddressDetail) <- ~lon + lat Then the relevant projection was assigned to the spatialpoint wgs84 <- ...


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 <- ...


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 = ...


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 ...


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 ...


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 ...


0

Package raster provides a function calc that executes a function over all raster pixels, returning the raster with results; it can deal with large rasters (out-of-memory). raster also provides cluster/multi-core functionality, look into function beginCluster; it's meant to work with calc.


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, ...


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 ...


0

Here is a suggestion. I create the circles with gBuffer and then reproject them into WGS84 for ggmap. To change the colors of the heat map use scale_fill_gradient(). library(ggmap) library(sp) library(rgdal) library(rgeos) # get the NY coordinates nyc <- geocode("New York") # create spatialPoint object coordinates(nyc) <- ~ lon + lat ...


-1

It really depends on how your files are structured, namely if the have multiple bands. Assuming that you have only one band in each file you can use the following code to make a rasterstack. Using the 'raster' library makes it easier. library(raster) datafiles <- Sys.glob("*.tif") #Or whatever identifies your files resultingStack <- stack() for(i in ...


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

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) { ...


0

One more description from Josh on how to implement it in QGIS Very good and detailed tutorial.


0

I would suggest using the MODIS reprojection tool (MRT - https://lpdaac.usgs.gov/tools/modis_reprojection_tool) to project and convert the data from HDF to TIFF. It's free (just need to create an account with NASA) and you'll be sure you're data is being transformed properly. Then you can work with your TIFF in R.


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) ...


0

A pesky feature of rasters is the multitude of "no data" possibilities, and how software doesn't consistently deal with it. What I expect that GRASS is doing is assigning minus infinity to the missing values. While not the best of options, you could try to go about it the other way around, eg. identify all acceptable values and then assign NA to the ones ...


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 <- ...


0

All models have an implicit (constant) covariate whose name is "1", which may be mentioned explicitly, as in y ~ 1. An example appears near the bottom of p. 12 at http://cran.r-project.org/web/packages/spdep/vignettes/sids.pdf.


0

I think you have to be more specific. Do you want something like this http://evansmurphy.wix.com/evansspatial#!random-forest-sdm/c1iuu?


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, ...


0

I'd split the area into square polygons (size is a challenge), compute count of points per each scenario. Remove ones that have none in both, check correlation between counts.


0

Use the subs function (also possible with reclassify, calc, etc.) I demonstrate the process with a small raster first: raster1 <- raster(ncol=99, nrow=99) # is the raster small enough to use the following functions canProcessInMemory(raster1) # TRUE, so the simple method of assigning values should work setValues(raster1, 100) raster1[10:50, 10:50] <- ...


0

I believe that you are making this much more difficult than it needs to be. You can, in fact, read ESRI format rasters (not filegeodatabases) using rgdal, which also means that the raster package can read ESRI rasters. Since your raster(s) is so large it will continue to plague you when operating in memory. By using the raster package, you can keep the ...


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

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 ...


0

install raster package, then load band by band. Then you can use landsat package. Ex: library(raster) setwd(L8 folder) b1=raster("L8b1.tif") b2=raster("L8b2.tif") image=stack(b1,b2)


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, ...


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, ...


0

Here, see if this works: # This sets the correct dependencies on CRAN setRepositories(ind = c(1,6)) # Then install package as usual. install.packages('rgdal')


0

The problem is the way you are addressing the raster with newraster <- raster1[is.na(raster2)]. Try is.na(raster2) by itself to see what you get! Trying to set values of a raster using vector terminology[] doesn't always work the way one expects. Use values() to set the actual data values of the raster: library(raster) # Create two small rasters and ...


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

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

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 ...


0

I dont have a very concrete answer for you, but here are some considerations that might help. I guess you have two separate problems here First you need to match the unsupervised classification values to the existing values. This also relates on how many classes your unsupervised classification will create. I think that for some classification algorithms ...


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 <- ...


0

It seems now you can use the shift() function in the {raster} package to move SpatialPolygons... http://www.inside-r.org/packages/cran/raster/docs/shift Actually it didn't work for me. So I used the elide() function in the {maptools} package...


0

Change geom_polygon(aes(x=long, y=lat, group=group), fill='grey', size=.2,color='green', data=data.shape, alpha=0) to geom_polygon(aes(long, lat, group=group, fill='grey', size=.2,color='green', data=data.shape, alpha=0)) m = map + geom_polygon(data=data.shape, aes(long,lat,group=group,fill="grey", size=2, color='green',alpha=0)) m


0

To expand on @Chris's answer, use raster::writeRaster to write to any of a large variety of raster formats. Depending on whether your build of GDAL, you should be able to write to these (and possibly more) filetypes (returned by raster::writeFormats()): # name long_name # 1 raster ...


1

If the square can be in meters instead of degrees, you can create a custom CRS in oblique mercator projection: +proj=omerc +lat_0= 21.04247 +lonc=88.38045 +alpha=-70 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs And create a vector grid with that CRS:



Top 50 recent answers are included