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["...
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 ...
3
Your reclassification matrix isn't right - I suspect the 0.990 on rows 3 and 4 should be 0.099 - as it is the rows make up overlapping classes so the function throws an error:
> reclmat
[,1] [,2] [,3]
[1,] -0.500 -0.251 1
[2,] -0.251 -0.101 2
[3,] -0.101 0.990 3
[4,] 0.990 0.269 4
[5,] 0.269 0.439 5
[6,] 0.439 0.659 6
[7,...
2
Starting with nc borders from the example(st_read) code, I'll first make a set of lines like in your data and then put them back together again. So we start with 100 North Carolina polygons in a spatial sf data frame:
> dim(nc)
[1] 100 15
First, to take them apart, convert to lines and intersect them so the overlaps are gone:
>
> ncl = st_cast(...
2
You could also do this using dplyr's group_by() and summarize() functions:
states %>%
group_by(group_var) %>%
summarize(geometry = st_union(geometry))
2
You were close; just cast the second argument in aggregate() as a list, like this:
x <- aggregate(states,
by = list(states$group_var),
FUN = mean)
x
plot(x)
2
For shapefile it has to be the name of the file without extension, so whatever you nominate as layer will be the file name, if you try to separate these by giving an explicit filename in dsn then the layer will be ignored.
u <- "http://geoportal.statistics.gov.uk/datasets/8edafbe3276d4b56aec60991cbddda50_2.geojson"
d <- rgdal::readOGR(u)
rgdal::...
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)
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(...
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")
...
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.
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 ...
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 <- ...
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,...
2
The R equivalent of the ST_SnapToGrid in PostGIS is in the lwgeom package:
# Snap to grid of 5000 m
lwgeom::st_snap_to_grid(x, 5000)
Works well to solve the non-noded intersection problem, and is quicker than applying a buffer of the same tolerance.
2
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
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 (...
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 ...
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
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
...
1
You can find the code from each chapter at https://github.com/Robinlovelace/geocompr/tree/master/code/chapters.
The code for the first animation is at https://github.com/Robinlovelace/geocompr/blob/master/code/08-urban-animation.R, and for the second one at https://github.com/Robinlovelace/geocompr/blob/master/code/08-usboundaries.R.
1
I don't think you need to create new objects if all you want to do is sample points in a raster according to a criteria. You can simply filter the cells you want using the which() function.
Using the raster library, you can call the values(r) function to get a list of values, or just use r[] (where r is the raster object.)
The which() function will tell ...
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
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 ...
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" "...
1
readOGR syntax can be a bit unfriendly at times. Here's a suggestion.
Assuming your shapefile path is - C:/Users/Merry/Desktop/ro_judete_polgion.shp, the following should work:
library(rgdal)
#try without specifying a layer
ro_judete_poligon <- readOGR("C:/Users/Merry/Desktop/ro_judete_poligon.shp")
1
Maybe I'm misunderstanding the problem, but to plot two sets of data using two different aesthetics and hence get two different legends in ggplot is done like this:
First some data, points will do:
> d1 = data.frame(x=runif(10),y=runif(10),z=runif(10))
> d2 = data.frame(x=runif(10),y=runif(10),z=runif(10))
> library(ggplot2)
Now make a plot of ...
1
You have an error in your code. SpatialPoints cannot be constructed from an atomic vector of coordinates. The following works fine:
library(sp)
library(rgdal)
download.file('http://download.massgis.digital.mass.gov/shapefiles/state/zipcodes_nt.zip',
f <- tempfile(), mode='wb')
unzip(f, exdir=tempdir())
p <- readOGR(file.path(tempdir(), '...
1
I found a solution on my own (thank you for your suggestion @spacedman).
The problem was due to the way I defined the basemap aestetics. Initially, I was defining the map aesthetics in ggplot()
basemap.110.gilbert <- ggplot(data = shape.110.tidy, aes(x=long, y=lat, group=group)) +
coord_map("gilbert", ylim = c(-60, 70)) +
geom_polygon(aes(x=long, ...
1
From the help file:
?gstat::idw
nmin = 0
nmax = Inf
maxdist = Inf
Only top voted, non community-wiki answers of a minimum length are eligible
