New answers tagged r
0
@Gagandeep Singh use raster::extract() to do so. ATTENTION: make sure you are using two datasets with the same crs via sp::identicalCRS()!
URL:https://www.rdocumentation.org/packages/raster/versions/2.9-23/topics/extract
raster1 <- raster()
# load specific raster-dataset (in your case global aridity index, ...)
coords <- readOGR()
# load the spatial ...
0
So I found a better solution with less code. All I had to do was use the ms_simplify function from the rmapshaper package after using the gSimplify function. The shapes load super fast now.
library(sp)
library(rgeos)
library(readr)
library(raster)
library(leaflet)
library(rmapshaper)
AUS <- getData("GADM", country = "AUS", level = 1)
CAN <- getData("...
0
As already posted in this answer, this now also works very nicely with sf:
require(sf)
fc <- sf::st_read("C:/path/to/your/filegeodatabase.gdb", layer = "some_featureclass")
But writing into a fgdb ist not implemented (yet?), you'd have to have a ArcGIS / ArcMap License as well as the R library arcgisbinding (see https://github.com/R-ArcGIS/r-bridge)
...
1
Let's make a reproducible example...
First use the raster::getData package to get a raster of the heights of Liechtenstein (because Liechtenstein is quite small so we can see what its doing, and also because I played in a band that had a number 3 hit in the country...):
> library(raster)
> library(sp)
> raster1 = getData('alt', country='LIE', mask=...
0
apologies for all these libraries, i'm not sure exactly which is needed... Good luck! this is kind of complicated. and not exactly sure how it works.
library(igraph)
library(RColorBrewer)
library(SpatialTools)
library(sp)
library(ANN)
library(rgeos)
library(maptools)
library(sp)
library(rgeos)
library(plyr)
library(FNN)
library(tidyr)
library(foreign)
...
0
Finally cracked this (I am now using Win10, ArcPro 2.2.4 and R 3.6.1, but that shouldn't matter):
library(reticulate)
arcpy3_dir <- 'C:/Program Files/ArcGIS/Pro/bin/Python'
system2(file.path(arcpy3_dir, 'Scripts', 'proenv.bat'))
Then you can
import('arcpy')
and proceed, or use {python} chunks in an R markdown document e.g.
```{python}
import arcpy
/*...
2
Make a reproducible example:
> set.seed(999)
> xy = data.frame(x=runif(24), y=runif(24))
> xy$ID = rep(c(1:6), rep(4,6))
> head(xy)
x y ID
1 0.38907138 0.8260703 1
2 0.58306072 0.8195141 1
3 0.09466569 0.5684927 1
4 0.85263123 0.6196068 1
5 0.78674676 0.8308805 2
6 0.11934226 0.4588336 2
Make the data frame into an sf ...
2
So close, use:
files = list.files("path_to_data",pattern="*.tif$", full.names=TRUE)
rs <- brick(files)
or
rs <- stack(files)
Don't use map as name, there're several R functions with the same name.
Then, to create a sum of all your rasters use calc function (from raster package):
rs1 <- calc(rs, sum)
4
The number of points above the 0.75 quantile is by definition always going be approximately 0.25. The variation you are getting is because if there's only a small number of points then the quantile computation is approximate or there's tied values. Compare:
> metrics(runif(1000))
[1] 0.25
> metrics(runif(12))
[1] 0.25
> metrics(runif(13))
[1] 0....
2
It was difficult to understand exactly what you were asking, but I think this solution using ggplot2 to view the rasters will better suit your needs.
I could only test using my own sample data, but hopefully this will translate to your data as well.
#dependencies
library(lidR)
library(ggplot2)
library(stats)
#specify file path
infile <- '...
1
Note where the error is coming from and step through your code. Its not even getting to the Savitsky-Golay line.
> v = as.vector(r)
> z = na.spline(v)
> s1.ts2 = ts(z, start=1, end=nlayers(z), frequency=12)
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘nlayers’ for signature ‘"numeric"’
The error ...
1
You can read the raster in R with the rasterpackage.
It can be loaded with map <- brick("file_path"). Then you can apply addition operation between your loaded objects. If you have 365 of them, you should read them through a loop.
2
I resolve with this script:
library(sp)
library(rgdal)
library(geosphere)
# example data from the thread
x <- c(-1.482156, -1.482318, -1.482129, -1.482880, -1.485735, -1.485770, -1.485913, -1.484275, -1.485866)
y <- c(54.90083, 54.90078, 54.90077, 54.90011, 54.89936, 54.89935, 54.89935, 54.89879, 54.89902)
# convert data to a SpatialPointsDataFrame ...
0
Make a grid by extending the bounds of your points out by half a cell size in all directions:
> cellsize = 0.5
> g2 = st_make_grid(
st_as_sfc(
st_bbox(points) +
c(-cellsize/2, -cellsize/2,
cellsize/2, cellsize/2)),
what="polygons", cellsize=cellsize)
Then your points are at the centres of ...
0
Try something like:
readandproject = function(filenames){
lapply(filenames, function(filename){
r = raster(filename)
crs(r) = "+init=epsg:4326" # fix broken CRS
r = projectRaster(r, crs="+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
return(r)
}
)
}
Then prasters = readandproject(c("file1.nc","file2.nc")) ...
0
Using these packages (MASS gets you the negative binomial distribution functions):
> library(MASS)
> library(raster)
Let's create a raster of coefficients that you have computed in some way from your landscape. Here its a 3x4 grid of numbers from 1 to 12:
> C = raster(matrix(1:12,3,4))
Then we can create negative binomial samples with a mu value ...
2
Loop over rows of the data frame and extract the i-th element from the vector:
for(i in 1:nrow(nsol)){
d = nsol$Date[i]
hillshade(..., nsol$Aspect[i], nsol$Azimuth[i],...etc...)
}
0
If you are working in R then merging is actually not required. You can use ggplot to display your data on a map. See, my example. This one liner works wonders with shape file bangladesh = fortify(p1, region = "ADM1_EN"). And, then you can use
ggplot() + geom_map(data = br2016a, aes(map_id = ADM1_EN, fill = HectareLocal),
map =...
1
After @Spacedman’s guidance I could finish my own function for what I wanted. I added the continuation to get the number of matches and mismatches and percentage of matching and made it to a function with the 1000 replicate. Here is how I did it. Probably there is a shorter and cleaner way to do it but I am a beginner so:
overlap<-function(x,y, shapefile)...
3
Using your test data, this:
> gDistance(sp_poly, sp_pts)
[1] 3.605551
returns the distance from the polygon to the nearest point, but if you add byid=TRUE:
> gDistance(sp_poly, sp_pts, byid=TRUE)
1
1 3.605551
2 5.000000
3 6.403124
4 7.810250
you get a matrix of distances - each row is one of your points. If you have more than one polygon ...
1
You would need to import the ASCII file as points, then interpolate a DEM. A quick google search will return many explanations. For example, See:
this answer to Creating DEM from Point Data using QGIS?
Then create a slope raster from the DEM. Again, google it. i.e.: Creating slope map from DEM using QGIS?
1
RQGIS only works with QGIS2, and I assume that you have already installed QGIS3. Hence, use RQGIS3, you can find it under https://github.com/jannes-m/RQGIS3. So far, it is not on CRAN due to an RStudio issue under Linux/Mac.
2
QGIS is warping from the projection of waterville into the UTM zone 16 you are using in your map - that's why it looks "curvy". See what the printout of raster("waterville.dem") provides as a "crs" (zone 17). Try looking at the dem in QGIS in its own map, you'll see the layer extent and the screen coordinates will match.
The analogous process with the ...
0
Your coords are:
[,1] [,2]
[1,] 14.650 40.350
[2,] 13.600 39.350
[3,] 0.033 -78.717
[4,] 37.017 -2.383
[5,] 46.837 -113.966
but coordinate matrices are always X then Y, ie Longitude then Latitude. Your NAs are coming from column 2 being more than 90 or less than -90, which are not valid latitudes.
extract(r, coords[,2:...
2
Lets do an example with some supplied sample data from the spData package. I know columbus isn't a grid but I guess your shape is a set of square polygons that make a grid.
library(rgdal)
shape <- readOGR(system.file("shapes/columbus.shp", package="spData")[1])
Now if your shape object doesn't have a row index, add one:
shape$row = 1:nrow(shape)
Let's ...
0
So I found a great blog post here that answers my question.
# Load the contour data and subset by layer required
readOGR("myfolder/contours.shp") %>%
subset(FEAT_TYPE =="ContourLine") -> c
# Create a target raster for the DEM
dem_bbox <- extent(c)
dem_raster <-raster(dem_bbox)
projection(dem_raster) <- CRS(projection(c))
# Set ...
0
Assuming your name is in a column called name and your spatial lines object is called d, you can do the same as with any data frame in R, and select rows by some logical operation. Here I use grepl to match a string starting NH:
> d = data.frame(name=c("M1","NH2","Not NH","NH3","NHx"),z=1:5)
> d[which(grepl("^NH",d$name)),]
name z
2 NH2 2
4 NH3 4
...
0
From documentation:
writeRaster
Write Raster Data To A File
Write an entire Raster* object to a file, using one of the many
supported formats.
When writing a file to disk, the file format is determined by the
'format=' argument if supplied, or else by the file extension (if the
extension is known). If other cases the default format is ...
0
i was able to accomplish this task by subsetting the data from the data frame which matched my classes. Here is the updates code section:
changeDet1 <- calc(stack(lc1,lc2), fun = change)
codes_ <- data.frame(ID = grid_$code,value = paste0('from ',grid_[,1],' to ',grid_[,2]))
logical_test <- which(grid_$change == T) # remove no change classes
codes_ ...
3
Finding the proper proj4-string can be a bit funny.
Such strings provide information on both the ellipsoid and how coordinates on the surface are arranged. Deciphering your example gives us the following two core elements:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
projection is UTM zone 18 with a WGS84 datum
ellipsoid is WGS84
When ...
2
Building off of genorama's answer above, you can also convert the output of bkde2D into a raster rather than contour lines, using the fhat values as the raster cell values
library("leaflet")
library("data.table")
library("sp")
library("rgdal")
# library("maptools")
library("KernSmooth")
library("raster")
inurl <- "https://data.cityofchicago.org/api/...
0
Your problem is that you seem to think that x[i] will get the i-th row of the data frame so that x[i]$year gets the year of that row. That's wrong. x[i] actually gets columns:
> x = data.frame(year=1980:1990, z=runif(11))
> x[1]
year
1 1980
2 1981
3 1982
...
> x[2]
z
1 0.7296205
2 0.3403549
...
To select rows in a 2-...
0
There is some useful information on Rotation of Coordinates Based On CORDEX Domains and you can contact to AgriMetSoft to get help.
Cheers
1
Let's work through a Haussdorf clustering of lines.
We'll use the sf package for spatial data and distance calculations:
library(sf)
starting with your final x, lets group everything by cyclone number, make line features, and keep the number of points in the group:
cyclones = x %>% group_by(CycloneNo) %>% mutate(n=n()) %>% summarize(n=mean(n),...
1
This might be a more appropriate question for stackoverflow.com, but I'll go ahead and answer. First off, you should probably combine all your different dataframes into one, with a separate column for the feature of interest:
library(dplyr)
Gcross_normal$feature <- "normal"
Gcross_tumor$feature <- "tumor"
Gcross_invasive$feature <- "invasive"
...
0
Look at the object
r1
or check if it has an associated filename:
filename(r1)
The source points to a filepath on someone's computer, and what happened was they created a raster object linked to a file i.e.
s <- stack("MOD13Q1_EVI_2008_097.tif")
and then saved it out
save(s, "myfile.RData")
All that did was save a summary of the object, like the ...
0
I imagine that your code is going sideways in the projectRaster function. You do not have the res argument defined so, it is likely that the cell resolution is not the same in the x,y dimensions. Since an ASCII Grid is a numeric array with a descriptive header, a symmetrical cell size is an absolute requirement. To address this just define the res argument ...
1
Construct a new ppp using the coordinates of the components. Create a marks parameter as a factor of the same length as the sum of both ppp sizes:
> p1 = ppp(runif(10),runif(10))
> p2 = ppp(runif(10),runif(10))
> pm = ppp(c(p1$x, p2$x),c(p1$y,p2$y),
marks=factor(c(rep("A",p1$n),rep("B",p2$n))))
> is.multitype(pm)
[1] TRUE
> plot(...
0
have you checked The markovchain Package: A Package for Easily
Handling Discrete Markov Chains in R?
the documentations has some examples that could be applied to LUCC.
have a look at the documentation here:
https://cran.r-project.org/web/packages/markovchain/vignettes/an_introduction_to_markovchain_package.pdf
2
If you want to do this you have to ceate a new class which is a subclass of RasterLayer:
> traster = setClass("traster", contains="RasterLayer", slots=c(t="numeric"))
Create a standard raster - it has no t slot:
> r =raster()
> r@t=99
Error in (function (cl, name, valueClass) :
‘t’ is not a slot in class “RasterLayer”
Convert to my traster ...
2
You are reinventing the wheel here. There is a function land.metrics in the spatialEco package that calculates point (radius) and polygon landscape metrics. The focal.lmetrics function facilitates moving window metrics. You could also take a look at the sample_lsm function in the landscapemetrics package. It is written in Rcpp so, is quite fast and efficient ...
1
As you can see, the error states that the file does not exist. The path you are passing as an argument to the raster function should not begin with \\ as the root is the D: disk.
Try removing the first two back slashes so your path looks like this:
D:\\Working\\Research work\\DDD\\Input_ParamFile\\DEM_Gil.tif
1
The reclassification matrix has to be three columns. The first two columns are the start and finish values. If you want to reclassify a single value then have it in the first and last column and do right=NA to specify open intervals on both sides.
For example this matrix will reclassify 4 as -2:
> rcl = matrix(c(1,3,-1,4,4,-2,5,20,-3),ncol=3,byrow=TRUE)
...
3
Code in that package uses gzip to decompress the .gz file:
wvar <- data.frame(data.table::fread(cmd = paste0("gzip -dc ",
files[[1]]), header = FALSE))
gzip is a standard utility with Linux but must be an optional extra for Windows that you don't have, or R isn't finding the gzip binary executable.
Its possible that R code should use the ...
0
This:
for(i in models){
rastert[[i]] <- predict(raster,i)
}
isn't going to work because you are trying to use i as a model (in predict) and as a list index (in rastert[[i]]).
To fix that, loop over the indexes (here i is 1 to the number of models) and extract the model you want inside the loop:
for(i in seq_along(models)){print(i);rastert[[i]]=...
0
I finally figured out the solution to this problem:
Using the same input from the question, and using the sp package and a workflow of Polygon >> Spatial Polygons >> Spatial Polygons dataframe (thanks R help; ?SpatialPolygons).
Prior to that I needed to find the concave hull for each cluster (using the concaveman package) and that required a matrix of ...
0
You can do something like that in FME Workbench. Here you have step-by-step instruction at knowledge.safe.com: It probably would need some modification for your purpose.
0
Now, I have rcount<-length(intersect(SpatialPoints(pts_s), r_list[[i]])) inside the first function and the problem is solved!!
SplitRas <- function(raster,ppside,save,plot){
h <- ceiling(ncol(raster)/ppside)
v <- ceiling(nrow(raster)/ppside)
agg <- aggregate(raster,fact=c(h,v))
agg[] <- 1:ncell(agg)
agg_poly ...
Top 50 recent answers are included
