New answers tagged raster
4
You can use Raster - Projections - Warp and set output resolution to 10:
2
You are wrong to assume that 1 grid cell (pixel) side is 1 meter in UTM or any other coordinate reference system. It can be anything (up to 60 km if you want to stay within UTM zone boundaries). With a longitude/latitude coordinate reference system (I assume that is what you erroneously refer to with WGS84) cell size is not constant, but for a very small ...
0
This is a fairly common occurrence. You should be able to fix is quite easily by running the "Fill Sinks" tool in SAGA. There are several ones to chose from. It literally just lifts low points in your DEM to allow for smooth flow of water. The resultant DEM is "hydrologically correct"
3
CRU is mm/month, ERA is m/day. So it depends on the number of days in the month. But for January you could divide CRU by 31000 or multiply ERA with 31000
Say you had an mean rainfall of 10 mm/day, ERA would report 0.01 and CRU would report 310 (for a month with 31 days).
310 / 31000
#[1] 0.01
0.01 * 31000
#[1] 310
As a function perhaps like this
era_to_cru &...
2
In the case of DN's (Digital Numbers) "rescaling" is not a simple transformation but, represents converting DN's to reflectance, representing a 0-1 floating point. Take a look at the landsat handbook or here for the equations for converting DN to reflectance.
For reflectance that is stored in a 16 bit range, you apply a scaling factor to return the ...
-1
You do not provide much information making it difficult to help. But I do know that the terra package handles time better than raster. So I would suggest trying the below
library(terra)
shp <- vect("path/shapefile.shp")
temp <- rast("path/temperature70-78.nc")
nc.mask <- mask(temp, shp)
writeCDF(nc.mask, filename = "path....
0
If you want to extract the values without having 0's, you need to use a different formula in the Raster Calculator.
"img1@1"/("img1@1">0)
0 values will be returned as no values.
0
(("Landcover_raster"=9)+("Landcover_raster"=10))*.21
should work. It multiplies X by .21 where X is either 0 or 1, depending on whether it is land cover (9 or 10).
It is helpful to remember that with 1 and 0,
OR is +,
AND is *
0
I am new-bee, my answer is welcome to be corrected.
Seems in ArcGis, you can open the aux file in 2 ways:
1. if your folder looks like:
in this case, the aux file acts like an entrance of the bundled-files, you can open this via:
step1. open ArcMap ( version > 10 ? seems this is a requirement to show the correct raster pyramids)
step2. Ctrl + N to open a ...
0
You will know from signal analysis that this is a very big convolution in the spatial domain. Can you transform to a wavenumber (frequency) domain? Then it's just a multiplication and probably much faster.
You could use gdal to read the data into a 2D raster, then
use Numpy to get the 2D Fourier transform of the data, and also of your filter (padded to the ...
0
I haven't figured out why GEE is coding masked pixels as NA when using the mean() reducer and as 0 when using the max() reducer but a work around can be found in this thread: https://stackoverflow.com/questions/36960974/how-to-replace-raster-values-less-than-0-to-na-in-r-code/49159943.
I edited my code to set the masked values to 20000 using the unmask() ...
0
Recently came across this and it ended up helping me, so I thought I would share an update.
Adding additional arguments help as well as adding raster constraints to the table you are trying to use.
ds = gdal.Open('PG:host='' dbname='' schema='' table='' user='' port='' password='' column='' mode=2')
SELECT AddRasterConstraints('raster_table'::name, '...
1
You can
use the join field tool to join your DBF to the raster attribute table, then
use the tabulate area tool with shapefile / country as the in_zone_data / zone_field and raster / type as in_class_data / class_field
5
You should use QgsHueSaturationFilter class:
huesat = QgsHueSaturationFilter()
huesat.setSaturation(30)
layer.pipe().set(huesat)
0
Consider trying the QGIS > Raster > Conversion > Polygonize (Raster to Vector) tool. Where adjacent raster cells contain the same value (in your example, the value 5), this tool dissolves out their shared boundary, resulting in fewer output polygons. Taken to its logical extreme, if all of the raster cells contain the same value, the vector output ...
1
I only see a small difference between raster and terra. That is not unexpected as the extent is not exactly defined when projecting a raster based on a different CRS alone. It is better practice to provide a raster template to project to. Here is what I see:
options(rgdal_show_exportToProj4_warnings = "none")
library(raster)
library(terra)
dtm &...
3
One option is using "Raster pixels to polygons" from the Processing Toolbox. It creates a new vector layer containing a polygon for each pixel. If you remove the symbol fill, you will get the boundary:
0
Try:
cpr_before = next(x for val, x in enumerate(change_point_raster)
if val < change_point_index)
Instead of
cpr_before = next(x for x, val in enumerate(change_point_raster)
if val < change_point_index)
Ithink you want to return an array based on an index, not the other way around. Same for ...
0
Install the geodata package to download some WorldClim data
# install.packages("remotes")
remotes::install_github("rspatial/geodata", dependencies=FALSE)
library(terra)
library(geodata)
wc <- worldclim_country("Iceland", "tavg", ".")
A fast and simple approach
bio1 <- mean( clamp(wc, 0, Inf) )
A ...
1
You can do mathematics with rasters with the usual mathematical operators.
There are assorted functions in the raster package for doing grouped functions, like monthly averages over stack of rasters.
There's also the terra package which works with rasters and should be faster than raster.
4
I think that you may want the 95% Confidence Interval and not the Standard Deviation, which is not terribly informative.
The Confidence Interval is the standard error * critical value for the T distribution P( t > cv ) = 0.025 = P( t < -cv )
Your example data was not quite large enough to effectively demonstrate a statistical solution. Sans plotting ...
1
You should not use nc_open. Just do
library(raster)
tmp1 <- brick("cru_ts_3_00.1901.2006.tmp.nc", varname="tmp")
Now subset
start <- 12*(1960-1901) + 1
end <- start + 30 * 12 - 1
start
#[1] 709
end
#[1] 1068
tmp <- tmp1[[start:end]]
But I would use terra like this (I am using a more recent version of the database):
library(...
3
106 years is exactly 12 * 106 = 1272 months, so that's the 1272. If it is in a sensible order (ie Jan 1901 to Dec 1901, then 1902, etc etc, a short formula can work out the numerical layer indices for a given year:
year.layers = function(year){
start = 1 + (year - 1901)*12
start:(start+11)
}
Quick test, it should return numbers from 1 to 12 for 1901, ...
2
This should work better with terra.
library(terra)
#terra version 1.2.4
r <- rast(nrow=10, ncol=10)
r[] <- 0
r[51:100] <- 1
r[3:6, 1:5] <- 2
r[1, 1] <- 3
coltab(r) <- c("#FF0000", "#FF9900" ,"#99FF00","#0000FF")
levels(r) <- c('Forest','Water body','City','Cropland')
x <- writeRaster(r, "...
0
After coming back to this problem, I was able to determine a cause, though I am not sure why. In the map properties I adjusted the Opacity to 50% and saved, it worked. A box will appear letting you know there will be some composition effects. I do not believe this causes any issues with GeoReferencing, but I could be wrong. I then changed back to 100% and it ...
0
Try this:
## 10-days mean
import pandas as pd
# get data
df = pd.DataFrame([{'lat': '82', 'long': '13', '20210102': 0.02, '20210103': 0.06, '20210104': 0.08, '20210105': 0.09, '20210106': 0.07, '20210107': 0.06, '20210108': 0.2, '20210109': 0.07, '20210110': 0.03, '20210111': 0.04, '20210112': 0.08}])
# transpose dates
df = df.melt(id_vars=["lat",...
3
Your stretch is set to custom. Stdev is not relevant to custom stretch so is disabled (you adjust the histogram manually).
Set Primary symbology to Stretch, change Stretch Type to Standard Deviation, change Statistics to Custom, specify your desired statistics:
1
That seems to be a low level error from libtiff library. By the error message the libtiff function cannot read data from a certain part of the TIFF file. Directory here means the internal directory that TIFF file has and offset means the number of bytes from the beginning of the file.
Either your copy of the TIFF file is corrupted or you have a different ...
0
Zeroing will disturb the mean.
Maybe a more complex solution, case-by-case is called for. If nodata is -9999
Cases
There are 7 cases:
1 case where all 3 exist
3 cases where 1 doesn't but 2 do exist
3 cases where 1 does but 2 don't exist
Test if all 3 exist. * acts like a logical AND if arguments are 1 or 0:
(A>-9999) * (B>-9999) * (C>-9999)
...
0
As user2856 commented, when opening a new file for writing you can just add "crs" argument in "rasterio.open":
with rasterio.open(f'/.../{lokalnyid}'+'_025'+'.tif', "w", crs='EPSG:2180', **out_meta) as dest:
dest.write(out_image)
An example from rasterio documentation:
profile = {'driver': 'GTiff', 'height': 100, 'width': ...
2
You can use terra like this
library(terra)
# terra version 1.2.5
x <- rast("3B43.19980101.7.HDF")
#Warning message:
#[rast] unknown extent
x
#class : SpatRaster
#dimensions : 1440, 400, 3 (nrow, ncol, nlyr)
#resolution : 0.0025, 0.0006944444 (x, y)
#extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +...
0
Depending on your raster data type, you can store them in different bands.
solution is here.
For instance, if you have single band information such as elevation, temperature, slope and you want to merge them in one layer. You can allocate each band to each layer.
Make sure you have the right naming, and also store the band information in a text file so it ...
1
Here is an approach with terra
Example data
library(terra)
#terra version 1.2.5
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
set.seed(1234)
pts <- spatSample(r, 2, xy=TRUE, values=FALSE, as.points=TRUE, na.rm=TRUE)
Solution:
b <- buffer(pts, 5000)
cnrs <- cells(r, b)
xy <- xyFromCell(r, cnrs[,2])
plot(...
3
If you are going to take your data to R, you can reclassify them there.
Example data:
library(raster)
r <- raster()
values(r) <- runif(ncell(r))
Solution
threshold <- 0.5
x <- r > threshold
x
#class : RasterLayer
#dimensions : 180, 360, 64800 (nrow, ncol, ncell)
#resolution : 1, 1 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ...
1
It seems that the slowest part of the code is the call to Moran. I am using raster::extent in lue of sf::st_bbox which seems to make a difference. I also moved the for loop into an lapply.
Let's add libraries and create some example data. I would recommend reading the data into R as a stack, and index out of the stack, rather than repeatedly reading a raster ...
3
I have a dem and some points:
> plot(dem)
> plot(pts, add=TRUE)
I'll make a 1km buffer for my data:
> d = 1000
> buf = st_buffer(pts, dist=d)
Then I can get all the grid points in the buffer by masking the DEM to the buffer and taking all the coordinates of non-NA pixels and making an sf points data frame out of them:
> gridpts = st_as_sf(
...
1
Run Zonal Statistics As Table (or the QGIS equivalent) on your vector layer to get the min and max per polygon. Select polygons where min is not equal to the max (i.e. these polygons cover multiple raster values). Clip the raster based on these polygons, and then polygonize this clipped raster. Merge the resulting polygons back with the original polygons ...
1
I would recommend to use Reclassify by table tool which does the same as r.class tool. However, it is more manageable as it gives you min, max, and value in an organized way to avoid any incorrect syntax when you write it as a text.
Select the input raster layer
Select the table and write the minimum, maximum, and new value. The default is min < value &...
1
Looks like mosaic.tif is compressed while smoothing.tif is not:
Image Structure Metadata: COMPRESSION=LZW
I don't have any experience with orfeo, but it seems you can pass GDAL creation options to create compressed outputs by using an "extended filename", see:
https://www.orfeo-toolbox.org/CookBook-6.6.1/AdvancedUse.html#examples
https://gdal.org/...
0
It is just a matter of writing a function that returns the proportions and using lapply to apply it to the list returned from exact_extract.
Add libraries and create example data
library(sf)
library(raster)
library(exactextractr)
nc <- st_read(system.file("shape/nc.shp", package="sf"))
nc <- st_cast(nc, "POLYGON")
r &...
0
Figured it out.
First, reclassify raster making agricultural values 1 and everything else 0:
agri_land_use_ras <- raster::reclassify(lu_raster, lu_matrix)
See here for how to set up a reclassification matrix.
Then, can just take the mean of the resulting raster:
Results <- exact_extract(agri_land_use_ras, counties, 'mean')
0
Please have a look at zonal statistics and zonal statistics as table tools.
Given the objective, you might want to use Zonal Statistics as Table tool. It will generate a table as output (in contrast to the Zonal Statistics which returns a raster as output). Once you get the table, you may simply join it to your input shapefile to get population for each ...
1
This:
raster::terrain(input$DEM$datapath, opt = 'aspect')
looks like you are calling terrain on the file path - you need to read it into a raster first. Perhaps:
raster::terrain(raster::raster(input$DEM$datapath), opt = 'aspect')
which doesn't remove the first error, which is because the field is empty. Once you upload a DEM you will get a plot:
you ...
0
I know this question is history now but in case anybody else comes by looking for help:
I'm not sure why to be honest, but running r.grow and r.reclass and then r.thin seems to work to prepare rasters for r.to.vect conversion.
The code should look something like:
# reclassify raster, mapping raster values that you want to turn into lines to 1, others to NULL
...
0
This is the additional code I added to my original code posted above
// projection of final raster: 250m is the resolution of the fire data
var reproj_params = {crs:'EPSG:4326',scale:250};
//map the reprojected new_ds fire image collection over the peatlands to see peatland fires only
var burntPeat = new_ds.map(function(image) {
return image.multiply(...
0
Create a difference layer with the Raster Calculator, to get the difference in glacier heights over time, ie. Raster2 - Raster1. This will create a new raster layer
Use Zonal Statistics on the new layer to compute the mean difference within each polygon.
Let us know if you get stuck.
1
The very first thing I would try is to bring your new TIN into a blank ArcMap session. The current data frame coordinate system may be the reason your TIN doesn't 'line up' with the imagery.
You may have done exactly this but these are the steps I need to take:
Copy a standard esri projection file (.prj) to a new location. For example if your data is in ...
0
I had the same problem a few years ago. Vector data came in fine but the LandXML was off. I went through some old notes and found how I got it to work, but I had to do it in Civil 3D. Maybe they could modify it for you. The scale will be different I'm sure.
In Civ3D:
Select the surface/contour (use QSELECT to query)
Type SCALE, enter
Click the surface/...
0
Replace "multiply an image collection by an image" with "multiply each image in the image collection by an image" and you already have the tool you need: .map().
var burnedPeat = burnedArea.map(function (burnImage) {
return burnImage.multiply(onlyPeat);
});
3
It must be the mean SCE for each pixel over the 10-day period.
Raster values may or may not indicate snow (they could be land/water etc.)
An approach could be:
Remove all non-snow values (water, land, etc). This can be done by Raster Calculator with an expression:
(raster1@1 > 100) * raster1@1
That should change any value below 101 to 0. 100 means "...
Top 50 recent answers are included
