New answers tagged raster
0
Let's make a reproducible example...
First use the raster::getData package to get a raster of the heights of Liechtenstein (because Liechtenstein is quite small so we can see what its doing, and also because I played in a band that had a number 3 hit in the country...):
> library(raster)
> library(sp)
> raster1 = getData('alt', country='LIE', mask=...
2
First step is to visualize tiles and their unique ID:
SELECT rid, rast::geometry
FROM yourRasterTable;
Then you can visualise specific tiles and group of pixel values:
SELECT (gv).geom, (gv).val
FROM (SELECT ST_DumpAsPolygons(rast) gv
FROM yourRasterTable
WHERE rid = XX) foo;
You can replace ST_DumpAsPolygons() with ST_PixelAsPolygons() if ...
1
In the Rasterize tool, notice the option to set a nodata value:
Assign a specified nodata value to output bands [optional]
You have 0.000000 set as the optional nodata value, so all polygons with 0 in their burnin field are converted to nodata values in the raster.
Clear the nodata value by clicking the X button in the nodata box:
1
Note the difference between band rendering done by QGIS and the actual raster extent and values. The values of your raster remain the same, regardless of the band rendering (Multicolor band, Singleband gray, Paletted unique/values) you select in Layer properties.
If you take your original layer (i.e. HDF4_EOS:EOS_GRID:"MCD19A2.A2018312.h24v06.006....
0
I don't think that there is a direct way of doing it (maybe with mapAlgebra?), but an easy way would be that once you get the value of your pixel, you should try using ST_PixelOfValue, then use the .x and .y and pass it to ST_RasterToWorldCoord (remember it gives you the upper left point of the pixel, not the center, I don't really understand why). That ...
0
Use "Export>Save as" from raster do you want.
Change format type to "Geopackage" and select the path and name of new gpkg file that will be created, and select tha name of Layer.
If you select the existing GPKG, be careful, could will overwrite and your will miss all data store.
1
I had the same error with QGIS 3, by writing multiple times to a vector layer with "qgis:zonalstatistics" processing function.
The layer wasn't in the registry:
registryLayers = QgsProject.instance().mapLayers().keys()
legendLayers = [ layer.layerId() for layer in QgsProject.instance().layerTreeRoot().findLayers() ]
layersToRemove = set( registryLayers ) - ...
1
OpenJUMP does not support PostGIS rasters at all as normal raster images. However, it can be used for visualizing vector data that is returned in the Well Known Binary (WKB) format by the PostGIS Raster function ST_AsBinary https://postgis.net/docs/RT_ST_AsBinary.html.
There is a tutorial about how OpenJUMP can be used for visualizing PostGIS raster data ...
2
The first warning means that the value of GeoASCIIParams tag is not read as it was written because the original image is having NULL character in the value of the tag. NULL can be used as a delimiter between strings http://freeimage.sourceforge.net/fnet/html/A633E9A9.htm but obviously GDAL takes just the first string.
The second error means that the writer ...
2
You can suppress the warnings (as long as you're sure there's no real issue with your data) with gdal.PushErrorHandler('CPLQuietErrorHandler'). If you do this any errors will also not get printed, so make sure you tell GDAL to raise a Python exception when an error occurs with gdal.UseExceptions().
E.g.
# Stop GDAL printing both warnings and errors to ...
1
You need to rescale the transform as well as the data, i.e.
from rasterio import Affine
from rasterio.enums import Resampling
scale = 10 # Reduce/upscale resolution scale factor
t = src.transform
# rescale the metadata
transform = Affine(t.a * scale, t.b, t.c, t.d, t.e * scale, t.f)
height = int(src.height / scale)
width = int(src.width / scale)
band = ...
0
I like to use gdal for these kind of operations. If you have access to the tools via the command line, it is very easy to resample your imagery at different resolutions like this example that would change "file1.tif" to a new file called "file1_0.5m.tif" which has been resampled at a resolution of 0.5m x 0.5m.
You can use gdalinfo to compare the resolution ...
0
If you're interested in using GRASS for this, here are the steps. To smooth the original, r.neighbors creates a new raster where each pixel gets the average value of surrounding pixels within a window. THen you can use r.recode to get a new raster with discrete categorical values replacing the continous population density. And finally, r.to.vect will convert ...
1
Note where the error is coming from and step through your code. Its not even getting to the Savitsky-Golay line.
> v = as.vector(r)
> z = na.spline(v)
> s1.ts2 = ts(z, start=1, end=nlayers(z), frequency=12)
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘nlayers’ for signature ‘"numeric"’
The error ...
1
Turn on snapping to vertices in the Snapping Toolbar.
Turn on editing for the orange layer.
Digitize a new polygon to fill in each grid square.
The snapping setting will cause the split tool to snap to the corners of adjacent grid squares, which will make the new polygon edges align exactly with the rest of the grid.
Tips:
Next time, don't delete the ...
0
As usual, it was a simple thing that was missed. I had neglected to set the output coordinate system and apply a transformation. I added these lines at the beginning:
env.outputCoordinateSystem = arcpy.SpatialReference("NAD 1983 UTM Zone 13N")
env.geographicTransformations = "WGS_1984_(ITRF00)_To_NAD_1983"
The rasters then rendered properly.
0
you can use this code....
// for graph/chart that contain mean NDVI of every year
var plotNDVI = ui.Chart.image.seriesByRegion(addnullimages, studyarea,ee.Reducer.mean(),
'nd',500,'system:index')
.setChartType('LineChart').setOptions({
title: 'NDVI short-term time series',
hAxis: {title: 'Date'},
...
0
// for graph/chart that contain mean NDVI of every year
var plotNDVI = ui.Chart.image.seriesByRegion(addnullimages, studyarea,ee.Reducer.mean(),
'nd',500,'system:index')
.setChartType('LineChart').setOptions({
title: 'NDVI short-term time series',
hAxis: {title: 'Date'},
vAxis: {title: 'NDVI'...
0
It seems to me that Matlab GeoTIFF reader can't deal with nodata value or with alpha channel https://se.mathworks.com/help/map/ref/geotiffread.html. Transparency is available only for PNG, CUR and ICO files https://se.mathworks.com/help/matlab/ref/imread.html.
Therefore the only way to remove the nodata areas from Matlab on QGIS side is probably to save the ...
0
Try something like:
readandproject = function(filenames){
lapply(filenames, function(filename){
r = raster(filename)
crs(r) = "+init=epsg:4326" # fix broken CRS
r = projectRaster(r, crs="+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
return(r)
}
)
}
Then prasters = readandproject(c("file1.nc","file2.nc")) ...
0
I wondered if your code was only applying itself to the first instance, as the indentation isn't correct for a for loop.
for index in index_list:
# set extent
env.extent = index
env.snapRaster = index
env.cellSize = index
This might make the code apply to all the indexes in your index_list.
1
Calculate Statistics on the Euc Dist raster (to populate maximum),
In the Raster Calculator:
"eucdistraster".maximum - "eucdistraster"
0
Multiple solutions possible to get your a very close result. If you can create 0.7 contour of raster, points on it are your best candidates.
Alternatively express your area in raster pixels and run focal statistics using relevant size window, e.g. 50*50 (not necessarily a square or rectangle) if your area is equivalent of 2500 cells. Use raster calculator:
...
0
Using these packages (MASS gets you the negative binomial distribution functions):
> library(MASS)
> library(raster)
Let's create a raster of coefficients that you have computed in some way from your landscape. Here its a 3x4 grid of numbers from 1 to 12:
> C = raster(matrix(1:12,3,4))
Then we can create negative binomial samples with a mu value ...
1
If all you need to do is visually represent "number of times an area has been burnt", there's a very simple solution with the polygon layer(s)
Simply make the polygon layers all the same color, and set their opacity to 30-50%. The more overlapping polygons, the darker they appear.
Here's what a layer looks like at 30% opacity.
1
Set the 'burnt' area of each of your 16 rasters equal to 1 and use r.null on each of your rasters to convert 'nodata' values to 0.
To use r.null, set "The value to replace the null value by [optional]" to 0 and set the "GRASS GIS 7 region extent" to anything larger than the total extent of all 16 rasters.
Then use raster calculator to add all 16 layers. ...
2
QGIS is out, it doesn't support raster attribute tables at all.
I don't know about GRASS.
GDAL would require you to (manually or in python code) handcraft VRTs with lookup tables from the attribute table fields for each "band", then composite the VRTs.
In ArcGIS, you can only use Lookup and Composite Bands. However, look at the versions of those functions ...
2
Loop over rows of the data frame and extract the i-th element from the vector:
for(i in 1:nrow(nsol)){
d = nsol$Date[i]
hillshade(..., nsol$Aspect[i], nsol$Azimuth[i],...etc...)
}
1
You normally would ST_ReClass() the useless values of your raster to nodata (NULL). Then you can convert the withvalue pixel zones to polygons with ST_DumpAsPolygons() or you could intersect them with other polygons with ST_Clip() or ST_Intersection().
2
have you tried to polygonize you raster? then you can clip the original raster with the right polygon and you should be good!
1
You would need to import the ASCII file as points, then interpolate a DEM. A quick google search will return many explanations. For example, See:
this answer to Creating DEM from Point Data using QGIS?
Then create a slope raster from the DEM. Again, google it. i.e.: Creating slope map from DEM using QGIS?
0
The image resolution of your exported geotiff won't affect the scale of the map you produce with it. When you produce a map in the print composer, QGIS will automatically zoom in or out on the image to the same scale as the map (that's the point of making it a georeferenced image rather than just a picture).
The resolution you choose just has to do with ...
0
I would use a vector mask or shapefile if the raster image is geospatially defined as part of a map especially if the raster is georeferenced in world coordinates (lat,lon, WGS1984, etc so forth) simply because it's easier for me to take a polygon of let's say Denver, Colorado where I live and extract satellite data using that polygon either via clip or mask....
0
In Qgis, try using using Raster -> Conversion -> Polygonize (Raster to Vector)
0
I had the same problem and solved it by crafting a syntactically valid JSON response by using FreeMarker. Of course, the GetFeatureInfo output type still is HTML (and so is the response header of the request) but luckily my client was not checking the header but just parsed the JSON-styled response as pure JSON:
header.ftl:
{"type":"FeatureCollection","...
0
Your coords are:
[,1] [,2]
[1,] 14.650 40.350
[2,] 13.600 39.350
[3,] 0.033 -78.717
[4,] 37.017 -2.383
[5,] 46.837 -113.966
but coordinate matrices are always X then Y, ie Longitude then Latitude. Your NAs are coming from column 2 being more than 90 or less than -90, which are not valid latitudes.
extract(r, coords[,2:...
2
A couple of things here:
Firstly don't set gdata to None until you are COMPLETELY done with that raster, data will return a broken dataset if you do. Please see for a more detailed explanation: https://trac.osgeo.org/gdal/wiki/PythonGotchas
Next thing uninstall gdal using pip then reinstall it using the .whl file from here:
https://www.lfd.uci.edu/~...
0
So I found a great blog post here that answers my question.
# Load the contour data and subset by layer required
readOGR("myfolder/contours.shp") %>%
subset(FEAT_TYPE =="ContourLine") -> c
# Create a target raster for the DEM
dem_bbox <- extent(c)
dem_raster <-raster(dem_bbox)
projection(dem_raster) <- CRS(projection(c))
# Set ...
0
From documentation:
writeRaster
Write Raster Data To A File
Write an entire Raster* object to a file, using one of the many
supported formats.
When writing a file to disk, the file format is determined by the
'format=' argument if supplied, or else by the file extension (if the
extension is known). If other cases the default format is ...
0
You should try two simple gdal tools: gdaltindex and gdalbuildvrt.
To create the mosaic, since your rasters are georreferenced, you just need:
gdalbuildvrt mosaic.vrt *.tif
QGIS will open this mosaic.vrt raster.
If you want to create the grid around your rasters, you can use gdaltindex.
gdaltindex mosaic.shp *.tif
Open both the mosaic.vrt and mosaic....
0
i was able to accomplish this task by subsetting the data from the data frame which matched my classes. Here is the updates code section:
changeDet1 <- calc(stack(lc1,lc2), fun = change)
codes_ <- data.frame(ID = grid_$code,value = paste0('from ',grid_[,1],' to ',grid_[,2]))
logical_test <- which(grid_$change == T) # remove no change classes
codes_ ...
1
You are over complicating this a bit. Certain types of data are suitable for interpolation whereas other data is more well-suited for simple binning. In the lidar realm, elevation and intensity are commonly interpolated but, other lidar attributes (class, return number, time-stamp, strip-id) are simply binned to a pre-existing grid. From an analytic ...
2
Use a value other than 0.0 for the burn-in value, such as 1.0.
Use the extent of the land cover raster or the 10km grid
Change output raster size units to map instead of pixel,
Specify 10000 (your map units are metres, i.e 10km) for the resolution.
2
If the rasters have correct location but one of them is rotated, then you can georeference incorrect raster based on the correct (georeferenced) one. You can either use Raster -> Georeferencer tool to do the job, or use Freehand raster georeferencer plugin which has the function to rotate the raster data.
You can download Freehand raster georeferencer ...
1
You should be able to use the Conditional tool, either standalone or in Raster Calculator. Set your conditional statement so that wherever Mask=Water will be 1, and everywhere else will be 0. Something like this:
Con(Water=Mask, 1, 0)
This will produce a binary raster marking the waterline. If you need to know the height of the waterline you could then use ...
0
Use the GDAL WMS/TMS/ArcGis driver and gdal_translate.
Create an XML file specifying the server parameters, CRS and extent, for example:
<GDAL_WMS>
<Service name="TMS">
<ServerUrl>http://tile.openstreetmap.org/${z}/${x}/${y}.png</ServerUrl>
</Service>
<DataWindow>
<UpperLeftX>-20037508....
1
Yes, gdal from out the command line should do the trick, and will be one of the fastest options! Now the question is in which format you want it. As it is quite big I would suggest using the format binary grid. To make such a translation, you need to make use of a commandline with gdal installed in it. Supposing you use windows: Press the windows key and ...
2
If you raster's projection unit is meter, then the output unit of the proximity raster will be meters if you specify georeferenced units. In case your raster's projection unit is not meter (e.g. is degree because it is on WGS84), you'll want to reproject your raster to another coordinate system (such as UTM, for example).
0
Look at the object
r1
or check if it has an associated filename:
filename(r1)
The source points to a filepath on someone's computer, and what happened was they created a raster object linked to a file i.e.
s <- stack("MOD13Q1_EVI_2008_097.tif")
and then saved it out
save(s, "myfile.RData")
All that did was save a summary of the object, like the ...
0
This can be done using the
MMQGIS plugin
Once added it will show up in the tool bar at the top of the screen. Select/highland the 'line.data' layer you want to convert to 'points.data'. click on MMQGIS and a drop down will appear. Then click on Modify and then Convert Geometry Type
Then a box appears as seen below -->
Change the 'New Geometry Type' ...
Top 50 recent answers are included
