Hot answers tagged

3

This error comes from rgdal, I believe, specifically as it calls normalizePath, which probably should be avoided for filenames starting with http or ftp. normalizePath("http://somesite.com/image.jpg", mustWork=FALSE) on windows leads to: #[1] "E:\\home\\http:\\somesite.com\\image.jpg" but on linux I see: #[1] "http://somesite.com/image.jpg" So on ...


3

You need to transform your points to the coordinate system of the raster. You could warp the raster to your points' coordinate system but warping rasters is a bit messy. First create a SpatialPoints object from your coordinates and tell R it is in lat-long coordinates (4326): library(raster);library(sp);library(rgdal) pts2 = SpatialPoints(points,...


2

In the future please provide information on what you have already tried and a reproducible example. Whereas this Q&A site is a great resource we are not a coding service. The rgeos library provides functions for overlay/intersection analysis. There are important nuances in the type of function (covers, within, contains) that is appropriate here. There ...


2

This is not a bug but an expected behavior. Note that rasterize takes a fun argument that handles grid cells with two or more values. By default it uses that last function, namely the value that appears last on the data data.frame is used. Similarly first will use the value that appears first. Other functions include count, mean, etc. Here is a short ...


2

+proj=longlat +datum=WGS84 is a coordinate reference system (crs). If you would create a map with that, by treating longitude as 'x' and latitude as 'y', you could call it the equirectangular projection. +proj=laea stand for "Lambert Equal Area". This is planar (Cartesian) crs (unlike angular longlat), and can be used with rgeos::gLenght. gLength computes ...


2

You might want to ask a more general question of the stats of this on Cross Validated. If you look for significant differences between multiple pairs you need to correct for family wise error, to avoid this problem: https://www.xkcd.com/882/. Basically, at a 5% significance you should expect a false positive every twenty times. For the specifics of your ...


2

R would be a good place to do this, if you're willing to copy a little code. I think in general you are looking for a raster function, essentially. If you need it as a polygon, you can convert the raster to polygon when you're finished. library(rgdal) library(raster) library(plyr) # load your points points <- readOGR("points.shp", "points") # create a ...


2

