Hot answers tagged

64

You can use the gdal.Dataset or gdal.Band ReadRaster method. See the GDAL and OGR API tutorials and the example below. ReadRaster does not use/require numpy, the return value is raw binary data and needs to be unpacked using the standard python struct module. An example: from osgeo import gdal,ogr import struct src_filename = '/tmp/test.tif' shp_filename ...


17

Assuming you are talking about a map of the world in Mercator Projection, Something like this: Then, You should know that not all pixels will represent the same area. There is a large amount of distortion as you move towards the poles, the area represented by each pixel will decrease. Tissot's indicatrices show a good visual representation of the ...


8

Yes. A raster pixel can hold one or more values in at least two ways that I know of: You can store multiple bands in a raster, so each pixel at a point can have as many values as you have bands; You can encode multiple values in the pixel of a single band using bit flags. Bitwise operators can be used to extract these values. Basically, for a raster band ...


6

Here's a quick script I put together to correct the cell size of a folder of GeoTiff rasters: import os, sys from osgeo import gdal from osgeo import gdalconst WorkingFolder = sys.argv[1] # first command line argument # change to a hard path like r'c:\your\path' or 'c:\\your\\path' # without the r if it suits your purposes.. for f in os.listdir(...


6

Try out my code (as script) at the Python Console of QGIS: from os import path import struct from osgeo import gdal def countRasterValue(val): layer = iface.activeLayer() provider = layer.dataProvider() fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'} my_path = provider....


6

The easiest way is to combine numpy's NaN functions (in this case nanmean) and ndimage.generic_filter, like so: import numpy as np from scipy import ndimage m = np.array([[13, 21, 13, 8], [ 5, 10, 22, 14], [21, 33, 9, 0], [ 0, 0, 0, 0]], dtype=np.float) result = ndimage.generic_filter(a, np.nanmean, size=3, ...


6

Per https://lpdaac.usgs.gov/dataset_discovery/modis, the viewing swath width of MODIS is 2,330 km, thus a large portion of the image is off-nadir in some way. https://modis.gsfc.nasa.gov/about/specifications.php The following forum post gives an explanation of how to calculate pixel size based on viewing position. (Note: still an estimate due to factors ...


5

Yes, it is possible to get pixel values without gdal; only PyQGIS. Easier way is with a QgsRasterBlock object. This is the code: layer = iface.activeLayer() #in my case, a 20x20 raster provider = layer.dataProvider() extent = provider.extent() rows = layer.height() cols = layer.width() block = provider.block(1, extent, cols, rows) for i in range(rows): ...


5

This a known anomaly - see anomaly #29 in the Sentinel 2 data quality report. https://sentinels.copernicus.eu/web/sentinel/missions/sentinel-2/data-quality-report


5

NDBI (like NDVI) is an index that can be useful to highlight built up areas, but it is not possible to use a single threshold to determine with 100% certainty if the observed land cover is built up. The use of indices should therefore be done carefully: while they can be very powerfull, in some cases it doesn't work and you need to use all spectral bands. In ...


4

Ground sampling distance is the distance between each measurement at nadir. Pixel size refers to what is delivered when data is purchased. Some (commercial) data providers use a smaller pixel size than the ground sampling distance. This means that the data is interpolated from the measured image grid into the delivered image grid. In the case of RapidEye, ...


4

I would generate a new raster by right-clicking on the image in the Table of Contents and Data > Export Data. Then change the cell size to something much less detailed, e.g. 10m or 100m etc. to suit. The new output raster will be much more pixelated. (Example below at 25m cell size)


4

According to this question, How to calculate distance in meter between QgsPoint objects in PyQgis?, it's easy. You just need to return a QgsPoint with your function pixel2coord() : def pixel2coord(x, y): xp = (pixelWidth * x) + originX + (pixelWidth/2) yp = (pixelHeight * y) + originY + (pixelHeight /2) pnt = QgsPoint(xp, yp) return pnt ...


4

I assume that by "lapse rate" you mean the regression slope. You are essentially describing a focal regression and since there is no multivariate version of raster::focal, there is no canned way to do what you are after. This could be accomplished through a series of focal matrix algebra functions. However, you can also leverage the read/write functions that ...


4

You are specifying the raster width/eight in terms of pixel, not the pixel width/height in terms of meters, so if the raster width is 10, it means the pixel size will be vector extent width divided by 10. As this number of pixels increases, the size of each pixel decreases. Assuming your data projection is expressed in meters, you can choose the option ...


