New answers tagged

0

It is quite unclear what you attempt to do. For instance summarize and feature are quite vague without further context. So I am guessing here. With your rounding of the coordinates you in effect moved the points in "square" cells to a central point for each cell. You don't need to. Just imagine the cells, and loop over them, calculating the aggregate you ...


1

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


0

I know this question is old, but I have just uploaded a Matlab toolbox that will compute width from either 1) input channel mask or 2) input banklines (your case). It has a number of other tools that might be helpful to you; check it out here. You will probably want to define a centerline (for parameterization of your along-stream variables i.e. width). If ...


0

I've tried out a bit and found that: RasterBricks, as mentioned by RobertH's answer, do work and are more user-friendly and easy to use; Rgdal methods like readGDAL also work, but with more parameters it's a little bit less user-friendly; So which option should one use? According to my tests (on my 420GB GeoTiff with dimensions of 18660x21592 and 374 ...


0

In R, this can easily be done using the mapview package. For instance, displaying Landsat 8 band 5 with a customized color scheme and breaks works as follows: library(mapview) ## custom color palette cols <- colorRampPalette(c("green", "blue")) ## visualize data m <- mapview(poppendorf[[5]], col.regions = cols(100), at = seq(5750, ...


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


0

As you have a single file, you should create a RasterBrick. That should make things faster as it could indeed benefit from the by pixel interleave. By creating a RasterStack you create a list of RasterLayers, i.e. you treat each "band" as a separate file.


0

The kde values from MASS::kde2d are stored as a matrix. The array coercion is to create a vector with the correct length. However, you cannot export an sp DataFrame object with an array in the @data slot and for some reason the function is not outputting a vector column. You can simply coerce the offending column into a vector using as.vector and overwriting ...


0

Here is one option: Use the Spatial Join method to join the polygon info to the points Create a new text field in the joined point layer Field calculate/concatenate the dcpWaterBodyType field values with the polygon name or id field Finally, right click on field created in step three and choose the Summarize option This should produce a count of the ...


0

Not sure what you have found, but one of the key players you want to look up is Roger Bivand. Couple examples below. There are many more, but if you poke around the references you will get there. Implementing functions for spatial statistical analysis using the R language. R Bivand, A Gebhardt - Journal of Geographical Systems, 2000 Approaches to classes ...


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.


0

I found I was able to solve the problem using the nncross() function in the package spatstat. I don't know if this was the simplest solution, since it required first converting my SpatialPoints to class ppp. In any case here it is: library('raster') library('spatstat') # Make the grid. Use raster() and then convert. # Not the most direct but it works. ...


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


0

This doesn't directly answer your questions, but some time ago I wrote a basic code to do what you're looking for. It extracts the orientation and inclination of RoofSurface polygons in CityGML, to estimate the solar irradiation of rooftops. The code is released on Github, so you might want to have a look. The part of the code that would probably interest ...


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.


0

I found a solution using the mutate function in dplyr package. Make sure you have a data frame with a column defining the months e.g m1, m2 etc df_Seasons <- df %>% mutate(season = ifelse(month %in% c(m1, m2, m3), "Season 1", ifelse(month %in% c(m4, m5, m6), "Season 2", ifelse(month %in% c(m7, m8, m9), "Season 3", "Error"))))


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


0

I created distances rasters to all 3 types of water bodies (cell size = 1000m*1000m) and used cell statistics (min) to find nearest distance. After that I ran zonal statistics as table and joined results back to polygons. This is how max distance looks like: What you call nuts is a set overlapping polygons, make sure to remove all with status < 3.


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


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


0

I really like using tmap for choropleth mapping in R. Here's the documentation. I think it's better for plotting maps. Here's an example of a function I once wrote to plot a choropleth map. function (function1_info, indicator, color_scheme, classes, classification, legend_title, map_title){ # The map plot will classify the data in quartiles, but the ...


0

Unfortunately, you cannot use spplot with base plotting. This is because spplot is a wrapper for lattice, which is a high level plotting engine. I have always found lattice obtuse and believe that it is becoming obsolete. If you are going to invest time in learning a higher level plotting engine I would recommend ggplot2. For you problem, there is no need ...


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

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


0

please provide a reproducible example. Perhaps this is because there are values in between 2 and 3 and 3 and 4? To sum the layers you should do sum(reclass) library(raster) s <- stack(system.file("external/rlogo.grd", package="raster")) s <- stack(s, s, s) ss <- reclassify(s, matrix(c(-Inf, 254, 0, 254, 255, 1), ncol=3, byrow=TRUE)) sum(ss) # ...


0

I am also currently looking for this functionality, using R and leaflet/mapview etc. In mapview, my current workaround is to convert the raster into a points layer then plot the raster and the points layer together, put full transparency on the points layer. When you hover over the raster cell you get the associated points value without even clicking (...


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

+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 angularlonglat), and can be used withrgeos::gLenght`. gLength computes the ...


0

Given that start with RasterStack you might consider writeRaster. Example data: library(raster) r <- raster(crs="+proj=utm +zone=32 +datum=WGS84", xmn=0, xmx=10, ymn=0, ymx=10) values(r) <- 1:ncell(r) s <- stack(r, r/2) This will write your RasterStack to a ncdf file: x <- writeRaster(s, filename="test.nc") It includes the CRS info. But ...


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


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/


0

There is certainly a neater approch, but this worked # convert SPDF to SLDF bound <- as(poly, 'SpatialLinesDataFrame') # overlay SLDF and SPDF bound <- raster::intersect(bound, poly) str(bound@data) 'data.frame': 43 obs. of 4 variables: $ ID_poly.1: num 6249 6249 6249 6249 6249 ... $ Code.1 : chr "33200" "33200" "33200" "33200"...


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


0

Mac Os X is an Unix system and the processing Python module search for an R executable because it uses the subprocess module to execute the script directly with R. The Python script involved is /Applications/QGIS.app/Contents/Resources/python/plugins/processing/algs/r/RUtils.py 1) it search for the R executable 2) it creates a R script with the command ...


0

That particular file, or at least that particular dataset within that file, is not a regular grid, and so can't be converted (easily) to a GeoTIFF file or read in as a raster data source. > long <- file["Soil_Moisture_Retrieval_Data"]["longitude"] > str(long) Formal class 'DataSet' [package "h5"] with 7 slots ..@ name : chr "longitude" ..@...


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

I also had this same error using gContains(). I didn't realise I had anything internal in my polygon. I dissolved the internal structure like this: spgeom1$rowNo = 1 spgeom1= unionSpatialPolygons(spgeom1, spgeom1$rowNo) and went on to use gContains() with no problem.


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


0

Addressing your exact output and expanding on the code you provided, you could query by date (the advantage of Posix) and index the rasters in a stack object. The use of "which" returns the raster stack/brick index meeting the query criteria. The use of "apply" is to sum rows, which represent pixels, to derive a seasonal sum to average. For large rasters ...


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

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


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


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


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


0

to loop through a list of raster file paths: library(raster) #load raster package rasterlist <- list("path 1","path 2",...) #create list of raster file paths outlist <- list() #create empty list to store outputs from loop for (i in 1:length(rasterlist)) { #for each raster in rasterlist r <- raster(rasterlist[[i]]) #read element i of rasterlist ...


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


0

After much searching, I think I no longer have any errors when running these tools. First, I deleted a few gdal installs on my computer and kept the 32-bit version. Then I followed the steps outlined here and set my environment paths to where my GDAL install was located. This fixed the error in rgdal but not when using gdalUtils. To get gdalUtils to work I ...



Top 50 recent answers are included