New answers tagged r
0
You can use the function crosstable in the package raster.
Assuming that you have the maps 2014, 2018 saved as raster landscape_2014 and landscape_2018, the code would be like that:
contingencyTable <- raster::crosstab(landscape_2014, landscape_2018, long = TRUE)
The argument long define if the results will be in data.frame or table format.
0
Just a short remark to the extract function from raster package when you work with lat and lon values - use SpatialPoints:
x<- 44.8386
y<- 0.5783
value<-extract(raster, SpatialPoints(cbind(x,y)))
Without it, you may just receive NA values.
1
This should get you started with reading in all of your feature classes, then you can perform a similar operation to what you've included in your question.
library(sf)
path<-"C:/Users/GIS Data/Some Folder/your.gdb"
gdb<-st_layers(path)
list_of_features<-purrr::map(gdb$name,~st_read(dsn=path,layer=.))
st_layers will list all of the names within ...
2
For a non memory-safe option you can efficiently use an apply type function on an sp SpatialPointsDataFrame object. This has the advantage of not having to use which, which replicates a vector the size of your raster (which could be huge), is much faster than a for loop and directly results in an sp point object.
Add libraries and create data
library(sp)...
0
check the data type,
class(shape)
if it returns sf "data.frame" you should use st_transform instead of spTransform
shape_proj<-st_transform(shape, CRS("+proj=gnom +lat_0=90 +lon_0=-50"))
1
I have a similar problem: I need to snap a point on a linestring. Even though I don't split the linestrings, I hope my solution would be useful for you.
My approach is this: Due to the rounding errors, the nearest point needs not to be precisely on a linestring. My solution is to add a new point to the linestring that is precisely the nearest point. (...
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")
...
0
While writing this up I had the idea ... perhaps just remove the alpha layer using the dropLayer() function? I do this using
Brick <- dropLayer(Brick,4)
And now tm_rgb() works fine, albeit it is significantly slower than plotRGB(). Perhaps there is a better solution out there that can speed things up?
0
Always make a tiny example before trying something on 325 layers of a zillion pixels:
> r1 = raster(matrix(seq(32,212,len=12),3,4))
> r2 = raster(matrix(sample(seq(32,212,len=12)),3,4))
> r3 = raster(matrix(sample(seq(32,212,len=12)),3,4))
> s = stack(r1,r2,r3)
> plot(s)
Then you can do arithmetic on raster stacks and it works on a pixel-by-...
-1
I found this code which might help other people as well:
https://github.com/environmentalinformatics-marburg/heavyRain/commit/30edea7e11a1fd41f98ed397fea32c4376560891#diff-f4c1e81d0b0b30884c2c20442238b2e4
1
Assuming you want to save to temporary shapefile the last line could be something like:
shapefile(shpFile, tempfile(pattern = "", fileext = ".shp"))
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 ...
0
I'd usually load my data into an sp object and use rgdal to project:
library(sp)
library(rgdal)
dfr <- read.csv('~/Downloads/NYPD_Complaint_Data_Current__Year_To_Date_.csv')
dfr <- with(dfr, dfr[!is.na(X_COORD_CD) & !is.na(Y_COORD_CD), ]
xy <- dfr[, c('X_COORD_CD', 'Y_COORD_CD')]
prj <- CRS("+proj=lcc +lat_1=40.66666666666666 +lat_2=41....
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 ...
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 <...
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)
# ...
0
Here's a solution using purrr
library(tidyverse)
library(sf)
df <- tibble::tribble(~ID, ~X, ~Y, ~prevXval, ~prevYval,
1, -0.121, 65.001, 0.067, 65.117,
2, 180.000, -50.039, -179.879, -50.156,
3, -0.075, 47.787, 0.038, 47.904,
4, ...
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 ...
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
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 ...
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
...
0
You could also do that by creating a virtual raster file (.vrt) in QGIS and editing it manually in any text editor.
All classes and colors will be interpreted from 0 on, thus I include NoData here as a hint to that.
<VRTDataset rasterXSize="109300" rasterYSize="89999">
<VRTRasterBand dataType="Byte" band="1">
...
<ColorInterp>...
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")
0
The solution to this was quite simple - I was able to open the .gdb file in QGIS and use the 'Select by Expression' function to select the fields I wanted. The formula I used was based on this post:
SCINAME IN (
'kent',
'essex',
'sussex',
'london'
)
I was then able to save the selection as a shapefile.
0
You are barely at the start of this analysis. As such, I would highly recommend trying to track down any resources as a starting point to learn the basics of spatial analysis in R. Honestly, Google is a good start and it is clear that you did not search this StackExchange very well because the fundamentals of your task have been covered numerous times herein....
0
Here's code that follows what you are doing.
First make a raster like your source:
source = raster(nrow = 467, ncol = 805,
xmn=-119, xmx=-85.45833, ymn=13.54167, ymx=33)
crs(source) = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
Fill it with zeroes plus a diagonal line of 99s so we can test alignments:
source[] = 0
source[seq(1,...
1
Please put the libraries in your code, or use ::; it's time consuming to have to infer what packages might be in your namespace.
Your data is irregular points. If you just want them in a dataframe, then covert them to sf or sp objects. Your data values are also nominal, so unless you want to calculate proportions there's no sensible way to convert that to ...
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
...
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 ...
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 ...
0
The argument to set this is legend.outside.size. The default depends on which style is used, but for the default style, the value is 0.3 (check tmap_options()$legend.outside.size).
Setting this to a smaller 0.1 value will make the legend take less space. Note however that the font of the legend is also reduced though, so there's a trade-off in the ideal ...
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
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)
...
0
Do your layers A and C just have NA where they don't have a value?
i'd suggest using raster::distance. Just run separately for rasters A and C and merge (if relevant)
Here's a dummy for you, ncells = 1 million, using an extremely crude boundary line;
require(raster)
# make 2 rasters (your A and C), with generic values
r1 <- raster(xmn=0, xmx=1000, ymn=...
0
The house data is POINTS so readShapePOLY will fail because you are trying to read points into polygons.
This should be obvious from the error message:
> readShapePoly("./test/house.shp")
Error in .Map2PolyDF(Map, IDs = IDvar, proj4string = proj4string, force_ring = force_ring, :
NULL geometry found: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, ...
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 ...
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
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 ...
0
As a test that you haven't got the rows in the wrong order, spCbind tests the equality of the row names of the arguments. If they don't match it will stop:
> d = data.frame(ZZ=runif(56),AA=runif(56))
> ss = spCbind(scot_BNG,d) # scot_BNG from example(readOGR)
Error in spCbind(scot_BNG, d) : row names not identical
If you want to override this and you ...
0
Make two sf data frames of circles from two sets of 10 points in a unit square:
set.seed(123)
p1 = st_cast(st_sfc(st_multipoint(cbind(runif(10),runif(10)))),"POINT")
p2 = st_cast(st_sfc(st_multipoint(cbind(runif(10),runif(10)))),"POINT")
b1 =st_buffer(p1, .15)
b2 =st_buffer(p2, .15)
b1d = data.frame(id=1:10)
b1d$geometry = b1
b1d = st_as_sf(b1d)
...
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")
...
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
...
0
You are out of luck, this is tremendously hard to do! If you are serious about it you should read about differential privacy because this is probably what you are after.
When you think of this problem, you should consider the case of a recluse person living at the end of long isolated road. Do you really think you can do something about their GPS coordinate ...
Top 50 recent answers are included
