New answers tagged r
3
In the best world the memory allocated to load a las file should be almost equal to the size of the file. But this is not possible especially in R. This is explained in this vignette.
It is important to note that R only enables manipulation of 32-bit integers and 64-bit decimal numbers. But the las specification states, for example, that the intensity is ...
1
R is reporting the size of the object in memory, not the size of the file.
If you run the example from readLAS you'll get las which is:
> las
class : LAS (LASF v1.2)
point format : 1
memory : 5.6 Mb
this size depends on which parts of the file you read in (via the select and other options) and how much the LAS file is compressed.
Most ...
0
You may consider using the R package hsdar as a starting point to carry out transformation of spectral response data as well calculation of vegetation indices. You will find additional information about this package in this article published in the Journal of Statistical Software. While the hyperSpec package might also be of use.
This tutorial provides ...
0
I can now duplicate your error with a test package I've made:
> d = st_read("/tmp/test.gpkg",query="select resolution")
Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
Query execution failed, cannot open layer.
In addition: Warning message:
In CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
GDAL Error 1: In ...
1
This was a bug fixed in version 2.1.4. In short in R when you subset a matrix you get a matrix (that makes sense) but when you subset a single row you get a vector instead of a 1 x n matrix. R's behavior is inconsistent. In rare cases you may have a single remaining triangle in pitfree(). This is why you got an internal error not a matrix. This case is now ...
1
gstat can do spatial-temporal kriging, as explained here: https://www.r-bloggers.com/spatio-temporal-kriging-in-r/ . Time is just a third dimension, and you could use "elevation" data instead of time, in the same algorithms, reaching the desired output.
As an alternative, although not completely native in R, Saga Gis 7.4 just implemented full 3d kriging. ...
4
Take a look at some of the identity functions in rgeos such as gTouches, gIntersects, gContains, gRelate, gWithin, ect...
At their most basic, these functions can return a Boolean that will allow you to perform the equivalent of a spatial select. To control output, take a close look at the function(s) arguments and experiment a bit with outputs from ...
7
Reprojecting rasters is usually a bad thing to do. It involves a non-reversible transformation from one grid system to another grid system that can have a non-linear relationship to the first. Hence the value in a cell of the new system can end up being some average of whichever grid cells in the source raster it overlapped.
If you have a raster and points ...
2
Lambert equal-area projection
laea_centered <- "+proj=laea +lat_0=-70.30744 +lon_0=43.15268"
call me crazy, but don't you have the lat and long switched around?
4
Yes, you can just loop through the individual features and save them like so:
library (sf)
dat <- read_sf ("~/data.shp")
for (i in seq_len (nrow (dat)))
{
fname <- paste0 ("~/split_", i, ".shp")
write_sf (dat [i, ], fname)
}
1
Mask raster function in R not working to extract values in one raster based on conditions in another
You have a mask with 0, 1, and NA values and you want to keep only the 1s. Hence you need to convert 0 and NA to the same value.
Here's a bit of a trick. Diving m by itself results in NA where m is 0 or NA (since 0/0 returns NaN, which acts like NA in mask) and 1 where its 1:
rmask = mask(r,m/m)
not sure if that's any quicker than the other solutions.
1
Mask raster function in R not working to extract values in one raster based on conditions in another
First change the values of 0 in wintPCP to NA:
wintPCP[wintPCP == 0] <- NA
then mask:
masked_ras <- mask(MAT2resampled, wintPCP)
1
An error in a factor of about ten in area corresponds to a scale error of about 3.16 in length (square root of ten), which I have seen happen when people use imperial foot units instead of metres (or yards for really old projection systems).
Oh look... https://epsg.io/2314
Remarks: Foot version of Trinidad 1903 / Trinidad Grid (code 30200) used by some US-...
1
Simple solution using the 'sf' package
library(sf)
library(ggplot2)
library(mapview)
library(lwgeom)
library(rnaturalearth)
# world data
world <- rnaturalearth::ne_countries(scale = 'small', returnclass = 'sf')
# Fix polygons so they don't get cut in ortho projection
world <- st_cast(world, 'MULTILINESTRING') %>%
st_cast('LINESTRING', do_split=...
2
First convert the rasters to variables in the same dataframe, then calculate the pairwise correlations and use the package 'corrplot' to display the results in a matrix.
#dependencies
library(raster)
library(corrplot)
#read in rasters
r1 <- raster("IMG_0003_1.tif")
r2 <- raster("IMG_0003_2.tif")
r3 <- raster("IMG_0003_3.tif")
r4 <- raster("...
2
The result is correct. Since the extent of A is larger, the function is correctly cropping the extent of B. To get the same "shape", try something along the lines of:
mask(B, crop(B, extent(A)), A)
However, this will still show a smaller extent in B. If you really need the two rasters to show the same extent of data then you will likely have to coerce ...
1
In QGIS you could use the "Extract layer extent" processing tool to get a bounding box for Raster B. Then use that as the mask layer in the "Clip Raster by Mask Layer" tool to clip Raster A
0
I wasn't originally going to answer my own question but a colleague of mine (who does not use this site), wrote me a bunch of python code to do what I am after; including limiting the cells to have the distance to coast for only terrestrial cells and leaving the sea based cells as NAs. The following code should be able to run from any python console, with ...
1
OSM data has a lot of different data types, tags, etc. It is probably best as a lot of different shape files. You could import it into a Spatialite file using either OGR2OGR or QGIS, and then extract each layer type that makes sense for your application. Eg. one for roads, one for roads, one for linear water features, another for polygon water features.
(as ...
1
Here's a reproducible example. Get county shapes from USAboundaries:
> library(USAboundaries)
> data = us_counties()
> dim(data)
[1] 3220 13
add a category column with about 300 in each category, and a varinterest
> data$category=paste0("cat",sample(1:10,nrow(data),TRUE))
> table(data$category)
cat1 cat10 cat2 cat3 cat4 cat5 cat6 ...
3
First you need a map of world polygons.
> library(rworldmap)
> m = getMap()
Then you need your points in a spatial format. I'll make some random points on the globe and create a SpatialPointsDataFrame using the sp package and making sure they have the same projection string as the map:
> set.seed(123); d = data.frame(long=runif(100,-180,180), lat=...
0
library(raster)
shapefile(tornados, filename='path/to/file.shp')
Make sure you have (or have given yourself) permissions to write to the directory.
2
I'd say this was a bug in lasmergespatial that occurs when the raster is not stored in memory.
The code eventually ends up here:
lidR:::lasmergeRasterLayer = function (las, raster)
{
cells <- raster::cellFromXY(raster, coordinates(las))
return(raster@data@values[cells])
}
but extracting the cell values like that only works for in-memory ...
2
This is my settings.
The location R is installed.
R scripts folder is accessed from Menu > Settings > User Profiles > Open Active Profile Folder.
User library folder is where I have R packages.
Then, RSX script for Ripley's K:
##Kazuhito Vector Analysis=group
##Layer=vector point
##showplots
library("maptools")
library("spatstat")
points=as.ppp(...
0
After I couldn't get @Spacedman's suggestion working I dropped it, feeling it was a bit out of hand; I later stumbled the simplest answer: I misunderstood load() and get().
After the file <- load(c:/data/regionX-yearY.rda) I had the filepath; this mislead me a bit as it's a ponderous load for these .rda files, actually "feels" like R has loaded something ...
0
If you have longitude, latitude and ellipsoidal height, you have coordinates over an ellipsoid (the WGS84 ellipsoid in your case) and related to a coordinates frame (the WSG84 geodesic frame in your case). Imagine that, a single ellipsoid among all possible, with its center, axes (ellipsoid axes) and origin of coordinates (center meridian) positioned and ...
0
You can also just get at the cell frequencies and proportions. If in a projected coordinate system then you can get the class areas by multiplying resolution by cell counts.
library(raster)
r <- raster(matrix(round(runif(25*25, 1, 4)),25,25))
( f <- freq(r) )
( p <- data.frame(f, p=f[,2] / sum(f[,2])) )
0
This can be done using dplyr with something like the following:
library(dplyr)
library(raster)
r = raster()
r[] = sample(-5:-2, ncell(r), replace=TRUE)
as.data.frame(r) %>%
group_by(layer) %>%
tally() %>%
mutate(area = n * res(r)[1] * res(r)[2])
# A tibble: 4 x 3
layer n area
<int> <int> <dbl>
1 -5 16260 16260
...
0
The ellipsoidal height does not affect the latitude longitude. so only extracting xy from xyz should be fine.
0
For this case, I would use min value between both rasters using overlay() function:
(reproducible example)
library(raster)
r <- raster()
r1 <- setValues(r,3)
r2 <- setValues(r,3)
r1[1:20000] <- 1
r2[44800:64800] <- 2
plot(stack(r1,r2), col=terrain.colors(255))
Minimum between both rasters:
r3 <- overlay(r1,r2,fun=min)
plot(r3)
If ...
0
Try to separate which are not NA and which are value 3
summer_indices <- which(propSummPCPmod4[] == 1,3)
summer_notNA <- which(!is.na(propSummPCPmod4[]))
newraster[summer_notNA] <- 3
winter_indices <- which(propWintPCPmod4[] == 2,3)
winter_notNA <- which(!is.na(propWintPCPmod4[]))
newraster[winter_notNA] <- 3
newraster[...
0
I solved the issue using a nested ifelse() construction:
myform<-function(PhenoDate,dPrecJan,dPrecFeb,[...],dPrecDec) { ifelse (
PhenoDate < 32, PhenoDate * dPrecJan, ifelse (
PhenoDate>=32 & PhenoDate<60, 31 * dPrecJan + ((PhenoDate-31) * dPrecFeb), ifelse (
PhenoDate>=60 & PhenoDate<91, 31 * dPrecJan+ 28 * ...
0
I had the same problem. You should install and the ncdf4 package.
install.packages("ncdf4")
Then be sure to run library(ncdf4) before running the brick function from raster. Below is my code:
library(raster)
library(ncdf4)
r <- brick("E:\\test.nc",varname="precip")
r
0
There seems to be nothing wrong with the file and I can't reproduce the error except by trying to access a file that doesn't exist. I have the file in the same folder that I start R from (ie ./) and get this:
> file.size("./liberia_concessions.tif")
[1] 7747782
> lc = raster::raster("./liberia_concessions.tif")
> lc
class : RasterLayer
band ...
1
Thanks for the responses. Indeed, the issue was with the CRS. The native CRS for the data I used was WGS 84 lat-long. The CRS of the mask I was clipping to was UTM14. As the CRS did not share the same extend, all of the values from the point data were excluded. When I imported the point data, I first set the CRS to UTM14. However, this transformation did not ...
1
This seems to be a bug that happens when there's less than 8 layers in your stack. Here's a test function that generates and runs on an ixjxk rows x columns x layers stack:
testkendall <-
function(i,j,k){
raster_stack = do.call(stack, replicate(k, raster(matrix(runif(i*j),i,j))))
raster.kendall(raster_stack)
}
And it works on 23x17x8:
> ...
2
Within QGIS you could try the GRASS plugin. As far as I know it manages the memory better than R, and I expect the other solution to fail on large areas.
the GRASS command is called r.grow.distance , which you can find in the processing toolbar. Note that you need to convert your line to raster at first.
One of your issue could be the size of the output, ...
1
Many options here:
Use LAStools
lidR is good for some tasks but weaker for some other. This is a typical case where las2las from LAStools is more suitable for what you want to achieve.
Use streaming filter
If you want to work within R anyway, do not load the full point cloud in R. Read only the points of interest. Your filter is very simple here.
Rplot &...
6
May be a solution to try:
Generate a grid (Type "point", algorithm "Create grid")
Calculate the nearest distance between your points (grid) and your line (coast) with the algorithm "join attribute by nearest". Be carefull to choose only a maximum of 1 nearest neighbors.
Now you should have a new point layer with the distance to the coast like in this ...
5
If you are trying to create a 2.5D digital terrain model (a height field) from a point cloud then there might be two possible problems. Multiple points in the cloud at the same location don't add any information to the data, so are dropped. Multiple points with the same XY location but different Z location (height) are inconsistent with the idea of a planar ...
1
Your example was not reproduceble so I tryed to create a small set imitating your data.
If I understood right, ou can do:
`library(raster)
propSummPCPmod4 <- raster(matrix(c(1,3,3,1,3,3,1,1,NA,NA,NA,NA,NA,NA,NA,NA),4,4))
propWintPCPmod4 <- raster(matrix(c(NA,NA,NA,NA,NA,NA,NA, NA,2,3,3,2,2,2,2,3),4,4))
newraster <- raster(matrix(c(NA,NA,NA,NA,NA,...
0
Here is an example, with 100 polygons and 150 points, using similar syntax of st_join. I'm never getting more than one match when using k=1 (there can be zero matches when the nearest feature is further than maxdist).
If you can please post a reproducible example where you are getting more than one match will be happy to look into it.
library(sf)
## ...
0
As commented, my rasterstack of 597 ndvi images is:
chitral_ROI_4
create a regular daily time series object by combining data and date information
(s <- bfastts(as.vector(chitral_ROI_4), dates, type = c("irregular")))
s.d.linear <- round(na.approx(s),0)
Create time series to data frame function
time.series.to.dataframe <- function(...
0
An additional answer for projected data. bfastSpatial::areaSieve() takes a threshold argument, such that clumps smaller than the threshold are removed.
https://www.rdocumentation.org/packages/bfastSpatial/versions/0.6.2/topics/areaSieve
0
Indeed, you have objects with different projection strings.
> identicalCRS(shape, point.rep)
[1] FALSE
Inspecting them makes it clearer:
> shape@proj4string
CRS arguments:
+proj=utm +zone=22 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
> point.rep@proj4string
CRS arguments:
+init=epsg:32722 +proj=utm +zone=22 +south +datum=WGS84 +...
9
With PyQGIS and GDAL python library is not very difficult to do that. You need geo transform parameters (top left x, x pixel resolution, rotation, top left y, rotation, n-s pixel resolution) and rows and columns number for creating resulting raster. For calculating the distance to the nearest coastline, it is necessary a vector layer for representing ...
0
I would try other way around. If you are using poligon of NZ then convert polygon edges to line. After that create buffer on the boundary for every 25 meters of distance from boundary(maybe centorid might help in determing when to stop). Then cut buffers out with polygon and then convert that polygons to raster.
I am not sure this would work but definitely ...
1
Answer as follows:
library(geosphere)
df_sp <- as_Spatial(df)
destPoint(df_sp, df_sp$distance, df_sp$bearing) %>%
as.data.frame() %>%
st_as_sf(coords = c('lon', 'lat')) -> df_sf
0
I used BFAST for my thesis some years ago, and this is the code I used. Given the age of this code, things may have changed significantly, to make specific elements not function entirely any more, but the process of creating a function to apply BFAST and then applying that function on the raster is still a good way to approach the problem. Do note that the ...
3
By design, the symbols in the legend are as close to the symbols in the map as possible.
If you want the legend to be different, one approach is to create a manual legend:
data(World)
tm_shape(World) +
tm_fill("black") +
tm_dots("footprint", size = .1, legend.show = FALSE) +
tm_add_legend(type = "fill",
col = c(tmaptools::...
Top 50 recent answers are included
