New answers tagged r
0
There is no "billinear" method. It should be "bilinear".
Also, this is probably not a good approach
Settlement_rprj <-projectRaster(Settlement, crs = '+proj=longlat')
Settlement_rsmpl=resample(Settlement_rprj, NDVI, "billinear")
Instead, you could do
Settlement_rsmpl <- projectRaster(Settlement, NDVI)
1
It seems like projectRaster did not work as you intended. The extent in your example is still not in the same units.
I tried to reproduce your problem with my own data and the code you provided also gave error messages (but differnt ones) for me. So, you might want to try out project and resample from {terra}, which worked for me.
"{terra} is very ...
0
Rather than working through the headache of accessing slots using sapply, I would recommend just using rgeos::gArea (or it's sf analog sf::st_area) and filter using the vector. In R ordering is your friend. In this case, the resulting vector from rgeos::gArea will match the SpatialPolygonsDataFrame so, you can just pull indices from the vector to operate on ...
0
The error Error in match.arg(type) : 'arg' must be of length 1 is associated with trying to pass more than one value to an argument that only accepts a single value. You are passing two values to the type argument. Try it with:
tiff("f0ex_plot.tiff", height = 10, width = 20, units= "cm",
type = "windows", res = 300)
#... ...
1
I'm not sure to understand what you are looking for. To my understanding you want to retrieve the plot (i.e. the file) where the trees were found. I'd go for something like that:
opt_output_files(ctg) <- here("output/trees/{ORIGINALFILENAME}")
Then you can extract the name from the filename:
read_and_add_name <- function(path) {
name <-...
0
The above answer works, just posting this to show how I got it working in R after creating the raster table in ArcMap (thanks to radouxju). Load the raster dbf in R and attach the categorical column and save it back and load the raster again in ArcMap.
# Import the dbf
library(foreign)
dummy_dbf = read.dbf("~/Dummy_Int.tif.vat.dbf")
...
3
A tif raster with categorical values will always store the information as Integer, then a look up table is used by different softwares to translate the information from the numerical value to its labels. In orde to optimize the size, Byte should be used if less than 255 categories.
In ArcGIS, the raster will be considered as a categorical raster if it is ...
1
With regard to the naming issue, I'm not sure what you are trying to do with gsub. You don't use any indexing in "fname" so the name isn't changing between iterations of your loop, which could be why you are getting this naming system. I'd suggest something like:
fname <- paste0("filename_", i, ".tiff")
0
I think the clearest way to do this is to use raster::extract in R, applied to files that contain binary indicators ("masks") of land cover created by raster::reclassify, for each of the land types. For example:
library(raster)
#Used only to generate random raster, not needed for your application
library(spatialEco)
#For creating sf shapefile
...
2
You trained the model on a data frame, DF:
randomForest(DF$training.class ~ ., data = DF,...
and that data frame had an ID column, and becuase you did ~. it used everything (except training.class) as variables in the model.
Predictions need values of the same variables used in the training. You are doing:
predict(rstr_stack, model1)
and the raster stack ...
0
If you are using a recent version of Postgis then you can output the whole query as geojson and maybe read that using st_read? Here is a query you can try: SELECT ST_AsGeoJSON(sub.*) AS geojson FROM dataset_test AS sub. Paul Ramsey's excellent blog explains this: http://blog.cleverelephant.ca/2019/08/postgis-3-geojson.html
0
Had the same issue but got a solution by using the file's absolute path. First, get the absolute path using the tools package.
library(tools)
r1_gdb <- file_path_as_absolute("/Volumes/LaCie/2019 PRISM GIS clean/Region1_Cleaned.gdb")
Then run the ogrListLayers function.
r1_lyrlist <- ogrListLayers(r1_gdb)
That should work.
1
If you get your route into a row-oriented order with names that match the columns in the municipalities:
> route = data.frame(matrix(c("08", "019", "26", "019", "26", "039"),ncol=2, byrow=TRUE))
> route
X1 X2
1 08 019
2 26 019
3 26 039
> names(route)=c("CVE_ENT","CVE_MUN&...
0
just use "Erase tool" if you are using arcmap for example.
Or use "intersect" than select the polygons that take an ID from the "water layer"
1
The margins are controlled by the margin parameter - ?rasterVis::levelplot tells you:
margin: A list or a logical. If it is TRUE, two marginal graphics
show the column (x) and row (y) summaries of the ‘Raster*’
object. The summary is computed with the function ‘mean’.
[etc]
So by default it shows the mean of each row and ...
1
Starting with a gridded set of X,Y,class coordinates:
> d = data.frame(expand.grid(x=1:5,y=1:4))
> d$class=sample(LETTERS[1:4],nrow(d),TRUE)
if your character column isn't a factor already, make one that is:
> d$fclass=factor(d$class)
then make a numeric index column of that factor:
> d$iclass = as.numeric(d$fclass)
> head(d)
x y ...
0
The answer by Jindra Lacko is specific to NUTS regions. Just a quick fix if you have an invalid geometry after updating sf to v1.0:
sf_use_s2(FALSE)
It does not solve the problem, but it brings back the old way of working.
I still don't have a general solution to make a geometry valid when using S2, but just if it could help someone...
-1
I would try using SF with RPostgreSQL or rpostgis. My understanding of Rgdal, as described here is that is uses SP, which as mentioned here, is usually, but not always, slower than SF. Using ST_Simplify may also speed up your query.
In answer to your second question, you could query the geom column with st_astext.
2
You can consider getting the NUTS regions via {giscoR} package - it serves formally valid NUTS regions.
nuts <- giscoR::gisco_get_nuts(year = '2016')
all(sf::st_is_valid(nuts))
[1] TRUE
The example is now:
library(sf)
library(giscoR)
NUTS_WGS84 = gisco_get_nuts(year = '2016')
y_coord <- c(46.976, 46.948)
x_coord <- c(7.483, 7.45)
coordinates_as_df &...
1
You can subset the string, create a date vector, sort it and use the resulting index to order the stack.
This vector can be the result of a list.files call. However, note that the subsequent code is structured around the file names only and do not include a full path.
r <- c("SENTINEL2A_20210513-153205-267_L2A_T18NUK_C_V1-0_CLM_R1.tif",
&...
2
I'm not familiar with Leaflet in R, but your problem is most probably use of projected coordinates in setView.
Regardless of map projection, Leaflet always requires unpprojected coordinates [lat, lng] as input to its functions/methods.
This worked for me OK in JS Leaflet:
var crs28992 = new L.Proj.CRS(
'EPSG:28992',
'+proj=sterea +lat_0=52.15616055555555 ...
1
If the only issue is to resolve NAs, and you are using {dplyr} anyhow, you may consider dplyr::coalesce() for a somewhat more concise code.
library(dplyr)
shp <- tibble::tribble(~FID, ~X1, ~X5,
1, 'VEG', 'PRPU',
2, 'VEG', 'PRPU',
3, 'VEG', 'PRPU',
4, 'VEG', '...
0
I just ran into the same problem while working with the same tutorial and solved it this way .
So the argument OrgDb doesn't take in "org.Hs.eg.db" which is a string. Instead , put in the name of the annotation database directly for that argument.
Example:
gse <- gseGO(geneList=gene_list,
ont ="ALL",
keyType ...
1
You were really close. You should use == instead of %in%.
shp %>%
mutate(class = ifelse(is.na(X5), X1, X5))
# FID X1 X5 class
#1 1 VEG PRPU PRPU
#2 2 VEG PRPU PRPU
#3 3 VEG PRPU PRPU
#4 4 VEG PRPU PRPU
#5 5 WTR WTR
#6 6 WTR WTR
#7 7 WTR WTR
#8 8 VEG PLSE PLSE
#9 9 VEG PLSE PLSE
#10 ...
0
The answer provided by @Spacedman doesn't work any in newer versions of sf. I think the issue is that aggregate automatically handles the union of the geometries and the FUN argument is intended to handle aggregating the non-geometry columns? I also think st_combine produces an sfc output, so the geometry column ends up as a list of sfc objects, not a single ...
1
There's some things in the function I don't quite understand, and I'd probably do this a different way, but here's my thoughts:
This line makes a stack by stacking cost on top of itself. This sets up a stack on which we are going to further stack the cost rasters for each of the points we have been given. I don't know why two copies of cost are stacked here, ...
2
If you ask a model to predict over a raster, it will predict over the entire raster. It doesn't know there's a river there.
Its prediction for location-types that you haven't trained on will be inaccurate. Maybe the nearest class for the river is "desert", and it will think rivers are deserts.
Either add some points on the river (or other water ...
0
I don't think that the CMIP6 data for 30s is available anywhere else at this time. If you are looking for modelled data with such spatial resolution, the only data available is SMIP5 via https://www.worldclim.org/data/v1.4/cmip5_30s.html
But I am also hoping the developers will deposit the data soon. Because it is already over a year past when the data is ...
2
In the specific context you describe - having top right and bottom left corner of your desired polygon - you may be able to get by with sf::st_bbox(). It returns the bounding box of an object, in your case of the Poly_Coords_df data frame (as intepreted in context of EPSG:32611).
poly <- Poly_Coord_df %>%
st_as_sf(coords = c("lon", "...
2
You can first resample one raster with respect to the other to match the resolution, and then stack the two.
library(raster)
# Load rasters
et = raster("path/et.tiff)
sum_raster = raster("path/sum_raster.tiff)
# Get required extent
extent(sum_raster)
# Store extent values (coordinates) in a variable
Extent = c(-160, 10, 30, 60)
# Crop
et_crop ...
1
You've only got two points so you need to construct five to complete a box.
I don't understand what all that pipeline is supposed to do, the simplest and most direct way is to use st_polygon and feed it a list of a matrix which has five rows and two columns by taking the elements from your data frame:
> pol = st_polygon(
list(
cbind(
...
0
This is the answer:
library(raster)
stacked <- stack("path/temp15.tif")
meanR <- calc(stacked, fun = mean)
sdR <- calc(stacked, fun = sd)
upper <- meanR + sdR * 2
lower <- meanR - sdR * 2
result <- stack()
for (i in nlayers(stacked)){
layers = stacked
layers[layers > upper] = NA
layers[layers < lower] = NA
result ...
0
The numbers in overlayKR extent are very large...
> overlayKR@extent
class : Extent
xmin : -293508
xmax : 822492
ymin : 1914749
ymax : 2929749
but the CRS says these are lat-long degrees:
> overlayKR@crs
CRS arguments: +proj=longlat +datum=WGS84 +no_defs
which is clearly not true.
overlayBRD looks plausible, because ...
2
As an alternative to the approach taken by elio-diaz I propose the following; the logic is built on on
identifying the "lakes on lakes" as polygons touching the multipolygon (the Fig.2 from original post) and
erasing the island area from the polygon layer via rmapshaper::ms_erase()
It is interesting that for some unknown (to me) reason sf::...
2
After taking a look at the data, and unless an osm expert points out a better filtering option, I suggest to filter out those polygons that "have" the islands and the islands as well, although they are not stored as polygons with rings (holes), in this case the big lake and the islands are three separate polygons. We may take them out using ...
0
See the answer to this question (as Erik pointed out) for how to do this with raster. Here I show something similar with terra
When asking an R question, always include some example data:
library(terra)
s <- rast(system.file("ex/logo.tif", package="terra"))
s
#class : SpatRaster
#dimensions : 77, 101, 3 (nrow, ncol, nlyr)
#...
1
Apparently, your loop part is "growing" objects. (https://www.burns-stat.com/pages/Tutor/R_inferno.pdf, page 12). When you use rbind(df_train, df) within a loop it has to duplicate the entire object on each iteration. You may speed things up by pre-allocating a vector, e.g. my_list <- vector('list', n) where n is the size of length(unique(...
1
Instead of looping through 0:20 you may use table function which gives you the count per each different value.
library(exactextractr)
library(raster)
# build a raster layer
rast <- raster::raster(matrix(sample(20, 100, replace = T), ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
# and a polygon layer
poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
...
0
If I understand correctly, there are point and polygon shapefiles in your path (shp can't have two geometry types). You may query with st_geometry_type after having all layers bound together with bind_rows.
library(sf)
library(dplyr)
nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)
nc_centr = st_centroid(nc) # just ...
1
I have solved it in case someone needed;
recode_poly <- function(x){
poly_recode = c()
for (i in (1:1000)){
if (x[i] == 1) {
poly_label = "cloud"
} else {
poly_label = "clear"
}
poly_recode <- c(poly_recode, poly_label)
}
return (poly_recode)
}
then you can convert poly_recode to a dataframe
but ...
1
That is odd; it suggests that something is missed in assessing the memory needs. You can check for yourself what is going on under the hood with
canProcessInMemory(NDVI_Stack, verbose=TRUE)
And you can allow raster to use less memory, for example like this
rasterOptions(memfrac=.3)
(and/or using smaller chunk-sizes), see ?rasterOptions. And see the ...
3
crop can only cut entire cells (rows, columns) so that won't help you here.
If you want to ignore the small difference, and the rasters have the same number of rows and columns, you can do
library(raster)
r1 = raster::brick("path/rast1.tif")
r2 = raster::brick("path/rast2.tif")
extent(r2) <- extent(r1)
s <- stack(r1, r2)
That is ...
0
Final solution:
library("sp")
library("sf")
library("raster")
library("parallel")
library("ggplot2")
df=read.table("~/centroids.csv",sep=",",header=TRUE)
df2=df
coordinates(df2)=~longitude+latitude
proj4string(df2)=CRS("+init=epsg:4326") # points are initially projected in ...
1
In your current code, what you are doing is converting the raster into a CRS object, that's why it fails in the next line. To correctly set the CRS you have to do the following.
country_template <- raster()
dimensions <- extent(124000, 463000, 2953600, 3127600)
country_template <- setExtent(country_template, dimensions)
res(country_template) <-...
1
In stead of
RGB = raster("path/RGB.tif")
You should do
RGB = brick("path/RGB.tif")
Or even better:
library(terra)
r = rast("path/RGB.tif")
plotRGB(r)
# or
RGB(r) <- 1:3
plot(r)
1
Here is an improved version of your code. As you do not provide example data I cannot say if it works better.
library(raster)
projection <- "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"
studyarea <- spTransform(studyarea, CRSobj=projection)
# make a template RasterLayer
template <- raster(studyarea, res=1000)
# use the template
...
0
Did you inspect filelist? It probably includes .zip or other files. To only select files that end on .tif do
filelist <- list.files(pattern = "\\.tif$")
$ means "the filename must end on these preceding characters", and \\. means "there really must be a point". If you did .tif$ the point symbol would be interpreted as a ...
0
I might not have the correct language to exactly explain this, but I just discovered a case where using sf::st_combine before calling sf::st_union made a big difference to me.
I was using dplyr::group_by and dplyr::summarise to group an sf tibble by a certain variable and then union each group to a single sf object:
sf_tibble %>%
dplyr::group_by(area) %&...
0
Check the following code:
# Load packages
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(mapview)
library(spdep)
#> Loading required package: sp
#> Loading required package: spData
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://...
0
I figured that there seems to be no alternative to osmdata which is why I decided to craft my own solution.
First, I formatted the query similar to the code I posted in my question:
query.osm.boundaries <- function(admin_level, region) {
osm_id <- tmaptools::geocode_OSM(region, details = TRUE)$osm_id
query <- paste0(
'[out:json][timeout:100];...
Top 50 recent answers are included