To create a convex hull you can use function chull from grDevices. NB: This function returns the indices of the points in the input, not the points itself. Furthermore I was using it to draw lines, so I had to add the first point again as last point. The code I used looks like this: hull <- chull(topo) hull <- c(hull, hpts[1]) hullpts <- topo[hull,...


2

I recommend Python scripting with QGIS. The library is named PyQGIS. It has Vector and Raster GIS functions. QGIS is also a desktop GIS application that you can manually work with to understand various workflows. Here is a nice online cookbook containing many code snippets - http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/


1

coordinates(new_sp_grd) returns the center coordinates of each cell. Here is a simple example: library("sp") grid <- SpatialGrid(GridTopology(c(0, 0), c(1, 1), c(5, 4))) # get coordinates of center coordinates of each cell grid_ctr <- coordinates(grid) plot(grid, axes=TRUE) points(grid_ctr, col="red", pch=16)


1

Edit this snippet into your code: library(rgeos) # to use dissolve polygons library(broom) # to replace fortify() with tidy(). # Select Affected states based on $T3. Subset and merge to one polygon. statesAffectedList <- as.character(unique(final.plot$NAME_1[final.plot$T3 == "State Affected"])) states.shp$STATE_NAME <- as.character(states.shp$...


1

0 (and 360) represents North, as long as a pixel "above" another pixel is indeed north of it --- that is almost always the case. However, if you have a rotated data set (a rare thing) such that N is not up, you would have to apply that rotation to the output to get the true direction. This is not stated in the documentation because it seemed obvious, but ...


1

It took me a bit to get this too. # shapefile path example: # C:/Users/User/GIS/MyShapefile.shp # OGR command: library(rgdal) readOGR("C:/Users/User/GIS", "MyShapefile") You have to leave off the ".shp" extension as well as the last forward slash on the path.


1

You could treat this as a (degenerate) interpolation problem, using a neighbourhood of one, try something along the lines of library(gstat) i = idw(var~1, old, ggg, nmax = 1, maxdist = 500) library(sp) spplot(i[1]) where var refers to the name of the variable you sampled.


1

Transform the point to the coordinate system of the raster. Make the point into a SpatialPointsDataFrame in the lat-long coordinate system: > e <- data.frame(x=148.1, y=-35.6, id=1) > coordinates(e)=~x+y > projection(e)=CRS("+init=epsg:4326") then transform to the coordinates of the raster: > et = spTransform(e,projection(r)) > plot(r)...


1

Load your shapefiles into R as SpatialPolygonsDataFrame objects - lets call them p1 for the year 2000 and p2 for the year 2010. Drop any rows in p1 with fips codes that aren't in p2 - these would be county fips codes that don't exist any more, so their geometry has changed. Also, drop any rows in p2 that don't have a fips code in p1. These are new fips ...


1

The warning message is quite clear here: 1: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") : Spatial object 1 is not projected; GEOS expects planar coordinates It requires planar (Cartesian) coordinates, i.e. in meters, miles, etc. You however use polar reference system (WGS-1984). You should re-project your SpatialPointDataFrame to a ...


1

Set the values via the index, once you tabulated it: r <- raster(ncol=5,nrow=5) r[] <- 0 vec <- c(1,1,1,3,4,5,6,7,8,8,8,9) tab <- tabulate(vec, ncell(r)) r[vec] <- tab[vec] That's wasteful for large rasters, but we need more details about how this needs to be done if that's a problem. This is also a great example of how painful some ...


1

This problem is inherently an NxN comparison and a Mantel or Partial Mantel are quite inappropriate here. A Mantel is a pairwise matrix correlation and is entirely dependent on distance (ie., ecological, geographic). I have a function "mwCorr" in the spatialEco package that implements a moving window version of the Dutileul modified t-test, accounting for ...


1

Great question you've written here. There are two problems with your code. The first and most crucial is shown below: states.spdf$frac.pop = states.spdf$total.pop / (states.spdf$ALAND10+states.spdf$AWATER10)*raster.res You are assuming that population is equally distributed in space, however you treat density wrong. Instead of using raster.res to ...


1

Example data (something you should provide) library(raster) n <- 100 set.seed(0) x <- runif(n) * 360 - 180 y <- runif(n) * 180 - 90 xy <- cbind(x, y) p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20)) hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20)) p1 <- list(p1, hole) p2 <- rbind(c(-10,0), c(...


1

There are a lot of different libraries for geoprocessing. Depending on your particular needs could be convenient to use one or another. Even though, it is possible to combine many of them using using a common format or datatype. Here is a small list of the libraries I'm using at the moment. Python/GDAL The so called "Swiss knife of GIS" https://...


1

I just solved the problem following these steps: open the Terminal, digit touch ~/.bash_profile; open ~/.bash_profile and add two rows: export R_HOME=“/Library/Frameworks/R.framework/Versions/3.3/Resources” export PYTHONPATH="/Applications/QGIS.app/Contents/Resources/python" save the file and close the Terminal open a new Terminal digit sudo nano /etc/...


1

This is caused by the matrix library used by gstat. Historically (gstat was released as open source code in 1997) it used the LDLfactor routine in the meschach library. Around 6 months ago I factored this out this code and replaced it with the BLAS/LAPACK which are native in R. LAPACK uses Choleski decomposition. LDLfactor allows for some non-positive ...


1

You can do library(raster) x <- shapefile('file.shp') crs(x) x$area_sqkm <- area(x) / 1000000 Assuming that your crs is longitude/latitude, or with meter as distance unit


1

I believe that a call to "data.frame" is all you need. If you have several rasters, just stack them. If you want the rasters as columns and the statistic as rows you can just transpose on the fly (see second example) library(raster) x <- stack(system.file("external/rlogo.grd", package="raster")) ( x.stats <- data.frame(x.mean=cellStats(x, "mean")) ) ...


1

Yes, this is correct. When you print the model by typing model.vari you'll see sill values, split up in a nugget component (the offset) and the exponential component. The sum of these two is usually indicated by "the sill value" (i.e., around 25).


1

What exactly are you wanting to compare? A point pattern is a representation of an explicit spatial process that is significant from a spatial random assumption. A raster does not meet the same criteria of a point process. You can test the similarity of values but, this does not at all demonstrate the equivalency of an underlying spatial process. Think of ...


1

I could apply an Elliptical Inverse Distance Weighting (EIDW) in R using the function gdal_grid of the gdalUtils package. I modified a little bit the gdal_grid help example in R documentation as follows: # Load packages library('sp') library('gdalUtils') # Load data data('meuse') coordinates(meuse) = ~x+y # Create a properly formatted CSV: ...


1

I found an answer on the following website on how to do this in R. The solution that is given above also works, but I am having memory problems to run it. http://article.gmane.org/gmane.comp.lang.r.geo/24010/match= ########################################################################## library(fields) library(chron) library(ncdf4) setwd("C:/Users/...



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