4
It's faster and easier to use st_intersects, no need to loop. The output is a bit obscure, essentially a classed list of feature IDs intersected, so we get the lengths() which is the number of points inside each feature.
First polygon has six points in it, in this example.
library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
poly <- ...
3
You were almost there. You missed to compute the pulseID with laspulse() and you missed that the scan angle is stored in ScanAngleRank
Since lidR 2.0.0 pulseID is no longer computed at read time. And since rlas 1.3.0 that introduced support of LAS 1.4 format the attribute ScanAngle is now ScanAngleRank. The name ScanAngle is reserved for angle stored in LAS ...
3
In your case the filter you are using is a simple one: Classification != 2. And you don't need the ground points at all. You are better to use a streaming filter and a streamed processing.
ctg <- catalog("/...2015TestGroup")
opt_chunk_size(ctg) <- 0
opt_chunk_buffer(ctg) <- 0
opt_output_files(ctg) <- ".../Outputs/2015nonground/{XLEFT}_{YBOTTOM}...
2
Working with an outline of Pakistan from the GADM data:
p = raster::getData("GADM",country="PAK",level=0)
p is a SpatialPolygons data frame in lat-long coordinates. You should convert this to a proper equal-area projection but EPSG:3857 should be close enough:
pp = sp::spTransform(p, "+init=epsg:3857")
Now chop using the code you pasted in your question. ...
2
Seems like you are passing a list to the read function instead of a single file
When you call
shp_paste<- paste(input_path, ifile, sep = "")
the ifile object can be a list of files, so shp_paste has a length greater than 1.
You might change this by changing the code to
shp_paste<- paste(input_path, f, sep = "")
since you are iterating over files....
2
lasfilter will only directly take objects of class LAS (and not LAScatalog) as per the package documentation.
One way to go is with catalog_apply:
This function gives users access to the LAScatalog processing engine. It allows the application of a user-defined routine over an entire catalog.
So, embed lasfilter within a user-defined function and pass it ...
2
Here is a very short sample of code that shows the proj4 library is functioning and your proj4string is correct:
library(proj4)
proj4string <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
# Source data
xy <- data.frame(x=1001557, y=217404)
# ...
2
This is not possible in the RSAGA package without a massive restructuring.
RSAGA works by running an external saga_cmd command, which starts a new process that does input and output via files.
To do otherwise could be done, but it would require dynamically linking to the SAGA C++ code, and writing R to C++ wrappers (which would probably be best done using ...
2
Your question is related to LAScatalog processing engine tuning. A topic not documented in the official documentation. The only one existing documentation at the time being (june 2019) is a wiki page that provide an example to change the drivers.
In short the drivers used to write objects to files are stored in the LAScatalog object. You can access to them ...
2
You're doing some unnecessary work. If you want a stack then just call rstk <- stack(files). If the data doesn't share extents then call agg <- sapply(files, raster). Some people always use stack--if that's you're intent fine--but just don't make a function that only calls another function.
It's not clear what your intent is with putting the ...
2
In R there is package cleangeo for geometry checking and celaning. You can use for example:
library(cleangeo)
report <- clgeo_CollectionReport(spatial_object)
report
#returns a table:
type valid issue_type
1 <NA> TRUE <NA>
2 rgeos_validity FALSE GEOM_VALIDITY
3 rgeos_error FALSE ORPHANED_HOLE
...
2
Found a solution using sf package:
So having a matrix of points, First I need sf object
pts_sf = sf::st_as_sf(as.data.frame(pts), coords = c(1,2))
Next I can perform buffer
pts_buf <- sf::st_buffer(pts_sf, 2)
then union
pts_com <- sf::st_union(pts_buf)
and ath end I can get normal polygon from it
pts_pol <- sf::st_cast(pts_com, "POLYGON")
...
1
ratify is the right option, but you should do an extra step. You need to create a dictionary to store desired values and create a numeric column to be used in rasterize process:
library('raster')
library('rgdal')
# Load a SpatialPolygonsDataFrame example (Brazil administrative level 2) shapefile
dat <- raster::getData(country = "BRA", level = 2)
# get ...
1
After posting this I realized that it cannot correctly convert the points when I'm not telling it which columns contain the points I need it to convert.
Here's a version that seems to work just fine:
####libraries####
library(proj4)
library(tidyverse)
library(leaflet)
####working directory####
wd <- ("MY WD")
setwd(wd)
####file location####
file <...
1
You are having an issue with only a small subset of your points actually intersecting your polygon grid as well as the need for a projection transformation from geographic to Mercator.
It is difficult to evaluate the results because all the attributes associated with the grid, that occur outside the intersection, will be nodata (NA). In just looking at the ...
1
Let's take a small sample of the house data for illustration:
> set.seed(123)
> h = house[sample(25357,100),]
Then you can use the deldir function to construct a Delaunay Triangulation:
> library(deldir) # install from CRAN if not got already
> hd = deldir(data.frame(coordinates(h)))
One element of this looks like this:
> head(hd$delsgs)
...
1
Your circles are okay, it seems the raster is in the wrong position.
Its bbox in its coordinates is:
> bbox(raster)
min max
s1 2101250 2318250
s2 315500 566500
and these are supposedly in this coordinate system:
> projection(raster)
[1] "+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25 +a=60 +rf=6378137 +lat_2=45"
and if you ...
1
Considering a case for instance the original raster has 30m x 30m resolution (pixel size), and we want it be 100m x 100m.
If we can find a common divisor between before-after resolutions (10m in this case) then;
disaggregate() to refine 30m grid to 3 x 10m grid, then;
aggregate() to upscale 10m grid to 100m grid.
And the code will be:
library(raster)
# ...
1
The data has coordinates in both geographical and projected CRS.
EPSG:4326 / LONGITUDE, LATITUDE
ESRI:102671 EPSG:3435 / X COORDINATE, Y COORDINATE
You have chosen LONGITUDE and LATITUDE as input, then just select EPSG:4326.
Try;
tif_data <- read.csv("~/Documents/code/ChicagoPackage/data/TifProjects.csv")
tif_data <- tif_data[!is.na(tif_data$...
1
You are converting from sf to sp in the first, and from from sp to sf in the second - you should avoid timing those conversions.
sf is sometimes slow with points because of the way they are stored, but what you gain is far greater consistency than with sp, and usually faster ops.
Here I think it is comparable, but will depend on your actual data.
Here'...
1
Tile services typically (not always, but almost always) use the Web Mercator projection (3857), because it has nice properties for tiling. It's a safe assumption.
In this case, you don't need to assume, since it's in the metadata for the service: https://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/
This tile pattern (z/x/y) is ...
1
If you row-bind together the list that you get when reading GeoJSON feature strings with read_sf together you get an sf data frame:
> do.call(rbind, lapply(d$geom,read_sf))
Simple feature collection with 3 features and 0 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
...
1
sapply returns a list, so in your case its returning a list of stacks. This code demonstrates your situation and reproduces your error:
> r = raster(matrix(1:12,3,4))
> s = stack(r,r,r)
> all_s = list(s,s,s,s)
> as.data.frame(all_s)
Error in as.data.frame.default(x[[i]], optional = TRUE) :
cannot coerce class "structure("RasterStack", package ...
1
If you read the documentation for spCircle, specifically the Value section:
Value:
A list with the following components...
spCircle : The "‘SpatialPolygons’" polygon object.
location : The "‘SpatialPoints’" point object.
you'll see it returns a list with two components. It looks like you want the polygon object. In your loop make a ...
1
Read the help for SpatialPoints, specifically the arguments:
Arguments:
coords: numeric matrix or data.frame with coordinates (each row is a
point); in case of SpatialPointsDataFrame an object of class
SpatialPoints-class is also allowed
So a simple two-column matrix (here with four rows) can be a set of four points:
> ...
1
You can use aggregate to summarise the table then merge the results back to the attribute table.
If you've read in 'd' using the sf package (or rgdal, though I haven't tested that):
agg <- aggregate(formula = Height ~ Location.ID, data = d, FUN = mean)
d <- merge(d, agg, by = "Location.ID")
1
Create a test data set to see what is going on:
> r1 = raster(matrix(1:6,2,3))
> s1 = stack(list(r1,r1*2,r1*3))
Convert to data frame this way:
> as.data.frame(s1)
layer.1 layer.2 layer.3
1 1 2 3
2 3 6 9
3 5 10 15
4 2 4 6
5 4 8 12
6 6 12 18
...
1
How about using *_join from dplyr package, for instance:
map_and_data <- dplyr::left_join(mymap, mydata, by='Country')
Here is one example you can test.
library(sf)
nc = st_read(system.file("shape/nc.shp", package="sf"))
my_table = data.frame(NAME= nc$NAME, INCOME= runif(100, 5000, 10000))
library(dplyr)
merged = left_join(nc, my_table, by= "NAME")
...
1
In order to answer, let’s put aside important, but broad issues:
The fact that identifying and segmenting trees is a very complex analysis which depends on many things (things related to the type of vegetation, and quality and amount of available data, for example).
That processing 'large point clouds' in R is a real concern (due to memory limitation), and ...
1
It is difficult to understand the memory problems you report as you do not show the code that causes it. Perhaps you do something wrong. It could also be useful to see the results of
show(big.raster)
and
canProcessInMemory(big.raster, 4, TRUE)
(this would look something like this)
#memory stats in GB
#mem available: 53.67
# 60% : 32.2
#mem ...
Only top voted, non community-wiki answers of a minimum length are eligible
