New answers tagged r
1
R thinks one or more of your polygons are holes, and a hole has to be inside a non-hole.
Hole detection is done in the Polygon function:
Polygon(coords, hole=as.logical(NA))
coords: 2-column numeric matrix with coordinates; first point (row)
should equal last coordinates (row); if the hole argument is
not given, the status of the ...
0
You are very close, although mixing code from spand sflibraries; you only attached the shp file, so I couldn't download your file, but made a dummy linestring object from your points so you can see how to plot a line on a raster. It is easier to use lower functions (par(mfrow = c(2,1)) to plot the line on both rasters of the stack
library(raster)
library(sf)
...
0
The filtering cannot be done in your call to list.files, you will have to read the files and check one by one and filter out those which do not meet the condition:
library(raster)
library(dplyr)
raster_list = list.files(pattern = "modis", full.names = T)
raster_stack = stack()
for(i in raster_list ) {
temp = raster(i)
if(dplyr::between(temp[][...
4
Function list.files() from base R allows to list files that match a pattern, in your case it would be .tif, or, if you only want those that start with wc it would be wc.*tif, using regular expressions. The argument full.names = T will give you complete names with directories, which helps to call these file names in functions; recursive = T tells the function ...
3
Aside from how the rule of large numbers come into play in that you no longer are working with a sample but rather a population. Conceptually, here is another good example on why what you are trying is not a stable approach.
I am simulating a very skewed distribution. As you will see, the mean and standard deviations are rather meaningless (as central ...
2
You can use na.rm = TRUE in normalize_height(). The problematic points are likely to be in the buffer anyway. Also in version 3.0.3 tin() enables to control on how interpolation is extrapolated outside the ground points. Check the doc. What changed recently is that it no longer extrapolates for points that are too far from any ground point (50 m) because it ...
4
You are referring to a simple type of point pattern analysis (PPA) called Quadrat analysis. Quadrat analysis can be used both for assessing goodness of fit (of empirical data to a theoretical model) and for assessing the independence of two distributions. A chi-squared statistic is commonly used in evaluating the hypothesis of a random point process, using a ...
0
Here is a function that ask for three parameters:
x --> sf_layer
n --> number of rings
d --> distance between rings (m)
Moreover this add and ID to each ring based on the distance transformed to km.
#-------------------------------------------------------------------------------
library(sf); library(mapview)
# Create a ringed buffer with a number ...
1
The Intensity int column has the ability to store all the treeID values. There is a slightly different format for copying rows to another and you can find the locations using str(las1). Additionally, some NA values were discovered from points not allocated, which I have remedied below as this initially failed to enable writing the .las successfully. Once ...
3
There is no select, mutate or any other dplyr verbs for LAS objects. LAS objects are not data.frame but... LAS objects.
A las file contains a set of core attributes. To save extra attributes such as treeID you must add extra bytes attributes. segment_trees() does it automatically. You can verify that by looking into the header:
library(lidR)
LASfile <- ...
2
The usual trick to turning a loop over i into an apply style loop is to apply over i and write a function that does the job for argument i. For example:
sapply(1:nrow(alaska_islands),
function(i){
min(
st_distance(
alaska_islands[i,],
alaska_all[alaska_all$area > alaska_islands$area[i],])
)
}
)
gives:
[1] 59657....
0
If you want to do all your computing on that region, you may use the function raster::crop(your_raster_object, your_vector_object); raster::extent(raster_or_vector_object) returns the bounding box of a geospatial object, which allows to add easily as if it were a buffer, e.g. extent(a_vector_in_utm) + 100 will add 100 m to the bounding box.
There's a manual ...
3
The prepared flag is essentially a hint to the underlying GEOS library that's doing the actual work behind the scenes.
If x is "prepared" it will be internally indexed so that repeated calls to the predicate (contains, intersects, etc.) will be faster. For example, if you want to test whether many points are within a single polygon, it's faster for ...
-2
It's a binary tag telling you about the output of a geoprocess without actually doing the process. X ought be the first entry in the commannd (ie 'g' in 'st_contains(g, pt3.sfg, prepared = TRUE').
Language implies that st_contains is meant to be used in a loop, so that either the first or second items in 'st_contains(g, pt3.sfg, prepared = TRUE)' could be a ...
2
layer 1 is Guadeloupe. In the metadata XML:
<Layer queryable="1">
<Name>1</Name>
<Title>U2018_CLC2018_V2020_20u1_FR_GLP_R_WM</Title>
which then goes:
<westBoundLongitude>-61.906490</westBoundLongitude>
<eastBoundLongitude>-60.904868</eastBoundLongitude>
<southBoundLatitude>15.735481</...
0
I think I solved it. The URL is correct but the layer name should have been "12" (not "1" as in my answer above). For those that may struggle to interpret these WMS XML files in future, I would recommend using QGIS (or similar) to get an easy-to-read list of WMS files e.g. http://www.qgistutorials.com/en/docs/working_with_wms.html
1
read.zoo does create the object. However some zoo methods do not work properly when there are more than one value per date. You can remove duplicate dates by averaging or by using the last observation or whatever that fits you.
# remove duplicated indexes by averaging
ndvi_bfast_mean<-aggregate(ndvi_bfast$mean, list(ndvi_bfast$date), mean)
# remove ...
2
clipped[i] returns a list with one element. You meant clipped[[i]] with double brackets. Also clipped[i] within paste0() is meaningless.
Actually your use case is already covered by the package. Assuming that your shapefile has an attribute PLOTID (replace with something else if needed)
ctg <- readLAScatalog("folder/plot_a3.las", filter = "-...
2
The fix for these images was a two part solution, one for the stripes and a further one for the CRS:
First: hiding .aux file from the search path, which had redundant information to the header (.hdr), as mentioned by @Michael Stimson, this can be done by renaming the file extension or by removing the file. This fixes the stripes on the picture.
Second: QGIS ...
1
I had this same problem. I created an R package that calculates the distance matrix on the fly, so it takes much less memory to calculate Moran's I. Its also quite fast. You can find it at github.com/mcooper/moranfast.
1
Here is how you return the four corner coordinates of the raster extent. If you use sp::coordinates on the extent object then the polygon centroid coordinate is returned.
Here is an example of your raster, note that I do not have to assign a projection because everything will remain in the same projection space without it defined.
library(raster)
r <- ...
1
The extent of that raster object is expressed in units according to the crs. So, to get "coordinates" (presumably you're looking for longitude and latitude) you need to create a spatial object with that extent information, and convert it to the desired crs. Based entirely off this answer, you can convert your extent of interest using the following:
...
5
One of the main difference you see between the two is that R is willing to go for non-square pixels, while Arc enforces exactly square pixels by default.
In your input dataset, your pixels appear to be square, but in fact, in a projected coordinate system, the geographical coordinates are not a square (except at equator). As such, your pixels should not be ...
2
Raster files, like tiff images, cannot be directly converted into vector formats. However, it is possible to convert each pixel into small rectangular polygon and save the polygons as vectors into some vector format like shapefile. Some utility programs, like gdal_polygonize https://gdal.org/programs/gdal_polygonize.html, can do a little bit more and merge ...
2
To use the first raster (let me call it Raster1) as the base layer to accommodate the values from the second one (Raster2), you can use Create Constant Raster Layer tool in the QGIS Processing Toolbox > Raster Creation.
(1) Click on the small ellipsis ... button and set the Raster1 as the basis for the extent calculation.
(2) Pixel size should be the ...
1
So you have a data frame with two columns. You should show part of that so we know what you have. I've created a random one:
df <- data.frame(
ndvi = runif(4321,0,1),
ts = runif(4321, 23, 57)
)
head(df)
ndvi ts
1 0.66967785 36.47674
2 0.89686408 27.06128
3 0.11688687 37.42891
4 0.07607912 43.13066
5 0.16757402 35.98790
6 0.72767844 ...
0
You can extract data from GRIB files using the standalone program CDO.
To extract from a single point, I've previously used this command in DOS / CMD:
cdo.exe -output -select,param=%PARAMETERS% -remapnn,lon=%LON%/lat=%LAT% inputGribFile.grb > outputfile.txt
Here, you need to specify which parameters from the grb-file that you need, which you can also use ...
0
the tutorial has been updated to take into account this error.
The call to gdal_translate was failing when providing the global bounding box. The call is working without those parameters.
1
To anyone who is struggling with MOD07_L2 product and was wondering how I solved it, I followed this answer from @fdetsch to Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal, but modified a little bit. So basically what I did is the following.
library(raster)
library(rgdal)
library(gdalUtils)
# hdf files directory
main_dir <- file.path(...
3
You can write a simple function using prop.table and table to return proportions of multiple classes. The catch is that you have to know what all of the classes are before hand so you can fix the number of expected classes.
Here is an example of what is going on.
Here we set our "known" classes and then set up a loop that randomly samples a vector ...
3
You can get the frequency of a single class by passing a custom summary function to exact_extract. For example, to get the fraction of pixels that have a value of 1, you could run:
exact_extract(rast, polys, function(value, fraction) {
sum(fraction[value == 1]) / sum(fraction)
})
If you have an arbitrary number of classes, here's a solution that will ...
1
This is a needed addition to this post and a good example of how "things are always changing". With the pending release of sf 1.0 the interface with the s2 geometry library and changes in GDAL and PROJ handle projections, the flat earth model will no longer be supported for geographic coordinate operations. A spherical geometry will now be used for ...
0
I managed to find a solution myself. When exporting the raster, the clipping extent has to be maintained (the box must be checked), so that the raster extent of the exported file does not jump back to the extent that logically follows from the resolution.
0
Very simple. Don't construct the URL yourself.
leaflet() %>% addProviderTiles("MapBox",
options = providerTileOptions(id = "bweinstein.CypressCity_03_25_2020",
accessToken ="########"))
1
You could try re-running the viewshed analysis and setting the "snap raster" environment variable to one of the other resultant grids, i.e. your R1-4 and see if the outputs match then.
1
It is possible to use basemaps in the non-interactive mode of the tmap package.
However, it requires two steps:
Downloading/preparing raster data with three layers (RGB)
Start making the map with tm_shape(my_rgb_data) + tm_rgb()
You can see some examples of using this approach at https://stackoverflow.com/a/56972352/2602477.
7
If you use st_read with a query that selects zero rows you get back an empty dataframe with all the columns intact:
> nc_sql = st_read(system.file("shape/nc.shp", package="sf"), query="select * from \"nc\" limit 0")
> dim(nc_sql)
[1] 0 15
> names(nc_sql)
[1] "AREA" "PERIMETER" ...
2
I have no experience with R, however the mentioned IDs are nodes belonging to ways which in turn are tagged as fast food:
node 2439442254 is part of way 615644218
node 2439828308 is part of way 235901529
node 2439828379 is also part of way 235901529
So I presume you are just looking at the wrong information in R.
Note that OSM has three basic elements: ...
1
Hi I hope this help you, follow steps, also you can reproducible it with the next data in this link
library(raster)
library(tmap)
library(cptcity)
library(tidyverse)
# reading a raster ---------------------------------------
dem <- raster('../Data/DEM.gpkg')
# definition some class for raster -----------------------
class <- c( -Inf, 3000 , 1,
...
0
Just watch out for TMS-Tiles. tiler does export TMS-Tiles but even the viewer seems to be unprepared for that.
I used leaflet-react and set the tms property with value "true" and it worked:
<TileLayer
attribution="© contributors"
url={
"/static/layers/" +
currentProject.projectname +
&...
2
I ran into the same issue a few days ago. What seemed to work for me was adding a projection before doing the conversion. I also noticed that using a SpatialGridDataFrame (instead of a RasterLayer) object worked much faster. So I suggest you try something like this (assuming you use the WGS84 CRS)
library(sf)
library(rgdal)
library(stars)
library(dplyr)
...
1
The Emerging Hot Spot Analysis tool does not use a Kendal statistic but, rather a space-time variant of the Getis-Ord Gi* to identify space-time clusters. The Kendal Tau is for analysis of temporal correlation and monotonic trend. If you want the Kendal Tau I have both a time-series vector and raster implementation in the spatialEco R package but, remember ...
0
The issue has been resolved. The data was being based to leaflet leaflet(data = hfxETA) which was overwriting the filter. Removing the data from the leaflet call fixed the issue.
0
If this is, in fact, your exact code, you are not assigning spTransform(train_ICC, CRS(ICC_m2_clip)) to a object (existing or new). It looks like you are running the function and then assigning the projection string back to the original non-transformed object.
I imagine that simply modifying your code to overwrite the original vector sp object will fix your ...
4
Without having tried your code, what I see is:
Don't do
route[i ,] <- ...
because in general modifying a data frame in place can be slow.
Do create a list upfront to hold each geometry, i.e.
vector("list", rsize)
and put each st_linestring() in that, wrap the entire list after the loop in st_sfc(), and do the st_sf() at the end with the sfc ...
0
You've got a minus sign before the Latitude:
cord.dec = SpatialPoints(cbind(mycsv$Longitude,
- mycsv$Latitude),
Which is giving you a negative Y coordinate:
"I get: 698714, -4308793 which is obviously not correct."
Fix that and you get:
> cord.dec = SpatialPoints(cbind(mycsv$Longitude,mycsv$Latitude),proj4string=...
0
I get
x <- cbind(-120.71, 38.906)
library(sp)
cord.dec = SpatialPoints(x, proj4string=CRS("+proj=longlat"))
coordinates(spTransform(cord.dec, CRS("+init=epsg:32610")))
#> coords.x1 coords.x2
#> [1,] 698568.5 4308838
## sf is the new and entirely different 'sp' plus life raft
sf::sf_project(from = "EPSG:4326", ...
0
You are looking for an spatial join, after converting your lat-longs to a spatial object. This is done easily with package sf, replace the filenames with yours on the following code:
library(sf)
library(dplyr)
lat_longs = read.csv(latlongs.csv) %>%
st_as_sf(coords = c("longitude_column", "latitude_column")) %>%
...
Top 50 recent answers are included
