Tag Info

New answers tagged

0

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


0

I agree with Micky T, it looks like you are using the wrong proj4 string. First, you will need to ensure which projection your numbers are in. WGS84 is the datum, not the projection. One possible projection is UTM, which defines your starting point for your easting and northing in meters. The x and y values that you give sure look like they are a ...


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


1

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.


0

Depending on whether you want a mean or a spatially weighted mean there are a couple things you can do using the rgeos and raster packages library(rgeos) proj4string(grid) <- proj4string(nuts) # I assumed these were the same projection? First find out which grid cells intersect your NUTS polygons grid_nuts <- gIntersects(grid,nuts,byid = TRUE) ...


0

In a SpatialPolygonsDataFrame, the polygons and the dataframe should automatically link up through the names of the polygons or the row.names of the dataframe. Assuming your SpatialPolygonsDataFrame is called spdf take a look at: row.names(spdf@data) names(spdf@polygons)


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

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


1

In the past I've used data generated using a Gibbs Sampler which is described at http://www.geog.leeds.ac.uk/people/a.turner/publications/archive/BrunsdonEtAl1999.html


0

This it's very easy using r package sp. But you should give an example of your data. I'm going to suppose that you have a table with the points coordinates. And you can save this table on CSV format. library(sp) test <- read.csv("test.points.csv") # We create a Spatial points data frame for our data, where X and Y are the fields with coordinates ...


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.


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


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


0

Super easy solution provided by @mdsumner: library(sp) library(raster) library(rgeos) library(spatstat) library(rgdal) library(maptools) setwd("C:/...") a<-readOGR(dsn=getwd(), layer="pol.a") b<- readOGR(dsn=getwd(), layer="pol.b") # use union in {raster} package ?raster::union ab<-union(a, b) resulted in : class(ab) [1] ...


0

All the g.. functions seem to discard the data by design, and you have to fiddle them back into a spatial dataframe (Don't ask me how...) Try using the union from the raster package. Haven't tried union, but intersect etc. worked as expected, producing spatial dataframes.


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


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


0

(1) Your plot shows the points as small circles. In the computation, they are considered as real, zero-area points, and the lines are considered as zero-width lines. If the plot suggests they intersect, this may be because the lines have a certain width and points a certain area, but that their mathematical representations do not intersect. You could compute ...


0

Others users might be interested to know that I managed to solve this. In order to extend the lines outside of the polygon (problem a): # Save centre point co-ordinates separately A.x <- Centre@coords[1] A.y <- Centre@coords[2] Using the great answer posted here: http://stackoverflow.com/a/7741655/4312169 I applied the code: # Generate points ...


0

Other users might be interested to know that I have managed to solve this. Rather than extracting the points directly, I used gIntersection in the rgeos package to truncate the lines to within the polygon, and then extracted the outer point. # Cut radial lines at polygon edge Centre_edge_lines <- gIntersection(Radial_spLines, Hull_poly, byid=T) # ...


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


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


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


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


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


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


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.


0

You should note that str is telling us that you have 25413 observations, not 200, in your "a" @data (data.frame) object. You could have a multi-part shapefile where numerious records are associated with each spatial location. If the results of dim(a) and dim(a@data) differ this is likely the case. This is why you provide as much information as possible in ...


0

This is how you can reproject a raster in R using the raster package. In this example, the input geotiff was in a NAD83 geographic coordinate system and I reproject to a NAD 83 UTM 15 projected coordinate system. A good reference for Proj4 format projections, which are used by RGDAL, can be found at spatialreference.org. library(raster) # Create ...


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.


1

You may want to consider to tell us why the above code didn't work, i.e. include the referring error message. I just tried to reproduce your issue using the German border from the countriesCoarse dataset included in the rworldmap package and everything worked out fine. Remember to set the projection and the desired extent of your raster template properly, ...


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


1

Use the function nowrapSpatialPolygons at the anti-meridian of your central longitude. This doesn't work with wrld_simpl from maptools, but it does work with countriesLow from rworldmap. library(rworldmap) data(countriesLow) library(maptools) prj <- "+proj=robin +lon_0=-198 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" x0 <- ...


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


2

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


1

Sounds like cleaning the geometry up will get you on the road to your non-raster solution ... herewith something of a kludge which does help with fixing bad geometry: # Load the library and problematic data library(rgeos) load("oneImage2_spdf.Rdata") >gIsValid(polysData, reason = T) Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, ...


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


1

Since you put GRASS into the tags, here's a solution based on GRASS: First, you need to know exactly what coordinate system the original data is in (as always with GRASS). I see that the *.prj file contains "TRANSVERSE MERCATOR" but it's not one of the standard UTM zones. Since you have not mentioned the EPSG code of this data, here's the proj4 string for ...


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


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


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


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


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


1

You need to perform the following actions: Read shapefile Create a function to reclassify your distance values Apply the function and write to a new field in the attributes "rank" Write the results to a new shapefile Of course you could simplify this if you only needed to read the shapefile attributes and work entirely within a data frame. ...


1

Depending on what you want to do with the data afterwards, you may want to have a look at the source code of dismo::gmap. As you will notice, the URL of the target image is stored in gurl, and after download, the temporarily saved *.gif is imported into R using raster::raster and inherits all the stuff that's required for subsequent spatial analysis (e.g. ...


5

This function will create a systematic point sample of a SpatialLinesDataFrame. The LID column in the results corresponds to the row.names of the source line object and can be used to relate the points to the lines. # Function for creating a systematic point sample of a SpatialLinesDataFrame # x SpatialLinesDataFrame # sdist ...


1

The SpatStat package has a function for creating points along lines at specified intervals. # maptools to read shapefile require(maptools) require(spatstat) # Load a line shapefile myLine <- readShapeLines("~/workspace/TEMP/myLine.shp") plot(myLine) # create a PSP object to use with spatstat myPSP <- as.psp(myLine) # Create points at 10 km ...



Top 50 recent answers are included