pyqgisqgisraster
After much searching, I've still not been able to figure out how to perform what seems like a simple function: I'd like to get the pixel size/cell size of a raster layer using PyQgis. Is there any native way to do this through QgsRasterLayer or a helper class? Or will I need to read the file in using GDAL functions, e.g.:
geotransform = dataset.GetGeoTransform()
if not geotransform is None:
print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'
Thanks for your help.
Alex
Following @Dango's idea I created and tested (on small rasters with the same extent and cell size) the following code:
import arcpy, numpy
inRaster = r"C:\tmp\RastersArray.gdb\InRaster"
inRaster2 = r"C:\tmp\RastersArray.gdb\InRaster2"
##Get properties of the input raster
inRasterDesc = arcpy.Describe(inRaster)
#coordinates of the lower left corner
rasXmin = inRasterDesc.Extent.Xmin
rasYmin = inRasterDesc.Extent.Ymin
# Cell size, raster size
rasMeanCellHeight = inRasterDesc.MeanCellHeight
rasMeanCellWidth = inRasterDesc.MeanCellWidth
rasHeight = inRasterDesc.Height
rasWidth = inRasterDesc.Width
##Calculate coordinates basing on raster properties
#create numpy array of coordinates of cell centroids
def rasCentrX(rasHeight, rasWidth):
coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)
return coordX
inRasterCoordX = numpy.fromfunction(rasCentrX, (rasHeight,rasWidth)) #numpy array of X coord
def rasCentrY(rasHeight, rasWidth):
coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)
return coordY
inRasterCoordY = numpy.fromfunction(rasCentrY, (rasHeight,rasWidth)) #numpy array of Y coord
#combine arrays of coordinates (although array for Y is before X, dstack produces [X, Y] pairs)
inRasterCoordinates = numpy.dstack((inRasterCoordY,inRasterCoordX))
##Raster conversion to NumPy Array
#create NumPy array from input rasters
inRasterArrayTopLeft = arcpy.RasterToNumPyArray(inRaster)
inRasterArrayTopLeft2 = arcpy.RasterToNumPyArray(inRaster2)
#flip array upside down - then lower left corner cells has the same index as cells in coordinates array
inRasterArray = numpy.flipud(inRasterArrayTopLeft)
inRasterArray2 = numpy.flipud(inRasterArrayTopLeft2)
# combine coordinates and value
inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))
#add values from second raster
rasterValuesArray = numpy.dstack((inRasterFullArray, inRasterArray2.T))
Based on @hmfly code, you can have access to desired values:
(height, width, dim )=rasterValuesArray.shape
for row in range(0,height):
for col in range(0,width):
#now you have access to single array of values for one cell location
Unfortunately there's one 'but' - the code is right for NumPy arrays which can be handled by system memory. For my system (8GB), the largest array was about 9000,9000.
As my experience doesn't let me provide more help, you can consider some suggestions about dealing wiht large arrays: https://stackoverflow.com/questions/1053928/python-numpy-very-large-matrices
arcpy.RasterToNumPyArray method allows to specify the subset of raster converted to NumPy array (ArcGIS10 help page) what can be useful when chunking large dataset into submatrices.
I'm not sure where you get the above-mentioned 'attribute' variable from, so I just assumed rownames(sp_df@data) (ranging from 0 to 193) as the desired output. Anyway, here is my approach that is probably not the fastest (or the most convenient in general), but it works.
Basically, the script starts to loop through all the cells of the generated raster template (100 m spatial resolution) and crops sp_df by the extent of the current pixel (which is extracted by setting all but the current cell to NA and then converting to polygons). After that, it checks whether no, one or more than one polygon IDs are to be found within the current tile, calculates the areas of the referring IDs and returns the ID covering the largest area. Subsequently, the thus extracted values are inserted into the referring raster cells.
# Required libraries
library(raster)
library(rgdal)
library(rgeos)
# Import shapefile
sp_df <- readOGR("data/", "stok_subs")
# Raster template (100 m spatial resolution)
r <- raster(extent(sp_df))
projection(r) <- proj4string(sp_df)
res(r) <- 100
# Per pixel, identify ID covering largest area
r_val <- sapply(1:ncell(r), function(i) {
r_dupl <- r
r_dupl[i] <- 1
p <- rasterToPolygons(r_dupl) # Current cell -> polygon
sp_df_crp <- crop(sp_df, p) # Crop initial polygons by current cell extent
# Case 1: no polygon intersecting current cell
if (is.null(sp_df_crp)) {
return(NA)
# Case 2: one polygon intersecting current cell
} else if (nrow(sp_df_crp@data) < 2) {
return(rownames(sp_df_crp@data))
# Case 3: multiple polygons intersecting current cell
} else {
areas <- gArea(sp_df_crp, byid = TRUE)
index <- which.max(areas)
return(rownames(sp_df_crp@data)[index])
}
})
# Write ID values covering the largest area per pixel into raster template
r[] <- as.numeric(r_val)
plot(r)
plot(sp_df, border = "grey45", add = TRUE)
Best Answer
You tagged your question with pyqgis. Here is a small solution using the Pyqgis modules