Tag Info

Hot answers tagged

5

You can use the reclassify function in the raster package to reclassify the DEM. The general idea is to generate a reclass matrix which provides the instructions on how to reclassify the continuous DEM elevation values. require(rgdal) require(raster) # Read DEM and convert to raster layer object dem = raster("C:/temp/dem.tif") # Generate a reclass ...


5

In R, the package CircStats is old and of rather limited scope and has been replaced by the more complete Circular package. There are tutorials and a book, Circular Statistics with R (2013, A. Pewsey, M. Neuhäuser, and G. D. Ruxton, Oxford University Press, 208 pp.) which explains how to use it (The R scripts can be downloaded from the resources site of the ...


4

Thanks for the clear question and reproducible example. Your understanding is correct, and this boils down to a bug in rgeos::over, which was fixed a month ago but has not made it into a CRAN release yet. The following is a work-around if you're only interested in the number of intersections: world.map$val = sapply(over(geometry(world.map), pointbuff.spdf, ...


3

In a nutshell, the problem lies in a mismatch between data behavior and some (strong) assumptions you are implicitly making. Diagnosis The strongest of these is that the data are one realization of a second-order stationary process. They clearly are not, as you can tell by comparing the region near (450000, 5075000) in the upper "neck" (which I will call ...


2

This looks like a multinomial regression problem to me. I.e, a logistic regression with more than two choices. Say, choice ~ income, age, zip code, distance to next bus stop, ..... and choice being one of bus / bike / car / walk Maybe the following posts are helpful: http://www.jameskeirstead.ca/blog/how-to-multinomial-regression-models-in-r/ ...


2

As handy as R is for so many tasks, it is important to remember that 1) R is not a GIS and 2) quality mapping is downright difficult compared to creating maps with QGIS or ArcGIS. The following example borrows heavily from two R-bloggers blogs (blog 1 and blog 2). Here, I simply mapped a polygon shapefile using Google Satellite Imagery as a basemap. ...


2

This is an alternative using the function cut to assign elevation values into classes of elevation, making possible to discretize colors in the map (color per class of elevation). #Generate reproducible example library(raster) f = system.file("external/test.grd", package="raster") #path to raster file DEM = raster(f) #import raster file DEM = ...


2

The y argument to extract must be 2D, i.e. a matrix or a Spatial* object (or whatever). Your y query is just an atomic numeric vector, and so is intepreted as two cell numbers (i.e. indexes), which is why you get two missing values since they are way out of bounds (cue debate about whether that kind of out of bounds should trigger a different kind of ...


2

you can convert your aspects into the sine and cosine, compute the mean of the sine's and the mean of the cosine's, then turn it back to aspect using atan2(sine,cosine). For more details, see Wikipedia


2

If you happen to work in R, you can use the vec2dtransf package (which is mine :) ). You would simply need to load your Shapefiles into R using rgdal and define your affine transformation to apply it on the data. After such process, you can export your data to a transformed Shapefile also via rgdal. In vec2dtransf, affine transformations can be defined from ...


2

potentially a multi-part question - 1) plotting grids with legends, 2) including shape files on grid, and 3) animate output images. each with multiple opportunities to accomplish the task. here's a quick run-down of at least 2-methods: using gdal, one should be able to read in the raster - perhaps something like (in a loop to get all rasters). raster = ...


1

In R: library(raster) library(animation) files <- list.files("path/to/asc", pattern = "asc$") saveHTML({ for (i in seq_along(files)) { r <- raster(files[i]) r <- plot(r) ## include additions like counties here } }) The animation package has other options for different output formats rather than HTML. The raster package has ...


1

Most textbooks suggest using atan2(Sigma(sin(x)), Sigma(cos(x))), however this is not always the right thing to do. For example, the average of 0, 0 and 90 degrees is atan( (sin(0)+sin(0)+sin(90)) / (cos(0)+cos(0)+cos(90)) ) = atan(1/2)= 26.56 deg, and not 30 deg as one may expect. Take a look at my article on CodeProject "Circular Values Math and ...


1

Yes, using a lapply, e.g.: x = c(0,1,1,0,0) y = c(0,0,1,1,0) bbox_matrix_sp = cbind(rep(x,13),rep(y,13)) require(sp) sp_re_alle = SpatialPolygons(lapply(1:13, function(x) Polygons(list(Polygon(bbox_matrix_sp[((x-1)*5+1):(x*5),])), paste0("reh",x)))) in case you are converting grids into polygons, there are direct conversion methods in sp, as in ...


1

All of the previous recommendations are solid approaches for species distribution modeling. However, an appropriate modeling approach really depends on your question and what you want. Do you want to draw inference from the model? Do you want a probabilistic estimate? Do you want to incorporate spatial process into the estimates? Do you want intensity and ...


1

Random Forests (RF) is a very powerful ensemble learning approach for regression (and classification) that is often used with spatial data. RF is well suited for spatial data because there are no parametric assumptions, which means that you can use binary (open/closed area, burned areas), categorical (soil type), continuous (distance from roads and rivers, ...


1

As per my comment in that question, it's to do with the spTransform. Simply remove the offending line: world <- spTransform(world, CRS("+proj=robin"))


1

Others may be able to specify a way to place the UKN point below the legend as you've specifically requested, but in the meantime you can plot the the point within the bounds of the map simply by specifying a coordinate: data$X1[data$id == "UNK"] <- -150 data$X2[data$id == "UNK"] <- 0 Then plot your map as normal. I appreciate this isn't exactly ...


1

Just a tiny error, as far as I can see. Your substrings are incorrect. This can be seen by comparing the result from a 'which(df$bits="0000100001000100")' with a number of observed unique values, which can be seen in ArcGIS when colouring the tif-file by unique values. 00001000 01000100 = 2116, and there are 3891233 of that number in both ArcGIS and R. This ...


1

As @Jeffrey states, readOGR from the rgdal library imports a CRS if there is one embedded in the shapefile. You can check by (example using a shp I've been playing with): proj4string(india) # from the sp package # [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" If this returns NA you can specify a CRS with: proj4string(india) <- ...


1

I have been able to replicate your problem using the provided shapefile (original). I've opened original.shp in R using both maptools and rgdal and have successfully plotted it. I've also been able to open this unedited file in QGIS 2.8 and ArcMap 10.1. In all cases the Japanese characters in the attribute table displayed correctly. However, when I edited ...


1

Tom Hengl told me this: Set the 'check.module.exists = FALSE' and 'warn=FALSE' -> this usually does the trick (http://www.rdocumentation.org/packages/RSAGA/functions/rsaga.geoprocessor). And Alexander Brenning told me that: did you notice the warning message, Warning message: In rsaga.geoprocessor(lib, module, param = list(h = ""), env = env, : This ...



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