Tag Info

Hot answers tagged

6

Nice and reproducible question. Personally, I'd expect that the reason for the difference is in the implementations of the bilinear reprojection. You can obviously look into source code for the two approaches, but I'd expect that to be a vast overkill. It appears that the R implementation introduces bigger "errors" / "changes" than the raw GDAL version ...


6

It sounds like you want to subset your data and plotting is secondary. Please keep in mind that it is not always necessary to create a new object. If you are only wanting to plot subsets of the data it is far more efficient to take @fdetsch advice and subset in the call to plot. Here are some examples of subseting and plotting sp class data using the meuse ...


4

You can use the merge function as there is now a default method for sp class objects. Add sp and the sp "meuse" data (SpatialPointsDataFrame) library(sp) data(meuse) coordinates(meuse) <- ~x+y Add "ID" column to "meuse" meuse@data <- data.frame(ID=1:nrow(meuse), meuse@data) Create a data.frame "df.new" with "IDS" (note different name) and "y" ...


4

Something very similar came up in the GIS Chat Room yesterday where @Kersten suggested reviewing Integrating external programs within ModelBuilder: This section describes how to integrate R, an external statistical package, within the ModelBuilder environment. In the example described below, there are two scripts needed to execute R functionality: a ...


4

Without your original data, I can't be sure this will work, but I thought it might help you out. I didn't bring it all the way there, this solution still likely needs some level of automation, but might give you a general way forward First, I create some spatial polygons polypoints1 <- matrix(c(1,2,2,1,1,2,2,1,1,2),ncol=2) polypoints2 <- ...


3

This is a reasonably simple problem to solve in ArcGIS. Open the Shapefile (it's not a "raster shapefile" btw; there's no such thin). Add two fields to your shapefile. One for Lat, one for Lon. Both should be of type Float or double. Populate one of these fields with the X, and one with the Y value for the point. (Using "calculate geometry"). Now, create ...


3

You want to do a left outer join on the shapefile's data.frame (data db) with the new table (table data). That will keep all the existing rows of your 'data db', join the appropriate fields together, and fill in missing data for rows that don't match. I'm assuming you're reading the shapefile as an sp object (using rgdal package in my example). I also ...


3

Take a look at the raster function in the raster package. It will let you create a raster with a specified extent, number of rows/columns and resolution. Here I will use characteristics of your data summary to create a 100x100 raster within the specified extent. I am passing an extent object to define the x and y limits. You can also use the specific ...


3

