# Tag Info

5

Create some simple test rasters in R: > m=matrix(1:9,3,3) > m2 = matrix(c(9,2,3,4,1,5,6,8,7),3,3) Then we can trivially compute the covariance between these matrices: > cov(c(m),c(m2)) [1] 2.125 and I would wager doing the computation by hand would get the same result. What does layerStats do? > D = stack(raster(m),raster(m2)) > ...

4

To read your shapefile, i recommend you to use rgdal package and its readOGR function, or eventually use readShapeLines from maptools package. These packages rely on the sp package as concerning how the geospatial data is structured in R. You can do easily this to convert your shapefile into data.frame (ie extract the attributes of the shapefile) ...

4

I use QGIS 2.2 so what I post may differ in certain areas. You can install the Processing plugin via: Plugins > Manage and Install Plugins... Once installed, you should see the Processing menu in your toolbar. Select Options and Configuration: Select the Providers list, find R scripts, check the Activate box and set the path of the R folders: Next ...

3

It would be very helpful for you to read up on some R basics, particularity pertaining to sp class objects. A very good starting point would be Bivand's ASDAR book and the sp vignette. Here are some other related R spatial analysis introductory material. As to your problem at hand. For one, you can easily create a variety of spatial weights matrices in ...

3

Not sure whether it is what you are looking for (it depends on the action you want to perform). For example, let's use a shapefile of populated places. library(rgdal) shp <- "C:\path\to\shapefile" pts <- readOGR(shp,'ne_10m_populated_places_simple') OGR data source with driver: ESRI Shapefile Source: "C:\path\to\shapefile", layer: ...

2

pardon any syntax errors...hope this helps. [edit: adding some more stuff to make clearer, as per suggestion] #import the numpy and gdal libraries import numpy as np from osgeo import gdal #an empty array/vector in which to store the different bands layers = [] #open raster ds = gdal.Open('raster.tif') #loop thru bands of raster and append each band of ...

2

I'm not sure where you get the above-mentioned 'attribute' variable from, so I just assumed rownames(sp_df@data) (ranging from 0 to 193) as the desired output. Anyway, here is my approach that is probably not the fastest (or the most convenient in general), but it works. Basically, the script starts to loop through all the cells of the generated raster ...

2

The gIntersection function from the rgeos library might help you. See the commented code below. library(sp) library(rgeos) pointB <- SpatialPoints(cbind(1,1)) pointC <- SpatialPoints(cbind(4,2)) distanceA2B <- 2 distanceA2C <- 3 # create the circle polygons around the points with the distances polyB <- gBuffer(pointB, width = distanceA2B) ...

2

I'd convert the output to a raster object. Then: require(spatstat) require(sp) require(raster) set.seed(1985) x <- runif(20) y <- runif(20) p <- SpatialPoints(coords = matrix(c(x, y), ncol = 2)) plot(p) Then compute densities: pp = ppp(x,y) # all points in a (0,1) default window d <- density.ppp(pp, sigma = 0.1) dp <- density.ppp(pp, sigma ...

2

The problem is that the data was never in UTM to begin with, and so by having a UTM projection, the file was ultimately being told to be something it wasn't. (Such is life) :) Reprojecting it doesn't fix the problem, because the transformation math is based on coordinates that don't match the assigned projection. To fix this I deleted the .prj file, and ...

1

Please read the sp vignette on spatial classes and methods. vignette(package="sp")[4] vignette("intro_sp") Since there is a slot (@data) that holds a data.frame related to the sp object, no coercion is required. class(foo@data) str(foo@data) ( df <- foo@data ) However, it is good practice to operate directly on the @data slot rather than pulling ...

1

If you would check help for as.im() you would notice that raster data type is not supported. We should thank @mdsumner for bringing up my old question (I do not remember the occasion that forced me to ask it - suppose it was just some testing and I wouldn't recall that question by myself). I used there a geostatsp package that provides needed functionality. ...

1

To figure out which file has a problematic projection and get additional clues about how to fix it, you can plot the maps on an ArcGIS basemap. This will get the shapefiles to appear on the same screen so you can see which is problematic. In this case you can quickly see that the SourceA file is plotted off in the Atlantic Ocean far from it's proper ...

1

To subset or operate on an individual object all you need is the "row" index, which corresponds to the slots. In this way it is independent of feature class type (i.e., point, polygon, pixel). # To iterate through the feature class for(i in 1:nrow(shp)) { p <- shp[i,] }

1

Here is an example with the 25-km dataset, for January 2014, followed by an example for the 4-km dataset, January 2014 (improvement are welcome). The 4-km version is read using the example given here. You need to download functions from the same website: val2col, image.scale and earth.dist. # Required librairies library(rasterVis) library(rgdal) ...

1

I've analysed the geometry issues in the attached data, and it seems it does not ONLY have orphaned holes but also geometry validity issues. It's true that an orphaned hole is somehow a geometry validity issue, but rgeos does not handle it in the same way, as for orphaned holes, an error is raised, instead of a simple warning. As you indicate, they are hints ...

1

A 32 bit Windows 7 Operating system can use at most 4 GB of RAM[1], combined for all activity on your machine. If your shapefile is very large and you also have many other programs running then you'll probably exceed 4 GB. So, to see which process uses up your memory you can open the Windows Task Manager (using ctrl-alt-delete) and then look at the ...

1

Here are some ideas library(raster) library(rgeos) DK1 <- getData('GADM', country='DNK', level=0) # extent coordinates e <- c(8.41, 9.04, 56.92, 57.16) DK_reg1 <- crop(DK1, extent(e)) # resolution to get about 1000 cells res <- sqrt((e[2]-e[1]) * (e[4]-e[3])) / sqrt(1000) r <- raster(extent(e), res=res, crs=crs(DK1)) # rasterize if ...

1

KevinMayall's suggestion of intersecting circles is the easiest. A good theoretical treatment can be found at Wolfram Mathworld, which may get into more detail than you want to get this working (but is nice background for writing a paper about your methodology). This StackOverflow answer has a broken link but outlines a basic approach. First ...

1

You need to iterate through the files, so your loop should be as follows: for file in $NPPFILES do r.out.gdal input=$file output=\$file.tif format=GTiff done You should add the semi-colons if you write the loop on one line.

1

Cant you just simply use text to fill your proj4string. So: First read in your shapefile calling it: pointshapefile. And then run the following line. pointshapefile@proj4string <- "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs" Afterwards you can use spTransform. I think R gives you a warning that you have ...

1

The NA's here refer to the fact that the min and max value are unknown (this is no longer the case in the current version of 'raster'). Your function should work if you use 'calc' instead of 'overlay'

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