Tag Info

New answers tagged

0

Change geom_polygon(aes(x=long, y=lat, group=group), fill='grey', size=.2,color='green', data=data.shape, alpha=0) to geom_polygon(aes(long, lat, group=group, fill='grey', size=.2,color='green', data=data.shape, alpha=0)) m = map + geom_polygon(data=data.shape, aes(long,lat,group=group,fill="grey", size=2, color='green',alpha=0)) m


0

To expand on @Chris's answer, use raster::writeRaster to write to any of a large variety of raster formats. Depending on whether your build of GDAL, you should be able to write to these (and possibly more) filetypes (returned by raster::writeFormats()): # name long_name # 1 raster ...


1

If the square can be in meters instead of degrees, you can create a custom CRS in oblique mercator projection: +proj=omerc +lat_0= 21.04247 +lonc=88.38045 +alpha=-70 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs And create a vector grid with that CRS:


0

Take a look at ?text. You can place the labels using the polygon coordinates which will return approximate polygon centroids. library(sp) Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2))) Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2))) Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5))) Srs1 = Polygons(list(Sr1), "1") Srs2 = Polygons(list(Sr2), "2") Srs3 = ...


2

It is very unlikely that this problem follows parametric assumptions. I would recommend exploring nonparametric group tests, available in R, such as: Mann-Whitney (wilcox.test), Wilcoxon Signed Rank (wilcox.test), Kruskal Wallis (kruskal.test), Friedman Test (friedman.test). You could also apply a Monte Carlo sampling approach, lbl_test in coin or ...


2