If I understand your question you would like specific values assigned to each resulting raster buffer. This can be done by simply assigning values via a index bracket assignment. First, create a raster with a polygon to buffer. library(raster) p <- spPolygons(rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))) r <- raster(ncol=100, ...


3

It is unclear on what you mean by columns. As far as I can tell you are subsetting and ploting rows, not columns. Your "ranks" variable is not a factor so, "nlevels" is returning 0. Replace "nlevels" with "length(unique())" or coerce the data to a factor using "as.factor" rather than "as.numeric" around your call to "lapply". The reclass function has a ...


3

Here's a suggestion: # Use rgdal, better library library(rgdal) coastline <- readOGR(dsn=".", layer="coastline") wets <- readOGR(dsn=".", layer="wets") Rasterize 'em: require(raster) # Create a generic raster, set the extent to the same as wetlands r.raster <- raster() Use extent (from raster package) to read bounds of vector and assign to ...


3

You can adjust the extent using the xlim and ylim parameters. library(rgeos) p1 = readWKT("POLYGON((0 1,0.95 0.31,0.59 -0.81,-0.59 -0.81,-0.95 0.31,0 1))") par(mfrow = c(1,2)) plot(gBuffer(p1,width=-0.2),col='black',xlim=c(-0.5,1),ylim=c(-0.5,1), xlab="Original") plot(gBuffer(p1,width=-0.2),col='black',xlim=c(-1,1),ylim=c(-1,1), xlab="Zoomed Out")


3

Provided I understood your question correctly, here is a suggestion. I made up some data for population, so if those are different for you the aggregate function call might need to be adapted. library(rgdal) library(sp) # shapefile tmpdir <- tempdir() ...


3

Here's a trivial example, which generates a single 'hotspot' around one region: # Generate coords extent of Tanzania x <- runif(1000, min = 30, max = 40) y <- runif(1000, min = -12, max = -1) # Assign some values, make data frame xyvals <- runif(1000, min=0.1, max=0.5) xydata <- data.frame(x, y, vals = xyvals) # Generate a hotspot centred ...


2

help(Snow.deaths) says: The scale of the source map is approx. 1:2000. The ‘(x, y)’ coordinate units are 100 meters, with an arbitrary origin. so I don't know why you think EPSG:3109 is going to help you, that seems to be Jersey transverse mercator... http://epsg.io/3109 To properly georeference these you'll need to find the lat-long of the arbitrary ...


2

For this to work, prop, Lat and dem need to be present in grd as well. Also, from your question we cannot infer whether dem is present in st.


2

So far I have found one fairly decent looking workaround: The packcircles R package may have been designed for another purpose, but it does a nice job pushing the points away from each other. I might not understand all of the inner workings of this package, but luckily, as you will find, the example file from the website can be used almost directly - all one ...


2

Yes, this is true. Your professor is right. Look for function fn_exponential in the source code if you want to be sure.


2

Do you want to change the extent of the actual data or just the plot? If it is just the plot you can use ylim and xlim to define the extent of the plot. Create some example data library(sp) x <- SpatialPointsDataFrame(SpatialPoints(cbind(runif(10, -115, -110), runif(10, 30, 45))), data.frame(ID=1:10) ) y <- ...


2

Merge is not doing what you want, probably because there are only two unique values in HJ or ZX. The merge columns should uniquely identify each record in oikismo and df1. Change your by.x and by.y appropriately.


2

If you do not need to apply a conditional subset, you can just use a bracket index on rows, without referencing the @data slot. library(sp) data(meuse) coordinates(meuse) <- ~x+y # Subset first observation and plot p <- meuse[1,] plot(meuse, pch=20) plot(p, pch=20, cex=1.5, col="red", add=TRUE) # display elev value for each observation for(i ...


2

Here are some functions to solve this issue. I use the "Cumulative Proportion" as a guide how many local principal components to keep. Just like global PCA, I define the percentage of variance and then select the local componets which cumulatively accounts for 85% and more variance on example data we would like to keep. cum..prop.var <- ...


2

I think another way of doing this is to just add the text file to ArcMap, display by x,y coordinates, and join to your to your original shapefile. You can then export the resulting shapefile (or just the table) as whatever you want.


2

You have to use the subset method (see ?subset.Spatial): subset(lines, X > 400 & Y=="YES") Alternatively you can use indexing operations via []: lines[lines$X > 400 & lines$Y=="YES", ] Your dplyr code filters just filters the data frame, but not the geometry.


2

I played with assign(), ls(), and mget() to accomplish something that I believe will improve your workflow. First I use ls() to get a list of all environment variables that start with "p": names_poly <- ls(pattern='^p.') I used combn() to find all the unique combinations of polygons combos <- combn(names_poly,2) I looped over combos using get() ...


1

There's a better library for reading really large shapefiles - fastshp. Doesn't seem to be available in repositories but the .tgz binaries are here. Here are the results for rgdal and fastshp with a 130MB shapefile with 32,545 features: library(fastshp) library(rgdal) system.time( test.fastshp <- read.shp("tz-landcover-ge.shp") ) # user system ...


1

Robin Lovelace has provided a nice little function to download a ggmap object and convert it to a raster. Using this you could do: library(ggmap) library(raster) library(rgdal) # courtesy R Lovelace ggmap_rast <- function(map){ map_bbox <- attr(map, 'bb') .extent <- extent(as.numeric(map_bbox[c(2,4,1,3)])) my_map <- raster(.extent, nrow= ...


1

If you do not need to merge topology, but just add new polygons, you can simply use: ab <- rbind(a,b) If you get a "non-unique Polygons ID slot values" error it means that the rownames of the objects are the same. To fix this you can use spChFIDs to change the rownames and associated slot relationships. Since the slots in the object use the rownames to ...


1

extract pulls raster values out based on an intersection with a vector. That could be helpful if you wanted to sample raster values given something like a points shapefile. In this case, however, you want a to identify the cells based on value, and then get the data directly from the raster. This may need to be modified if performance in reading the data is ...


1

If I understand your desired outcome correctly, you would like a count (richness) of species for each grid cell in the defined raster. I cannot speak to the differences between R and QGIS but I came up with a much more optimized and faster way to conduct your analysis. I leverage the raster package and use a raster stack to accumulate species. The workflow ...



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