4

In the ArcGIS raster calculator you can use the Con, IsNull and FocalStatistics functions: Con(IsNull("YourRaster") & (FocalStatistics("YourRaster", NbrRectangle(3, 3, "CELL"), "VARIETY", "") < 5), FocalStatistics("YourRaster", NbrRectangle(3, 3, "CELL"), "MAJORITY", ""), "YourRaster") Reformatted for reading: Con( IsNull("YourRaster") & (...


4

You can use ST_PixelAsPoints to get a geometry for each point and then to a "Cartesian join", ie, a full join between the points and the pixels, eg, WITH pixels (x, y, geom) AS ( SELECT x, y, geom FROM (SELECT (ST_PixelAsPoints(rast, 1)).* FROM rast_table WHERE rid = 1 ) SELECT gid, x, y, ST_Distance(points....


4

The unit is degrees longitude and latitude, which is a result of your georeference system being epsg:4326. Pixel size and pixel spacing are used interchangeably, as the concept of overlapping pixels is generally ignored. The best way to convert the numbers to meters would be to reproject your data into a more suitable coordinate system, such as the ...


3

Upsample the image to the suspected resolution and Downsample it to the original. If the resulting and original are the same (or close enough), it is very likely that you have a lower resolution image.


3

Once you did your analysis you can Copy your output with arcpy.CopyRaster_management. Where you can specify bit depth through the pixel_type arg. For example - arcpy.CopyRaster_management("inrast","outrast","","","-9999", "","","8_BIT") You can then delete the earlier output to cleanup.


3

By using 'identify' method of QgsRasterDataProvider (PyQGIS), it is not very difficult to do that. I tested my approach with point and raster layers of next image. For simplicity, I only have considered the first feature of point layer; but the procedure can be easily extended for all features. The code is: mapcanvas = iface.mapCanvas() layers = ...


3

You should check your Environment Settings and adjust the Cell Size to an appropriate cell size: ArcGIS Cell Size


3

Simpler and faster with gdalnumeric (gdal+numpy): import gdalnumeric raster_file = gdalnumeric.LoadFile("/path/to/file.tif") pixel_count = (raster_file == 1).sum() # for pixel value = 1 print(pixel_count) Inside Qgis (for the active layer selected): import gdalnumeric layer = iface.activeLayer() provider = layer.dataProvider() raster_file = gdalnumeric....


3

As discussing the possible reasons with @ahmadhanb, I have realised the only issue was my ignorance. If the help page linked to question is read properly, one can see the only bit depth option for Float data is Floating-point 32 bit, therefore ArcGIS promotes the pixel depth to 32-bit if any float raster is selected as input/parameter. Possible resolution ...


3

There were some missing pieces of code. I used the solution by @YoLecomte for obtaining a working solution: from qgis.core import * from osgeo import gdal, ogr, osr import numpy def pixel2coord(x, y): xp = (pixelWidth * x) + originX + (pixelWidth/2) yp = (pixelHeight * y) + originY + (pixelHeight /2) return QgsPoint(xp, yp) pntLayer = ...


3

I used your code (slightly modified): from osgeo import gdal import numpy def pixel2coord(x, y): xp = (pixelWidth * x) + originX + (pixelWidth/2) yp = (pixelHeight * y) + originY + (pixelHeight /2) return QgsPoint(xp, yp) #pntLayer = QgsVectorLayer("/Data/points.shp","pointLayer",'ogr') pntLayer = iface.activeLayer() feats = [ feat for feat ...


3

Based on @Michael Stimson in above and QGIS document (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/raster.html#query-values). In below example it is possible to get pixel values without gdal. from qgis.core import * def pixel2coord(x, y): xp = (pixelWidth * x) + originX + (pixelWidth/2) yp = (pixelHeight * y) + originY + (...


3

The raster package has a pixel area calculation function for unprojected rasters. Given a RasterLayer named r library(raster) area(r) * r


3

Seninel 2's Technical documentation Cloud Masks It appears to be calculated based on reflectance, however there is a known issue as mentioned Two products have been found affected by this anomaly. The products have very small data coverage and are completely cloudy. The cloud mask is accurate but the cloud coverage metadata is reported as zero. The ...


3

You can try Zonal Statistics as Table: Summarizes the values of a raster within the zones of another dataset and reports the results to a table. I belive one of the outputs will be pixel Count (per zone).


Only top voted, non community-wiki answers of a minimum length are eligible