There are several ways you can tackle this in R, including spDists in sp and gDistance in rgeos. An efficient way, that is expandable to multiple kNN ID's and distances, is to use spdep. require(spdep) data(meuse) coordinates(meuse) = ~x+y meuse <- meuse[1:10,] meuse@data$IDS <- 1:10 # Neighbor row indices and add neighbor attribute ID's ( ...


-2

EDIT Should have been posted as a comment...Didn't want to delete the other comments to this comment...but the selection of tests should be directed towards the nature of the data itself, even @MikeRSpenser suggests in the second last line that the t-test is unlikely the best test...hopefully the OP read the comment ORIGINAL bearing in mind that landuse is ...


1

The R software is a good solution. You can read your rasters, extract values for the sub areas and then run an ANOVA, t-test or whatever else you might fancy. Remember to check whether your data are normally distributed so you can pick your test. Here are the basics: # read rasters slope = raster("~/dir/dir/file.tif") land = raster("~/dir/dir/file.tif") # ...


3

There is a border argument in the plot-method: plot(x, border="grey") R> args(sp:::plot.SpatialPolygons) function (x, col, border = par("fg"), add = FALSE, xlim = NULL, ylim = NULL, xpd = NULL, density = NULL, angle = 45, pbg = NULL, axes = FALSE, lty = par("lty"), ..., setParUsrBB = FALSE, usePolypath = NULL, rule = NULL)


2

The short answer is that you can aggregate SpatialPolygons by one or more fields with raster::aggregate. The long answer is that this is how I'd approach your entire problem with either base plotting, or lattice: Data preparation: library(rgdal) library(raster) # Read in data on major cities cities <- read.csv("ciutats.csv", header=FALSE, ...


0

it's good, I managed to do it on Ubuntu Virtual Machine thanks to @Pascal code: library(ggplot2) library(gstat) library(sp) library(maptools) library(rgdal) library(zoo) library(xts) library(RPostgreSQL) library(spacetime) library(foreign) library(RODBC) library(rgeos) library(Rcartogram);install_github('chrisbrunsdon/getcartr',subdir='getcartr') ...


2

If I understand you correctly, you want to test sample variation. That is to say, how well your sample distribution matches your population (raster) distribution. This is not a correlative relationship and is commonly done by comparing the sample mean and variance against the population. Here is an example where I calculate mean, variance and quantiles for ...


0

It turns out that sp::disaggregate can be used to pull apart the singlepart polygons. pts <- SpatialPoints(cbind(c(1, 1, 2, 3), c(1, 2, 1.5, 2.5))) b <- gBuffer(pts, width=0.6) over(pts, disaggregate(b)) # [1] 1 1 1 2 in_hole <- SpatialPoints(cbind(1.375, 1.5)) over(in_hole, disaggregate(b)) # [1] NA (raster::aggregate can be used to combine ...


1

Just coerce "Points" to a ppp object using "w" as the window. w <- as.owin(StudyArea) pts.ppp <- as.ppp(coordinates(Points), w) plot(pts.ppp)


2

It seems like your data has some topological issues (like overlapping polygons, multi or micro polygons inside the largest one). I tried to do it with the original data in QGIS and it had the same result, so I ran a v.clean on a data subset (because original datasets are too large for testing purposes) and using this topology cleaned layer the intersection ...


1

jbaums solution is good, but in cases where you cannot make a stack you can do something like the below. Just a single loop (not three or four!). Pau's solution (first making a table) is needlessly complex, and involves too much manual labor, and cannot be easily applied to other datasets. # read files setwd("D:/Data/LANDSAT") library(raster) library(rgdal) ...


1

"Locating pixels" can assume many meanings so, I am just taking a guess at what you are after. I recommend using the raster package as this will keep the problem memory safe and provide a number of functions specific to operating on raster arrays. Let's create some data (note; you can use "raster" or "stack" to read raster format files on disk). ...


1

You can load data from postgis with the rgdal package into R. library(rgdal) library(sp) dbname = "yourdatabase" host = "yourhost" user = "AUser" pass = "ThisUsersPassword" name = "ASpatialTable" # Postgis table dsn = paste0("PG:dbname='",dbname,"' host='",host,"' user='",user,"' password='",pass,"'") res = readOGR(dsn,name) plot(res) Write yourself ...


1

I would recommend using the raster or rgdal package(s) to read tiff files. This will also set you up for success in your second question, which not only needs to be a new question but also needs more than "how to locate pixels". In the raster package you can use "raster" for single band or stack for muntiband (e.g., RGB) tiff's. In rgdal you can use ...


1

Your description is a bit confusing. It sounds like you would like to aggregate the mean of values, located at points, based on the intersection with a polygon object. Even if this is not exactly the case, this code should at least point you in the right direction. require(sp) # Create example data data(meuse) coordinates(meuse) = ~x+y sr <- ...


0

Some background: Those messages come from libtiff, which is widely used for reading TIFF images. R is just one of the users. TIFF is a flexible format, and additional information can be added using extra "tags" (hence the name of the format), not all of which are necessarily public. The specification will give you more detail if needed. Different tools ...


1

This is almost a duplicate of this post, but you have an additional cropping step, so I'll post a new solution. Given your .img files all have identical extent and resolution, you can save a lot of hassle by stacking them from the start (you can pass a vector of file names to raster::stack). You can then crop the stack in one shot, and write them all out ...


2

I am sympathetic, this is a bit of a pain. You have to explode the slots, of the gBuffer polygon object, into individual polygons. require(sp) require(rgeos) pts <- SpatialPoints(cbind(c(1, 1, 3), c(1, 2, 3))) b <- gBuffer(pts, width=0.6) over(pts, b) ########################################################### # explodes slots into individual ...


0

I recommend you create a table (excel, txt, csv, etc) with the name of the raster you want to process and the folder where you want to save the output raster in different columns. You can use this script: #load library library(raster) library(rgdal) # #read table TABLE=as.data.frame(read.table("D:/table.csv", sep=";", header=T)) #initial value id=0 ...


1

Thanks to @gene and https://geoscripting-wur.github.io/AdvancedRasterAnalysis/ I can now answer my question (copied and modified): library(raster) # create some raster data r <- raster(ncols=12, nrows=12) set.seed(0) r[] <- round(runif(ncell(r))*0.7 ) r[r==0]<-NA # extend r with a number of rows and culomns (at each side) # to isolate clumps ...


1

#reproducible example r <- raster(ncols=12, nrows=12) set.seed(0) r[] <- round(runif(ncell(r))*0.7 ) rc <- clump(r) #extract IDs of clumps according to some criteria clump9 = data.frame(freq(rc)) clump9 = clump9[ ! clump9$count < 9, ] #remove clump observations with frequency smaller than 9 clump9 = as.vector(clump9$value) # record IDs from ...


1

This can be achieved with raster::extract in R: library(raster) s <- stack(setNames(replicate(4, raster(matrix(runif(100), 10))), c('elev', 'aspect', 'slope', 'curv'))) Create a matrix of x,y coordinates: xy <- expand.grid(x=seq(0, 1, 0.15), y=seq(0, 1, 0.15)) Sample the layers of s at each of the points in xy: z <- ...


5

Any truly general-purpose effective method will standardize the representations of the shapes so that they will not change upon rotation, translation, reflection, or trivial changes in internal representation. One way to do this is to list each connected shape as an alternating sequence of edge lengths and (signed) angles, starting from one end. (The shape ...


1

A good method to compare these polylines would be to rely on a representation as a sequence of (distances, turn angles) at each vertice: For a line composed of points P1, P2, ..., PN, such sequence would be: (distance(P1P2), angle(P1,P2,P3), distance(P2P3),... ,angle(P(N-2),P(N-1),PN), distance(P(N-1)PN)). According to your requirements, two lines are ...


1

You're asking for a lot with arbitrary rotation and dilation! Not sure how useful Hausdorff distance would be there, but check it out. My approach would be reducing the number of cases to check via cheap data. For example, you could skip expensive comparisons if the length of the two linestrings is not an integer ratio (assuming integer/graduated scaling). ...


2

Your question is not at all clear. I am assuming that you want to operate on multiple shapefiles and not just the contents of a single shapefile. In R you can write a for loop for reading, modifying and writing shapefiles. Here is an example where I read each shapefile, apply an ifelse statement to classify a defined field and write the results. I am ...


1

From section 10.13 of R for Mac OS X FAQ: When executing system commands (for example directly via system or indirectly via functions that call other programs such as install.packages) the locations in which the shell is looking for programs is governed by the PATH environment variable. That variable may be set differently for R started from an ...


1

If you have a list of raster objects with equal extent and resolution, it's probably easiest to stack them first and then use writeRaster with bylayer=TRUE. For example: library(raster) Create some dummy data L <- setNames(replicate(3, raster(matrix(runif(100), 10))), c('A', 'B', 'C')) L # $A # class : RasterLayer # dimensions : 10, 10, ...


1

It looks like there are some issues with the path variable, i.e. the shell opened by R doesn't know the path to the gdal binaries. There are two ways to fix this: Specifying the full path You can always use the whole path to gdalinfo in your system call to make it work: path <- "/path/to/gdal/bin/gdalinfo" system2(path, "--version") This may be the ...


1

I think this should be possible to do using ArcGIS for Desktop by: Using Create Fishnet (Data Management) to make label points that are located where you want them i.e. perhaps offsetting their origin coordinate by half a fishnet cell height and width. Using Extract Values to Points (Spatial Analyst) which: Extracts the cell values of a raster based ...


0

I found that by increasing the maxmemory and chunksize that I got modest speed increases. I had to do this conservatively, because the processing sometimes can take up 3-5 times as much memory as the maxmemory setting (it must be writing multiple files at that max memory setting). Here's the code for changing those memory settings: ...


4

Have a look at nlayers(s). The returned number of layers will equal 28 - at least for the above example with 4 multi-layer objects encompassing 7 layers each. Applying stack to multiple multi-layer files results in one huge 'RasterStack' object, i.e. all the single multi-layer objects are appended to one another. If you would like to have separate stacks ...


3

There's no need for any fancy tools here at all. Just a decent text-editor (notepad will do if necessary) will let you turn it into a ASCII text file in a few easy steps. Process Export your data as a CSV file (Excel/LibreOffice). (If you can export it with space as the delimiter, skip step 4). This gives me for my subset of 6 squares: ...


2

Sorry, The solution was very simple, I needed to convert to character before converting to numeric: commune_5m@data$NUMERO2 = as.numeric(as.character(commune_5m@data$NUMERO)) ### Create numeric attribute from factor


0

In R you do not need to use loops for things like this. Please use a reproducible example like the one below to ask questions. library(dismo) set.seed(1) s <- stack(system.file("external/rlogo.grd", package="raster")) p <- randomPoints(s, 100) d <- extract(s, p) boxplot(d) Or, what I think you want: boxplot(t(d))


2

The extent of the raster is equal to the extent of the cell centres, expanded by half the resolution. Here's an example: Create a dummy raster with extent c(0, 1, 0, 1) and resolution c(0.1, 0.1): library(raster) r <- raster(res=0.1, xmn=0, xmx=1, ymn=0, ymx=1) Extract cell centres: p <- rasterToPoints(r) head(p) # x y # [1,] 0.05 ...


1

You are using R syntax, why not just do this in R? Anyway, I have no idea how GME is parsing to R but from an syntax standpoint your backslashes either need to be double ("\") or single forward slashes ("/"). I believe that your actual issue is that you need a separator argument in the paste function. The way you have it now is creating a space between the ...


3

You can do this with the Intersect tool. Normally, performing an intersect with polygons will only return the overlapping area. But if you change the output_type from INPUT to LINE, then you can get just the collinear borders between polygons. If this is our input: And we change the output_type parameter: We get the green lines as output: The ...


1

You can do this with python if you are at least at ArcGIS 10.1. If you have ArcInfo, you can use the Feature to Line tool. Otherwise, you can use this simple script. This script does not support true curves though. The resulting lines, if topologically correct should overlap then you can run an intersect of this line fc on itself to find the boundaries ...


1

The split command requires a data frame or vector and you are giving it a SpatialPointsDataFrame. All you have to do is to do the split first and then create your SpatialPointsDataFrame. Below is the complete code for this. library(sp) library(dplyr) SB <- read.csv("SB_040914.csv", header=TRUE, sep = ",") SB <- SB %>% group_by(Fecha, CodPar, x, y) ...


0

The difference in your shapefiles is that the first shapefile is in UTM coordinates and the second is in lat/lon coordinates.


1

You can do a spatial intersection between your two datasets ("Al", "output") with the following code: library(rgeos) inter = gIntersection(Al, output, byid=T, drop_not_poly=T)


0

The regular RGDAL from the R package repository has had issues in the past - you can use the one provided by kynchaos.com: Download the DMG from the Frameworks page. Mount the DMG. In R, install the package from the included .tgz file: install.packages("/Volumes/rgdal/rgdal_0.9-1.tgz", repos = NULL) (substitute /Volumes/rgdal/ for another package if it's ...


1

The behaviour of gIntersection is not to pass any intersected data by design: Since there are no general matches between intersected spatial objects, any arbitrary operations on attributes require assumptions about unknown user intentions. This is why no data slots should be passed through ... ... The design of gIntesection() is inentional, ...



Top 50 recent answers are included