New answers tagged raster
0
This should work:
mystar = st_set_dimensions(mystar, 3, values = dates, names = "date")
where dates is the Date vector.
2
It sounds like you have Scale Dependent Visibility set on your raster layer. Have a look under the Rendering tab of the Layer Properties dialog of your raster layer and make sure that the check box for Scale Dependent Visibility is unchecked.
0
You can try the function 'blur' in the package 'spatstat' to apply a kernel based blur (Gaussian or otherwise) to your image. You can perform multiple iterations or adjust the parameters of the function until you achieve the desired result.
library(spatstat)
DEM_blurred <- blur(DEM, sigma=0.5, kernel="gaussian")
#adjust sigma value to change the ...
3
This error is resulting from an error check within the function.
if (length(dim(x)) > 2) {
warning("data must be grayscale image")
}
Given this condition length(dim(x)), a single band rasterLayer class object will always return 3 representing row, column and nlayers dimensions.
Whereas this function, in theory, operates on a matrix object it is not ...
0
The code above is incorrect. See this edit to avoid an 'expected raster or raster layer' error.
if i == 0:
out2 = out1
i += 1
else:
out2 = out2 + out1
0
If you have a polygon of interest you can:
1. Create a grid of points inside your polygon by using Regular Points from your toolbox and setting the input extent to your polygon's extent and the point spacing as you wish.
2. Take the values from your raster by using the Sample raster values with input point layer your regular points layer.
3. Add geometry ...
0
Think I solved it. If you add the WV2 data to a mosaicked dataset and use "Pansharpen and Multispectral" Processing Template it looks like it sharpens all of the bands
0
QGIS has a simple and effective plugin named "Serval" that I believe will help you out. To install Serval, select Plugins > Manage and Install Plugins > All > and scroll down to Serval, and click the Install Plugin button. A Serval toolbar section is thus created.
Using Serval is pretty straightforward. It has 3 modes:
Probing: Use this mode to change ...
1
If you are using ESRI products and you have point data with the proposed final elevation values of your project you can create a TIN from those points. Convert the TIN to a raster. Then replace that portion of the DEM with the raster version of your berm and pond. There is a lot of help to get you there online and these steps use standard tools.
Some ...
2
If you have very many raster files and you create overviews for individual images, then the rendering software must open lots of files for constructing a small scale output. In extreme case if you want to show the hole site all of the images must be opened, even just the smallest overview with not much data gets read from each file. I would say that opening ...
0
Try referencing the selected layer by name:
def get_statistic(self):
selected_layer = self.dialog.chooseCombo.currentText()
rlayer = QgsProject.instance().mapLayersByName(selected_layer)[0]
provider = rlayer.dataProvider()
5
Use Clip Raster by Mask tool from the Toolbox and click on the green iterator to clip the raster image by each smaller boundary inside the shapefile:
Input Raster:
Using the above tool with iterator enabled, here is the result:
0
There is a Python code for this on the ArcticDEM page:
"Strip DEM files are provided at 2-meter spatial resolution in 32-bit GeoTIFF format. Elevation units are meters and are referenced to the WGS84 ellipsoid. Strip DEM files include metadata text files describing the xyz offsets to filtered IceSAT altimetry data, although these translations have not been ...
0
You might be interested in rioxarray.
import rioxarray
xds = rioxarray.open_rasterio(...)
https://corteva.github.io/rioxarray/html/examples/reproject.html
xds_lonlat = xds.rio.reproject("epsg:4326")
https://corteva.github.io/rioxarray/html/examples/clip_geom.html
clipped = xds.rio.clip(features, features_crs)
0
You could use rioxarray https://corteva.github.io/rioxarray/html/examples/clip_geom.html
geometries = [
{
'type': 'Polygon',
'coordinates': [[
[12, 17],
[15, 17],
[15, 18],
[12, 18],
[12, 17],
]]
}
]
xds = rioxarray.open_rasterio(...)
clipped = xds.rio.clip(...
0
Create the mosaic with gdal_merge (documentation here):
gdal_merge.py -o merged.tif raster1.tif raster2.tif raster3.tif ...
If the first image (raster1.tif) have a resolution of 0.05m, the mosaic will have a resolution of 0.05m. If not, you can define the pixel size with the -ps parameter.
If there are two or more contributing pixels somewhere, the ...
7
Reprojecting rasters is usually a bad thing to do. It involves a non-reversible transformation from one grid system to another grid system that can have a non-linear relationship to the first. Hence the value in a cell of the new system can end up being some average of whichever grid cells in the source raster it overlapped.
If you have a raster and points ...
4
Yes, you can write a one bit raster with rasterio*.
You need to:
write to a format that supports a 1bit dataset, such as GeoTIFF;
ensure your numpy array is np.uint8/ubyte so rasterio doesnt raise the TypeError: invalid dtype: 'bool' exception; and
pass the NBITS=1 creation option to tell the underlying GDAL GeoTIFF driver to create a one bit file.
import ...
4
If you call rasterio.dtypes.check_dtype(np.bool_) you'll see that it's not a known dtype, because gdal doesn't support a true 1-bit dtype. GDT_Byte is the smallest. The list that rasterio is checking against is:
dtype_fwd = {
0: None, # GDT_Unknown
1: ubyte, # GDT_Byte
2: uint16, # GDT_UInt16
3: int16, # GDT_Int16
4:...
2
If you want to take an average of all 12 images, the process is very simple. Within a raster calculator (this can be the calculator in the raster drop-down menu or the GDAL/GRASS/SAGA calculators in the toolbox), create a tool-specific formula that describes an averaging equation:
(Pjan + Pfeb + ... + Pdec) / N
Where N is your number of observations (12).
...
0
I found a solution for my case:
I was using the given "minimum" value -3.40282e+38 from the Layers overview in QGIS.
However, using "Raster layer statistics" from the Processing Toolbox in QGIS (3.4.12), is showing more fractional digits: 'MIN': -3.4028234663852886e+38
Using this number to "add additional no data value" via the transparency tab in QGIS ...
0
If anyone is ever wondering the same thing, I deleted them and there were no issues.
1
We can carry out a small experiment that proves that the closest value is taken when reprojecting a raster by exporting it by default.
We can create a 5x5 pixel raster, with some simple values, and georeference it to a grid of 10 degrees of longitude by 10 degrees of latitude somewhere in the world.
Let's reproject on the fly to some system that ...
1
Yes you can pass an array. The documentation specifies:
rasterio.plot.show(source, etc...)
Parameters
source (array or dataset object opened in 'r' mode or Band
or tuple(dataset, bidx))
Yes, it's the same. Demo using rasterio.plot.adjust_band that show uses to do the adjustment:
import rasterio.plot as rp
import numpy as np
def norm(...
1
You can take you LAS files and add them to an LAS Dataset in ArcGIS. Then, use the LAS Dataset tool to filter the point cloud for different returns. Finally, use the LAS Dataset to Raster tool to create surfaces from your LiDAR data. See this link.
1
Mask raster function in R not working to extract values in one raster based on conditions in another
You have a mask with 0, 1, and NA values and you want to keep only the 1s. Hence you need to convert 0 and NA to the same value.
Here's a bit of a trick. Diving m by itself results in NA where m is 0 or NA (since 0/0 returns NaN, which acts like NA in mask) and 1 where its 1:
rmask = mask(r,m/m)
not sure if that's any quicker than the other solutions.
1
Mask raster function in R not working to extract values in one raster based on conditions in another
First change the values of 0 in wintPCP to NA:
wintPCP[wintPCP == 0] <- NA
then mask:
masked_ras <- mask(MAT2resampled, wintPCP)
0
For those of you interested in the future, divinding by the quantification value as suggested by @Zac Wang produces the right result
Another option is to change the data type from uint16 to float. Using the following code, the issue is solved:
red = rasterio.open(S_files[2]).read().astype('float')
nir = rasterio.open(S_files[3]).read().astype('float')
ndvi ...
0
In case someone is interested in the future, I managed to make it work with the following code.
First, open an empty list. Then import the shp with geopandas. Next, we access the geometric information and extract it. The output of that is a dict containing the information required by rasterio.mask.mask
geom = []
StudyA = gpd.read_file('.../...
0
Using no data value -9999 the black background was removed. Thank you for the help.
0
So I have currently done this by just making a JOIN on the raster geometry and this seems to work reasonably well (given the size of the rasters). If I am missing something basic and there is a much faster way to do this it would be good to know:
-- this works and seems to get us correct values
WITH solarpot AS (
SELECT fid,
(ST_Intersection(T....
0
Rearrange bands tool (QGIS version 3.4 +) can delete the 4th band easily.
Go to the Processing Toolbox: GDAL > Raster conversion and start Rearrange bands tool.
Select your multiband (RGBI) raster as the Input layer.
Click on a small three-dots icon at the right of Selected band(s) option.
Select first three bands (Bands 1, 2, 3) and hit OK, then Run to ...
0
I got a similar weird result when I used complete jp2 images (with DN -digital numbers- values without transforming in reflectances) for calculating NDVI index for a Spanish region.
import rasterio
b8 = '/home/zeito/pyqgis_data/jp2/T30STF_20170422T110651_B08.jp2'
b4 = '/home/zeito/pyqgis_data/jp2/T30STF_20170422T110651_B04.jp2'
NIR = rasterio.open(b8).read(...
1
After using Python read the data, you may have ignored the no data value such as 0 or -999 from the numpy array.
When no data value gets involved in the calculation, it could affect the result but I think QGIS will ignore no data value while calculating.
0
This thread is a little old. I thought it would be good to share my method for solving the issue presented. I have GEOS 3.8.0 installed so the ST_Union is pretty snappy. I haven't tried testing it on large datasets.
I use a function to determine the bounding tiles, merge those into one raster, and then apply a mask with a buffer returning the masked raster ...
0
You can find the solution here. For my opinion the best solution is second. Work on vactor or raster layer.
Under Vector menu go to Research tools and you will find polygon from layer extent. Also under processing you will find that function.
In QGIS 3 use ctrl+k and start typing extent.
0
In QGIS, you can use the Tile Index tool from:
Raster / Miscellaneous / Tile Index
1
Try converting it using the Raster To Polygon tool which:
Converts a raster dataset to polygon features.
then you can use it as a Mask.
2
First convert the rasters to variables in the same dataframe, then calculate the pairwise correlations and use the package 'corrplot' to display the results in a matrix.
#dependencies
library(raster)
library(corrplot)
#read in rasters
r1 <- raster("IMG_0003_1.tif")
r2 <- raster("IMG_0003_2.tif")
r3 <- raster("IMG_0003_3.tif")
r4 <- raster("...
2
The result is correct. Since the extent of A is larger, the function is correctly cropping the extent of B. To get the same "shape", try something along the lines of:
mask(B, crop(B, extent(A)), A)
However, this will still show a smaller extent in B. If you really need the two rasters to show the same extent of data then you will likely have to coerce ...
1
In QGIS you could use the "Extract layer extent" processing tool to get a bounding box for Raster B. Then use that as the mask layer in the "Clip Raster by Mask Layer" tool to clip Raster A
0
You may be prescribing the new height and width incorrectly. What happens when you call out_image.shape before writing the raster? It may be the shape is actually 1 x 6094 x 4000 (or whatever the width is) and so when you prescribe shape[0] as the height you are actually giving it a height of 1. The value of 1 in the zeroth element of the shape tuple is ...
0
You do not specify the software?
In arcgis the "From TIN" toolset will convert the tin to another geometry from which you can get the attribute table.
While it is possible to "Build a Raster Attribute Table" with the tool of the same name in ArcGIS this is for integer format tables because these are likely to have a limited number of values for each of ...
0
It's a bit of a pain, but you need to write the resampled numpy array to a rasterio Dataset (either on file or in memory) and adjust the transform to match the resampling..
Here's an example of both:
# Example licensed under cc by-sa 4.0 with attribution required
from contextlib import contextmanager
import rasterio
from rasterio import Affine, ...
1
You'll want to use the translation, rotation, and scale tools to transform the data manually. For translation use the "shift" tool, for rotation use "Rotate" and for scaling use "rescale".
1
It is possible by using windowed reading and writing available in rasterio. First, you need to find each raster bounding box and intersects them by using shapely python module. Intersecting bounding box is used for calculating indices of column, row for first and last pixel in each raster. So, these values are used for windowed reading both raster as array ...
1
You can use Zonal Statistics tool to calculate the average raster cells within polygon grid. Just search for Zonal Statistics and you will find it as in the image below.
0
I tried this solution of mosaicing to a folder instead of FGDB, it worked for some of my projects but still had gaps in others. After further experimentation, I found that recalculating the pyramids to 0 for all the raster tiles first, then mosaic to new raster, was successful in creating the full mosaic without gaps. Possibly a combination of this plus ...
2
Use the SAGA tool Clip raster with polygon. This tool can be found in the Processing Toolbox.
Click the round green arrow button next to the Polygons layer. This button toggles the option to
Iterate over this layer, creating a separate output for every feature in the layer.
0
I was able to remove that "noise" by using the "Clip (Data Management)" tool.
Top 50 recent answers are included
