New answers tagged r
0
One approach is to write a function and apply it across all of your datasets using lapply. For example:
require(raster)
# Define the input and output workspaces
inws <- '/path/to/input_workspace'
outws <- '/path/to/output_workspace'
# Get a list of tif files from inws
files <- list.files(inws, pattern = "*.tif", full.names = TRUE)
# Function to ...
0
My mistake!
I forgot to include p1 into a list.
This code works
st_polygon(list(p1))
POLYGON ((0.4240843 0.6460535, 0.316801 0.6083381, 0.394125 0.851623,
0.4240843 0.6460535))
2
As I said in comment the error comes from your code and it means that somehow a raster with 0 cells has been created
This appended here
metrics = list(
Prom_class = integer(),
TreeCov = numeric())
if(length( z[z > 0]) <= 5){return(metrics)}
You returned nothing and internally grid_metrics tried to build a raster out of nothing. You can return ...
0
As I am not completely understanding what the exact output is you want (you want new "features" with every time the point and the five nearest (=500 new features), or do you want just to eliminate the points that are never the nearest, or just a dataframe with indicating the index of the 5 nearest point for every point, or....?), I will tell you how it is ...
0
Create a factor column before plotting, then you have full control over the order.
E.g. for your example:
vals <- meuse$zinc
# compute breaks (equivalent to `do.log = TRUE`)
breaks <- c(min(vals),
tail(head(exp(seq(min(log(vals)), max(log(vals)), length=6)), -1), -1),
max(vals))
categories <- cut(vals, breaks, dig.lab=4, ...
2
This is not how you replace elements of a vector:
if (myshapefile$varA== "1") {
myshapefile$d1 == myshapefile$income
} else if (myshapefile$varA != "1") {
myshapefile$d1 == "0"
}
An if test has to evaluate to a single true or false value to decide which block of code to execute. You are trying to run it on a vector of your varA values.
I guess you really ...
0
See here a reprex, the key functions are st_make_grid, st_nearest_feature and st_distance. One thing, to make all this work, I need to project the shapes. I used for this example EPSG:3857 for no particular reason, but you have to pay attention to the units on st_crs(yourshape), as the grid would use the same unit.
For the reprex I used countries as ...
2
You are reading your files in with rgdal. Read them in with sf (st_read) then you can apply an sf operation like st_join to it. Or change your objects from Spatial Point/Polygons dataframe to sf objects using st_as_sf(deciles) and st_as_sf(regions) and then apply st_join to it.
2
What "they" didn't teach you about R is that instead of "mutate":
deciles = mutate(mydata, quantile_rank = ntile(mydata$income,10))
deciles
which only works sometimes, you can always create the new column directly:
deciles$quantile_rank = ntile(mydata$income,10)
deciles
Your input is an sp class "SpatialPolygons Data Frame", which has been read in from a ...
3
What about st_sample(shp, n_points_wanted, type="regular" )? It gives you approximately 200 points, as it has some degree of randomness:
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 3.5.3
#> Warning: package 'ggplot2' was built under R version 3.5.3
#> Warning: package 'tibble' was built under R version 3.5.3
#> ...
7
STOP USING SO MANY PIPES!!!
n_grid <- (area_box / (area_box - area_shp ) )*n_points %>% sqrt()
is doing
n_grid <- (area_box/(area_box - area_shp)) * (n_points %>% sqrt())
in other words its square-rooting n_points. Fix that and I get exactly 200 points.
The pipe operator has high priority:
> 3 * 9 %>% sqrt()
[1] 9
> 3 * (9 %>...
1
I was able to find the algorithm name using find_algorithms()[grep("histogram",find_algorithms()] which led me to native:zonalhistogram instead of the expected qgis:zonalhistogram. And, it appears than a raster extraction that took 52 minutes using R functions, took only 94 seconds with zonalhistogram!
0
This is untested but here is how I would approach it:
library(dplyr)
library(purrr)
library(tibble)
library(sp)
#Get list of files
mydir="C:/path/to/textfiles/"
temp = list.files(mydir,pattern="*.shp$",full.names = TRUE)
fn = list.files(mydir,pattern="*.txt$")
temp=enframe(temp)%>%
mutate(savedir=fn)
temp = temp %>%
#load in your table
mutate(...
0
Thanks @Spacedman for you help -- I managed to accomplish what I was after. The data is indeed long-lat coordinates, but I wanted to use UTM as my understanding is that this is more accurate when dealing with distances and areas.
After setting the CRS of the raw data to long-lat, I converted it to UTM prior to running my home range code and converting the ...
1
You need to export the dates variable as well because it is a global
beginCluster()
z1 <- clusterR(ndvi, calc, args=list(fun=f), export=c('a','dates'))
endCluster()
or
c1 <- getCluster()
clusterExport(c1, c("a","dates"))
z2 <- calc(ndvi, fun=function(x){ t(parApply(c1, x, 1, f))} )
0
That should be as easy as:
filename="vourlionas"
setwd("mydir")
d = read.table(paste(filename,"txt",sep="."),
header = FALSE,
sep="",
col.names=c("x", "y"),
fill=FALSE,
strip.white=TRUE,
stringsAsFactors = FALSE)
some data manipulation here
writeOGR(obj=test1, dsn="mydir", layer = ...
0
This data looks like lat-long coordinates:
data<-data.frame(x=c(-50.3, -49.9, -50.0, -50.6, -55.3, -55.4, -55.5, -55.3, -54.9, -54.4, -51.5, -51.2, -50.8, -50.3),y=c(50.3, 48.8, 48.1, 47.4, 48.2, 47.4, 50.1, 48.1, 47.5, 50.7,50.4, 50.7, 50.5, 48.3))
but you set its CRS to a UTM 21 region:
#set CRS
proj4string(data)<-CRS("+proj=utm +zone=21 +ellps=...
2
If the focal value is not specifically indexed in the function it looks like the focal function is assigning the incorrect values and positions back to the matrix. However, this is expected behavior. For these types of operations, you have to specifically index the focal value to get the correct results. One would think that the focal value would ...
0
Looks like your raster datatype may be set to integer.
Try: rasterOptions(datatype='FLT4S')
0
Don't do the join. st_nearest_feature(a,b) will get you the index (row number) of the nearest feature in b to each feature in a.
EG using data p and l from ?st_nearest_feature made into sf data frame:
> (nearest = st_nearest_feature(p,l))
[1] 1 2 2 3
Then use st_distance to get the element-wise distances between each element of p and the corresponding ...
0
I would recommend using the lhs package documentation for reference
You would need to install the released version of lhs from CRAN
install.packages("lhs")
or via devtools
devtools::install_github("bertcarnell/lhs")
in your case a example to create a random LHS with 30 samples and 3 variables would be:
X <- randomLHS(30, 3)
1
A little more information would prove useful to help answer your question. What language are you using, and what function are you using to import your data? I will base my answer on the assumption you are using R with the raster library and the raster function to import your .tif files. If you want your tif files to come in as a stack, use the raster::stack()...
1
I did not manage to solve why grid_metrics fails to create the .vrt file when applied to a big catalog but I managed to be able to create the same output. I solved my problem by usig the gdalUtils tools from QGIS instead of Anaconda as the error said "CreateProcess' failed to run 'C:\PROGRA~3\ANACON~1\Library\bin\gdalbuildvrt.exe" and also the ...
1
If you have na.rm=TRUE then movingFun will use a smaller window at the edges:
> movingFun(1:12, 3, max,na.rm=TRUE, type="around")
[1] 2 3 4 5 6 7 8 9 10 11 12 12
So the two 12s at the end are the max of 10,11,12 at the 11th location, and the max of 11,12,NA at the 12th location.
If you don't want those rasters in your output, subset them off.
2
Best thing to do with questions here is to make a representative simple example.
Here's a 3x4 raster with some values:
r = raster(matrix(c(1,1,1,2,1,3,55,12,3,3,2,1),3,4))
Now build a data frame - find the unique ID values:
ids = unique(r[])
ids
# [1] 1 2 55 3 12
And make a data frame like yours with a couple of data columns:
d = data.frame(ID=...
0
Fast without fancy programming: gContains from package rgeos
download.file("http://gis.ices.dk/shapefiles/ICES_ecoregions.zip",
destfile = "ICES_ecoregions.zip")
unzip("ICES_ecoregions.zip")
Read eco region shapefiles
ices_eco <- rgdal::readOGR(".", "ICES_ecoregions_20171207_erase_ESRI", verbose=FALSE)
Make a large data.frame (361,722 ...
0
Sometimes parallel is not the way to go. If the function is not written for it you get all the overheads and none of the benefits.
Often, just reducing the demand for memory allocation solves performance issues.
I have found that if the number of points is high it is much quicker to hit the raster layers one at a time (or in smaller stacks) rather than ...
1
You could use exactextractr::exact_extract for this. It will work with multi-part polygon inputs.
library(raster)
library(sf)
library(exactextractr)
semana21 <- raster("Semana21.tif")
poly <- st_read("shapefile.shp")
poly$sum_semana21 <- exact_extract(semana21, poly, 'sum')
Note that these results will differ from raster::extract because the ...
2
You can get a table of cell values, coverage fractions, and center coordinates with the include_xy argument to exactextractr::exact_extract. Here's an example:
library(raster)
library(sf)
library(exactextractr)
# Pull municipal boundaries for Brazil
brazil <- st_as_sf(getData('GADM', country='BRA', level=2))[1:10, ]
# Pull gridded precipitation data
...
2
This question have already been answered (partially) here. Individual tree segmentation in broad coverages is not currently (v2.2.2) an easy task in lidR. There are several reasons for that:
lidR brings 3 algorithms for making a CHM which each have various tweaks and various parameters + 4 algorithms for ITS including a point-cloud based methods which is ...
1
I don't think this is possible (yet). GetFeatureInfo requests have not been integrated in r leaflet. This was an issue in the leaflet.extras package, but not (yet) implemented. See https://github.com/bhaskarvk/leaflet.extras/issues/84
3
you can use the raster::xyFromCell function to get a matrix of coordinates from a vector of cell numbers, which you can cbind onto the extract (as the xyFromCell result is in exactly the same order as the vector provided)
I've used sf but of course this works with readOGR etc.;
require(raster)
require(sf)
# make 2 dummy rasters for a dummy stack
r <- ...
1
You can easily modify an sf object using a matrix instead of a numeric vector, whether something you've built directly or loaded from a geopackage or shapefile. The only difference is you have to declare the scale or transform factor as a matrix .
# Using your example but not assigning WGS84
p1 <- rbind(c(1000,1000), c(1000,2000), c(2000,2000), c(1000,...
1
In your example, just multiply earlier in the polygon construction;
library(sf)
p <- st_polygon(list(rbind(c(1000,1000), c(1000,2000),
c(2000,2000), c(1000,1000)),
rbind(c(500,500), c(500,100),
c(1000,1000), c(500,500))))
p.shift <- st_sfc(p * c(2,3))
If you want to do this on more complex data you ...
2
I'm dealing with the same issue here.
The 'problem' is originated because the library comes from the GEOS library, which considers that coordinates in a two-dimensional, flat, Euclidian space.
BUT....For longitude latitude data, this is not the case. For more info on this, please check... https://cran.r-project.org/web/packages/sf/vignettes/sf6.html#...
1
outlist is a list. The output of catalog_apply() is thus a list with one element per chunk and each element is a list. The engine is not able to merge such complex output with nested lists. So automerge option fails and the output is the regular, unmerged, list + a warning to tell you that automerge did nothing.
You are trying compute mean elevation for ...
2
After plotting the the points you need to call text():
plot(elev)
plot(pts_dotcoords, add = TRUE, pch = 16)
text(pts_dotcoords, labels = c("locA", "locB"), pos = 4, offset = 0.7)
See also ?text:
pos: a position specifier for the text. If specified this
overrides any adj value given. Values of ‘1’, ‘2’, ‘3’ and
‘4’, respectively ...
0
Edit
Version 0.9.0 of the Resource Sharing plugin was made available on http://plugins.qgis.org today (https://plugins.qgis.org/plugins/qgis_resource_sharing), and the good news is that it includes support for R script collections.
See https://gis.stackexchange.com/a/341008/22646 for more details, and for information about available collections and ...
1
You need to test if any of the park NAME attributes are in the list of reserves:
> natparks_sub <- natparks[natparks$NAME %in% marinereserves,,drop=TRUE]
> dim(natparks_sub)
[1] 7 29
Using == here tests each element of NAME with each element of marinereserves in turn. You would have got a warning if not for the chance that you have 7 marine ...
3
This seems to be Zone 4 rather than Zone 3. If I do:
> proj4string(coord_gk) <- CRS("+init=epsg:31468")
I get:
> spTransform(coord_gk,CRS("+init=epsg:4326"))
SpatialPoints:
GK_R GK_H
[1,] 12.7424 48.2167
Coordinate Reference System (CRS) arguments: +init=epsg:4326
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
which ...
3
You can use the dissolve argument, which dissolves all features with same value into polygon (or multi-polygon)
library(raster)
test <- raster("sinks.tif")
regions <- clump(test)
regions_shp <- rasterToPolygons(regions, na.rm = T, dissolve = T)
regions_shp
result of regions_shp:
class : SpatialPolygonsDataFrame
features : 887
extent ...
0
Thanks @ElioDiaz, @Spacedman - here is code (could be simplified with chained operation) based on your suggestion that does exactly what I intended with trying to emulate result of union tool in ArcMap.
outline <- read_sf('outline.shp') # black circle in picture
other_polys <- read_sf('other_polys.shp') # purple polygons
within <- ...
1
The nb object is a ragged list, so you can either get them as elements of the list:
> nb[[1]]
[1] 0
> nb[[2]]
[1] 3 5 6
> nb[[3]]
[1] 2 4 5 6 7
showing the first polygon has no neighbours, the second one connects to 3, 5, and 6, and so on.
Or convert to a full binary adjacency matrix:
> mm = nb2mat(nb,zero.policy=TRUE, style="B")
> mm
[,...
Top 50 recent answers are included
