New answers tagged raster
0
Let's say your DEM layers have the names DEM1 and DEM2.
The following expression (in the QGIS Raster Calculator) returns 0 for the differences that have an absolute value that is less than 0.005, and the calculated difference larger absolute differences:
"DEM1@1" - "DEM1@1" * (("DEM1@1" - "DEM1@1" > 0.005) OR ("DEM1@1" - "DEM1@1" < -0.005))
0
Please update your installation of the GRASS GIS addon r.recode.attr to the latest version. I have updated it to be Python 3 compliant (https://github.com/OSGeo/grass-addons/pull/91 - xrange compatibility fix).
The addon now works.
2
Alum.DESKTOP-EEJENVF looks like a bad choice for an output path. The dot in the folder name may throw off some tools. Try something simple like C:\temp\cut1.tif instead.
1
Your indentation is off, but I'll assume that's just a copy/paste issue.
Your issue is that 'F:\\prep\\converted\\*.tif' will try to output a file called literally *.tif which won't work. You need to provide an actual file name. I suggest using os.path.join and os.path.basename
Try this:
from osgeo import gdal
import glob
import os
input_path = 'F:\\...
0
From the author's comment:
I have figured this one out, if you remove the line
pipe.set(renderer.clone()) , it will make an exact duplicate of your
DTM and not as a 4-banded rendered image. Hope this is of some use to
anyone trying to batch save raster tiles like i was.
I can also confirm that this is the solution for the issue where single band ...
3
The h.binCount is a fixed value that represents the length of h.histogramVector. In other words, it is the number of possible classes of grouped values between h.minimum and h.maximum.
It can be observed with following code and three examples.
l = iface.activeLayer()
p = l.dataProvider()
#regardless third parameter
p.initHistogram(QgsRasterHistogram(),1,...
1
Without more information on the CRS of your coordinates, and their construction., it will be difficult to create a working solution, but I'll make assumptions along the way and see if you can adapt it to your solution. What you'll want to do is map a function to your list that will take each element of the list (the polygon coordinates) and make a geometry ...
0
We cannot see what CPC_rename and folder7 are set to, and consequently it is possible that CPC_lines is an ill-formed or otherwise inappropriate workspace name.
I believe that I have often seen ERROR 999999 when the name of a workspace being written to is either ill-formed or otherwise inappropriate.
1
You need to export the dates variable as well because it is a global
beginCluster()
z1 <- clusterR(ndvi, calc, args=list(fun=f), export=c('a','dates'))
endCluster()
or
c1 <- getCluster()
clusterExport(c1, c("a","dates"))
z2 <- calc(ndvi, fun=function(x){ t(parApply(c1, x, 1, f))} )
0
British Geological Survey did an OpenLayers WCS 2.0.1 demonstrator, you can see it in operation at:
http://ogcdev.bgs.ac.uk/ogcclient/WCS/GetCoverage_v2_0_1.html
Also note the OneGeology portal (which uses OpenLayers) supports WCS services.
2
Take a look at np.where and maybe do something like:
y_true= np.ravel(img_r)
y_pred= np.ravel(img_pre_r)
y_true = y_true.astype('int')
y_pred = y_pred.astype('int')
indx = np.where((y_true != -9999) & (y_true != 205) & (y_true != 210) & (y_true != 215))
y_true = y_true[indx]
y_pred = y_pred[indx]
Or apply the np.where call to both true and ...
0
I know it's a bit late but hopefully this answer helps someone else. The raster layer in contained within the Geodatabase file provided by NRCS. However, at the time of me writing this QGIS cannot open raster files contained within Geodatabases. It has to do with ESRI's API. Here's more on that topic Opening Raster stored in Esri Geodatabase using QGIS
...
0
I had the same problem with my DEM merging.
Finally, I launched the QGIS desktop 3.8.3, where everything was fine.
Alternatively, you can use SAGA for it. Both ways are under the link below:
QGIS Digital Elevation Model merging problem
0
While I can't test or verify without information on how the original data was processed, my guess is that there are default CRS values that differ slightly in QGIS from the software where your original raster was created.
I've had issues between FME and ArcGIS before where reprojections of the same data using the same settings yielded very slightly ...
0
This sounds like the Raster Calculator's output extent, resolution, and CRS do not match your input layer's. QGIS is therefore carefully resampling the input layer to generate the output layer with the "wrong" lattice structure.
Make sure you invoke the the Raster Calculator with the input layer, in your case bathymetry, selected; not another layer with ...
0
When you are using the Raster Calculator tool, make sure you are setting the Cell size, Output extent, and Output CRS to the same as the reference layer. If you are unsure of the cell size of the reference layer, you can leave that blank, but make sure to set the other two settings to the same as the reference layer. This should provide you a 1:1 match ...
1
At Arc Desktop 10.7.1, this works for me (note that it requires the Spatial Analyst extension):
Open the Raster Calculator and enter the following in the Map Algebra Expression pane:
Con("your_raster" == "your_raster".maximum, "your_raster")
[Note: this selects the cell with the maximum value from your_raster.]
Enter the name of the output raster that ...
3
Not sure if transforming vector is worth it, but here's the trick.
According to docs, .transform() is ony applied to an ee.Feature, not ee.FeatureCollection. So you have to map the procedure over a FeatureCollection.
// Define an arbitrary region in which to compute random points.
var region = ee.Geometry.Rectangle(-119.224, 34.669, -99.536, 50.064);
// ...
2
You can choose which bands to display under Symbology. Either singleband or multiband (RGB):
You can add the layer many times to table of contents to display all bands on top of eachother.
1
This is quite easy with gdalwarp (2.4 and newer):
gdalwarp -of GPKG inputfolder\*.tif output.gpkg
gdaladdo -r cubic output.gpkg
If your source data is single band 8 bit grayscale, then the tiles in the GeoPackage will be likewise.
However, the GDAL GeoPackage driver will by default expose a GeoPackage dataset as four bands RGBA. This means that even ...
1
Check what gdalinfo says about the image. Notice that you did not get an error but just a warning and it says Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel.
The warning means that the image has, for example, 4 samples for each pixel but the Protometric TIFF tag is set to RGB and there is no ExtraSamples tag ...
0
not sure if this is for GISN08 but I hit the same error on one of this module.
The way to resolve this is to put in the full projection of the coordinate system for input dataset as an argument to ProjectRaster_management:
arcpy.ProjectRaster_management (ifpn, ofpn, "4326", "NEAREST", "#", "#", "#", "PROJCS['...
1
You're missing just a little bit of detail in your screenshots that makes me think you're going to have some other issues crop up.
Can you please let us know how the symbology is being assigned in the first (working) example you gave. Is the field being symbolised the "Label" field, or something else?
For example, here I am symbolising DonutType
In your ...
2
If the focal value is not specifically indexed in the function it looks like the focal function is assigning the incorrect values and positions back to the matrix. However, this is expected behavior. For these types of operations, you have to specifically index the focal value to get the correct results. One would think that the focal value would ...
0
Looks like your raster datatype may be set to integer.
Try: rasterOptions(datatype='FLT4S')
2
You can set in the layer styling the Layer Rendering Blending Mode to Multiply for the coloured layer:
1
I took another approach as I could not figure out the above solution. The above is probably easy but somehow I failed.
My approach was:
Perform a union with [Polygon] and [Study_Area]. Output: [Polygon_StudyArea_Union].
Edit values of the polygon in [Polygon_StudyArea_Union] to 0 and the resulting to 1.
Convert [Polygon_StudyArea_Union] to raster. Output: [...
0
I solved it by another way with earthpy library. You'll need a dataframe with daily and hourly azimut and altitude.
Here's the script:
def hillshade (index):
altitude = df.altitude[index]
azimut = df.azimut[index]
return ep.plot_bands(
es.hillshade(elevation, azimuth = azimut, altitude = altitude),
)
plt.show()
for i in df.index:
...
1
A little more information would prove useful to help answer your question. What language are you using, and what function are you using to import your data? I will base my answer on the assumption you are using R with the raster library and the raster function to import your .tif files. If you want your tif files to come in as a stack, use the raster::stack()...
1
This is a problem with handling EPSG:3857 in GDAL 3.X. You can run transformations with gdal 2.X and geoserver will read this file.
For example, you can use gdal_translate from folder "C:\OSGeo4W64\apps\gdal2\bin" if you use qgis from OSGeo4W.
Detail information on Geoserver`s JIRA:
https://osgeo-org.atlassian.net/browse/GEOS-9475
3
Yes, of course.
Actually the power of GRASS really shines when you need to perform repetitive tasks in a loop.
The specific answer depends on how your CSV is formatted, what operating system you're working on, and what command shell you prefer. For example:
If I had a CSV that looks like:
day,month,year,hour,min
01,01,2019,10,30
02,01,2019,11,30
03,01,2019,...
2
Best thing to do with questions here is to make a representative simple example.
Here's a 3x4 raster with some values:
r = raster(matrix(c(1,1,1,2,1,3,55,12,3,3,2,1),3,4))
Now build a data frame - find the unique ID values:
ids = unique(r[])
ids
# [1] 1 2 55 3 12
And make a data frame like yours with a couple of data columns:
d = data.frame(ID=...
0
I finally solved the problem by using numpy.memmap to
create a memory-map to an array stored in a binary file on disk
and then processing the input rasters in windows and blocks. It might be slower and but it works and I'm happy with the result (need to thank user @Thomas that helped me in some steps).
The code I am using is taken and modified from the ...
0
Sometimes parallel is not the way to go. If the function is not written for it you get all the overheads and none of the benefits.
Often, just reducing the demand for memory allocation solves performance issues.
I have found that if the number of points is high it is much quicker to hit the raster layers one at a time (or in smaller stacks) rather than ...
1
You could use exactextractr::exact_extract for this. It will work with multi-part polygon inputs.
library(raster)
library(sf)
library(exactextractr)
semana21 <- raster("Semana21.tif")
poly <- st_read("shapefile.shp")
poly$sum_semana21 <- exact_extract(semana21, poly, 'sum')
Note that these results will differ from raster::extract because the ...
2
You can get a table of cell values, coverage fractions, and center coordinates with the include_xy argument to exactextractr::exact_extract. Here's an example:
library(raster)
library(sf)
library(exactextractr)
# Pull municipal boundaries for Brazil
brazil <- st_as_sf(getData('GADM', country='BRA', level=2))[1:10, ]
# Pull gridded precipitation data
...
2
The description of the -a_nodata parameter from the documentation of gdal_translate https://gdal.org/programs/gdal_translate.html
-a_nodata [value]
Assign a specified nodata value to output bands. It can be set to none to avoid setting a nodata value to the
output file if one exists for the source file. Note that, if the input
dataset has a ...
2
For me, the clue/solution was in the metadata linked here:
http://www.cacpd.org.s3-website-us-west-2.amazonaws.com/climate_normals/NA_ReadMe.txt
Coordinate Systems
------------------
X,Y:
Projection: Lambert Conformal Conic
False Easting: 0.00000000
False Northing: 0.00000000
Central Meridian: -95.00000000
Standard Parallel 1: 49.00000000
Standard ...
1
Thank you!
But i solved my problem in google engine. here is my solution.
// create new raster with the required information
var sd_inputraster = inputraster.reduce(ee.Reducer.stdDev());
var mean_inputraster = inputraster.reduce(ee.Reducer.mean());
var vco_inputraster = inputraster.expression('vco = sd / x',{'sd':sd_inputraster,'x':mean_inputraster});
//...
3
you can use the raster::xyFromCell function to get a matrix of coordinates from a vector of cell numbers, which you can cbind onto the extract (as the xyFromCell result is in exactly the same order as the vector provided)
I've used sf but of course this works with readOGR etc.;
require(raster)
require(sf)
# make 2 dummy rasters for a dummy stack
r <- ...
0
any advice that can be offered?
According to the Circuitscape web site
Oct 2013- Circuitscape can now be called from an ArcGIS toolbox. No more converting input grids to ASCII format! Just install the ArcGIS toolbox from our downloads page when you install Circuitscape.
So you don't need to convert to Arc/Info ASCII Grid at all.
1
It seems that there is a problem with the elevation values in this tile and other tiles in that area of San Bernardino. I got the same issue as you describe when opening the tile, the elevation values are way off. However, tiles in other projects/cities appear perfectly fine.
If you divide the values by 0.3048 with raster calculator, the values match ...
1
The window describes the relative location of the block within your file.
If CRS of both files are identical, you can use the bounds of the input window to compute its relative location within output window
import numpy as np
import rasterio
from rasterio import Affine
from rasterio.windows import bounds, from_bounds
# adopt as needed
input_file = "/some/...
1
.asc is definitely there for you to choose.
0
You can use numpy and scipy for that: the following script shows the basic principle how to access rasterdata with numpy via gdal:
import numpy as np
from scipy import stats
from osgeo import gdal
'Active Layer must be a raster layer
raster = iface.activeLayer()
rl = raster.dataProvider().dataSourceUri()
ds = gdal.Open(rl) #open raster with gdal
data = np....
4
(using QGIS 3.8.2)
I was able to merge 4 adjacent tiles from the link provided in the comment from @Erik using:
Warp (gdalwarp), in the QGIS toolbox, outputting them to .tif
Merge (gdal_merge), in the toolbox again.
I didn't have to change any settings, I left everything as default.
3
this is not an answer to your question, but may solve your problem anyways. the rasters that you are trying to mosaic are also available as a WCS under this link: https://isk.geobasis-bb.de/ows/dgm_wcs
You can just add the WCS to QGIS and use it for your processing.
1
I can think of two solutions:
Dump your raster as polygons using this function ST_PixelAsPolygons(). Then, do your calculations for the centroid of each polygon (first calculate the centroids of the polygons). Now, you should have a table of points that each of them has a value that you calculated. Finally, build an array of your points (a geomval object) ...
3
Your images have a positive N-S pixel resolution in the geotransform. This is
quite unusual and means the origin coordinates for the image are in the
lower-left rather than the upper-left of the image. In the usual upper-left
case the geotransform N-S pixel resolution would be negative.
Are you sure your image georeferences are correct?
Try exporting ...
1
You may try and convert those tiles to a virtual raster. Export this virtual raster then into a new .xyz file and you should be good!
Top 50 recent answers are included
