New answers tagged r
0
The csv export from that site gives the extent of the bounding box in the order (xmin,ymin,xmax,ymax) so you can create a data frame with the vertices of that box, turn it into a spatial object, and then transform to your desired CRS:
library(rgdal)
bbox = c(5.9559,45.818,10.4921,47.8084)
vertices = as.data.frame(rbind(bbox[c(1,2)],
bbox[c(3,2)]...
1
Try this:
recalc<-function(ndvi,pv){
ifelse(ndvi<0.2,0.92,ifelse(ndvi >= 0.2 & ndvi <= 0.5,0.004 * pv + 0.986,ifelse(ndvi >0.5,0.99,NA)))
}
r<-stack(ndvi,pv)
res <- overlay(r, fun = recalc)
As per updated comment the solution that worked was instead using calc:
res<-calc(r,fun=recalc)
The recalc function is basically ...
1
There's a couple of things going on here, I think. First, fasterize doesn't have ... enabled so you can't pass na.rm = TRUE through your summarising function, and all of those default to FALSE.
The other issue is that as the documentation warns, the output of st_combine() is often not valid under the simple features spec so may perform unexpectedly. If you ...
1
I was prompted to have a look closer at the GPS time stamps which the code uses to assign pulse IDs and then thin when use_pulse=TRUE. Out of my 5,000+ tiles, it turns out a few tiles had the GPS time stamps removed from the whole tile and thus, wasn't thinning because it only assigned one pulse for the whole tile (4+million points!).
I ended up writing a ...
0
I finally used rf<-writeRaster(px, filename="predict.tif", datatype="INT1U", format="GTiff", overwrite=TRUE) to export a integer type tiff file. And then I imported it in ArcGIS, shut down the background processing, and then I can create an attribute table for this raster layer using Arctool :)
0
You can just use raster::extract to extract values for a polygon (buffer), or an area around a point, and then pass the values to a table function to tabulate.
Create some data
library(raster)
library(sp)
library(rgeos)
r <- raster::raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540,
ymx=4453690, resolution=270, crs =...
0
Simplest way is to extract coordinates for the set of cells over the threshold:
> r = raster(matrix(runif(50*50),50,50))
> coords = coordinates(r)[r[]>0.9,]
This should show the raster with some dots on the greenest points:
> plot(r)
> points(coords)
2
so I played around a little creating an sp objet
rs <- raster(file_path, band=2)
sp <- SpatialPoints(rs)
res <- extract(rs, sp,
fun=max,
df=T,
cellnumbers=F)
then kept those with high intensity
res <- res[res$All_Cells > 0.9,]
coord <- xyFromCell(rs, cell=res$ID)
plot(coord, cex=0.05, col="...
0
Here's one approach that will give you the coordinates of each pixel that meets a criteria:
# Generate a raster to play with
library(raster)
library(RStoolbox)
mat = matrix(c(1,2,4,4),nrow=2,ncol=2)
ras=raster(mat)
# Use ggplot2::fortify() to turn it into a data.frame
ras.fort = ggplot2::fortify(ras)
# Select the raster values you're interested in (here, ...
3
If these values:
bbox: xmin: 313368.6 ymin: 4832795 xmax: 314211.3 ymax: 4833629
(which are the range of the data location values) are latitude-longitude:
proj4string: +proj=longlat +ellps=clrk66 +no_defs
then something is terribly wrong with the shape of the earth.
Lat-long coordinates should only ever be (-180,180),(-90,90) - your ...
1
To calculate the product (or any in-formula calculation) use the I() function:
mod <- lm(A ~ B + C + I(B*C))
B and C should be included by themselves (as noted above). If they are not significant and the interaction term is, still leave them in.
1
Firstly, you need them to be the same CRS and resolution in order to stack them. For this. Assuming your 5k raster also has a bigger extent for the purposes of this example, you can use raster::projectRaster():
bigger_raster <- projectRaster(raster5k, smaller_raster)
Now, if your extents are still off - you'll need to use raster::crop() to fix this.
...
3
The lidR package relies on the rlas package to read and write las file. The rlas package has a recent support of LAS 1.4 files (v1.3.0 release date: 2019-02-03). Moreover the point record formats >6 are a bit different than former point formats. Your code is correct and you actually found a bug in function write.las from rlas that occurs with point format 6 (...
0
In order to answer this question, I first had to understand the output of sf::st_coordinates which returns a matrix with the following rows: X,Y,L1, L2, L3.
X: X-coordinate
X: Y-coordinate
L1: differentiates between convex hull and hole coordinates for each polygon
L2: differentiates between polygons in multipolygon
L3: differentiates between ...
2
There isn't an assignment for x in the code above so I can't generate your exact output, but you should be able to plot non-0 entries in your histogram by using hist(r[r!=0]).
Here's an assignment for x (just a placeholder) and some artificially-crafted 0's for the raster values that demonstrate the difference between hist(r) and hist(r[r!=0]):
x = seq(-1,...
3
In addition to @Spacedman's answer, you can set up a loop in R to compress tiffs using LZW with the writeRaster function in the raster package, which still uses GDAL. The options argument allows you to apply LZW compression.
In this example, the file is not overwritten but rather has an "_LZW" appended to the original name. To just overwrite the original ...
4
You can do this using the command line tool gdal_translate. This is available for Windows, Linux, and Mac OSs (you don't state your OS).
Running:
gdalinfo none.tif
will show the info on the file, including the compression type and the locations:
Driver: GTiff/GeoTIFF
Files: none.tif
Size is 204, 228
Coordinate System is:
GEOGCS["WGS 84",
DATUM["...
0
It is also possible you don't have even read permissions over the folder you are trying to open a layer from. This problem is especially pesky if running stuff from a docker container. In order to give yourself permission just type:
chmod -R 777 /path/to/folder
in terminal.
0
It sounds like 'Extract by Mask' is what you need. Here's info... https://support.esri.com/en/technical-article/000013000
0
Use the extract() function to find all values associated with a particular location, and visualize with hist(). Your raster's projection and extent will determine what points you want to use for PointOfInterest. A reproducible example:
library(raster)
r1 = raster(matrix(c(1:25),
nrow=5,ncol=5))
r2 = raster(matrix(rep(5, 25),
...
1
The addLayer function works by stacking all its arguments - to add a layer to an existing raster or stack make it the first argument. This is shown in the help:
file <- system.file("external/test.grd", package="raster")
s <- stack(file, file, file)
r <- raster(file)
s <- addLayer(s, r/2, r*2)
Running with only one argument produces that ...
0
To find the countries of a data.frame with lat-long coordinates, convert them to an sf object with sf::st_as_sf(mydf, coords=c("lon","lat"), crs=4326).
Then intersect it with a map of polygons. This is vectorized, so no loops are needed.
Here's a complete working example:
# Dataframe with latlong coordinates:
d <- read.table(sep=",", header=TRUE, text=
"...
1
Your terminology is bit unclear but I think it boils down to you trying to save spatial features vectors and expecting each element to have some attributes. You can't do that with vectors, you have to make spatial features data frames.
If you have an object that is an sfc class:
> class(poly3)
[1] "sfc_GEOMETRY" "sfc"
then saving it as a ...
1
As Spacedman sujested, you should share with us more than that, almost a piece of your tibble data. In any case, and without knowledge about your data sources, projection, etc, you could find those "cells" where your points are just by searching the closest cell center. This is just valid (and logic) if you tibble represents a regular raster grid. As I said, ...
2
First I'll set up some data like yours:
Extents and coordinate systems:
> eland = extent(-79.39998, -79.32998, 43.61302, 43.63498)
> cland = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
> ewater=extent(-79.39958, -79.32958, 43.61292, 43.63458)
> cwater = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
Make land ...
0
It is possible if done via LandViewer. Once you locate an area of your interest you will have an option to filter the images based on cloudiness level on the top right and download them with ease via download button on the right. An example of such imagery can be seen here or you can view It on the platform directly here.
0
Have you tried using alternative software to download those images? For example, MODIS data can be accessed via LandViewer. All you need to do is locate an area of your interest and select MODIS imagery in the filter section. Once done you will have the option to download the data on the right.
An example of such imagery can be seen HERE or view this ...
-1
I think this video is helpful for getting data from MOD16A2 HDF File
0
I found a really simple and good solution with the rasterize function from raster package
library(raster)
library(sp)
Create polygons
x_coord <- c(0, 1, 6, 4, 2)
y_coord <- c(10, 5, 5, 12, 10)
xym <- cbind(x_coord, y_coord)
p1 = Polygon(xym)
p2 = Polygon(xym - 3)
ps1 = Polygons(list(p1), 1)
ps2 = Polygons(list(p2), 2)
sps = SpatialPolygons(...
2
I managed to convert the coordinates using the rgdal package:
temp <- data.frame (x = c(598223, 598812, 598824, 598232, 597614, 597629),
y = c(7095460, 7095426, 7094827, 7094227, 7094821, 7095433))
temp <- SpatialPoints (temp, proj4string = CRS ('+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'))
cord.WGS84 <- ...
0
in your particular context, your point divides the original line in the longitude axis. Thus, you could create two other lines, each of which would have all points of the original one to the left and right of the longitude of the point of interest. A clumsy way of doing so is to filter the coordinates of the SpatialLinesDataFrame out of a data frame
library(...
1
You can read in irregular data if you know the max number of columns by giving column names and a fill argument:
> data = read.table("./lines.txt",col.names=c("x","y"),fill=TRUE, stringsAsFactors=FALSE)
> data
x y
1 ID1 NA
2 3285.48 -63.32
3 3285.14 -64.14
4 3284.67 -63.56
5 3285.00 -62.77
6 ID2 NA
7 3299.84 -76.82
...
0
In the end, I figured out a method with the kknn package. If somebody is still interested in this answer, I recommend reading my blog post on the topic.
0
R and raster crops to a bounding box so you will need to also mask to only retain the area of interest. Mask works by setting values inside or outside of a shape to a value usually NA in this context.
See raster::mask command
2
Your data must be projected to WGS84 before plotting:
adm_proj <- spTransform(adm, ‘+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs’)
leaflet() %>%
addTiles() %>%
addPolygons(data=adm_proj, weight = 2, fillColor = "yellow", popup=popup)
0
a good and fast way to start is to create your fishnet using the raster package to create your grid as a raster and then transform it into a polygon.
library(raster)
fishnet <- raster() # this creates a raster with the extent of the whole world.
res(fishnet) <- .5 # sets the resolution of your grid -- here .5 x .5 degrees
crs(fishnet) <- CRS("+...
0
Now we know enough to reproduce your error message, I'll do that with a data frame of lat-long numbers:
> library(sp)
> s_gis = data.frame(long=runif(5), lat=runif(5))
that should be equivalent to your s_gis with lat-long selected. Continue with your code to make sps:
> p = Polygon(s_gis)
> ps = Polygons(list(p), 1000)
> sps = ...
2
Thanks a lot, @Spacedman and @JepsonNomad for your comments. The comment by @Spacedman was very helpful. This line of code worked for me:
soil_raster2 <- rasterFromXYZ(moisture_data,res = 0.125,
crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",
digits = 3)
plot(...
0
Use SpatialPointsDataFrame with data.frame to coerce into the desired object.
library(sp)
x <- SpatialPoints( rbind(c(1.5, 2), c(2.5, 2), c(0.5, 0.5), c(1, 0.25),
c(1.5, 0), c(2, 0), c(2.5, 0), c(3, 0.25),
c(3.5, 0.5)))
class(x)
Here we create a data.frame within the call to SpatialPointsDataFrame ...
2
sf::st_intersection() will work with the last version of sf (0.7-3). I don't know why it was not working in the previous version.
0
In v2 of Q you can add it in as a File/System vector layer:
I would suspect you can do this in v3 too.
2
Given x a SpatialPoints object:
> x
class : SpatialPoints
features : 50
extent : 0.0006317429, 0.9926516, 0.02675848, 0.9901886 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
...you can convert to SpatialPointsDataFrame with as:
> as(x,"SpatialPointsDataFrame")
...
0
There is now a sf::st_crop function available but, you may need the development github version of sf. You can download the development version using the devtools/remotes package(s).
install.packages("devtools")
remotes::install_github("r-spatial/sf")
For specific irregular geometries, the raster::intersect function can clip data, and retain attributes, ...
0
Luckily, I was able to figure this out. It is necessary to set the GDAL_DATA directory in the same terminal session, open R, then install rgdal.
Set the GDAL_DATA environmental variable:
GDAL_DATA=/usr/share/gdal/
Location will vary by system; another good guess is /usr/local/share/gdal/. You should use find to locate your directory.
Then, in the same ...
1
The following code works including the load_output option. Just out of curiosity, are you running RQGIS3 from within RStudio? Because in my case running RQGIS3::open_app() crashes the RStudio R session. This was also confirmed by several other users (see the github issue tracker of RQGIS3). Since RQGIS3 works when run from the CLI, I am not really sure what ...
0
geom_polygon() is what you are looking after. The line
ggplot(Coordinates) + geom_polygon(aes(long, lat, group = grid))
should do the trick.
ps: next time, make your data available so that we can be more precise on addressing you issue :)
1
There's no mention of what happens with layer names in the documentation for calc so I suspect the answer is "just because". If you rely on layer names in your code then you should probably explicitly set them any time you think they might change.
Note that arithmetic can change layer names - even though both operands here have the same name, the output is ...
1
For me, run_qgis creates a bundle of .sdat files in the working directory:
> list.files(wd,pattern="*.sdat$")
[1] "aspect.sdat" "crossSectionalCurvature.sdat"
[3] "flowLineCurvature.sdat" "generalCurvature.sdat"
[5] "longitudinalCurvature.sdat" "maximalCurvature.sdat"
[7] "minimalCurvature.sdat" "...
0
Only half an answer: maxent() returns an object of the class MaxEnt. The fact that this is an S4-type object doesn't matter much here, but if you want top know more about that you can read on here. The actual "problem" is that the MaxEnt class does not provide a method for the generic residuals() function (resid() is just a shortcut for residuals()). I other ...
0
Convert the bnd object to an sp object using bnd2sp from the R2BayesX package.
Join data from the LeukSurv data frame using the region variable as the row number in the sp object.
You don't need to do anything with coordinate systems. The data are scaled to an approximate unit square and there's no need to georeference them for an analysis.
There's more ...
Top 50 recent answers are included
