17

With PyQGIS ras = QgsRasterLayer("raster.tif") pixelSizeX= ras.rasterUnitsPerPixelX() pixelSizeY = ras.rasterUnitsPerPixelY() print pixelSizeX 2.11668210081 print pixelSizeY 2.11685012701 With GDAL from osgeo import gdal raster = gdal.Open('raster.tif') gt =raster.GetGeoTransform() print gt (258012.37107330866, 2.11668210080698, 0.0, 163176.6385398821, 0....


16

Using rasterio you could do import rasterio file_list = ['file1.tif', 'file2.tif', 'file3.tif'] # Read metadata of first file with rasterio.open(file_list[0]) as src0: meta = src0.meta # Update meta to reflect the number of layers meta.update(count = len(file_list)) # Read each layer and write it to stack with rasterio.open('stack.tif', 'w', **meta) ...


12

The Python API method that supports the rio-sample command is documented here: https://rasterio.readthedocs.io/en/latest/api/rasterio._io.html#rasterio._io.DatasetReaderBase.sample src.sample() takes an iterator over x, y tuples, so do: for val in src.sample([(x, y)]): print(val)


10

Reading Rasters by block can be done in rasterio and I'd argue it's easier than in GDAL. There is even a tutorial on windowed read/write over at GitHub. Let's take a look at the read functions arguments, which allows you to set a window to read data from: def read(self, indexes=None, out=None, window=None, masked=False, boundless=False): """Read ...


10

If using GDAL 2.1+ it's as simple as gdal.BuildVRT then gdal.Translate: from osgeo import gdal outvrt = '/vsimem/stacked.vrt' #/vsimem is special in-memory virtual "directory" outtif = '/tmp/stacked.tif' tifs = ['a.tif', 'b.tif', 'c.tif', 'd.tif'] #or for all tifs in a dir #import glob #tifs = glob.glob('dir/*.tif') outds = gdal.BuildVRT(outvrt, tifs, ...


9

I added that recipe to the rasterio documentation. Since it was such a simple shape, in this case I just unzipped the coords in the single record contained by the shapefile. That is, x, y = zip(*features[0]['coordinates'][0]), and then just plot. More generally, I use PolygonPatch from descartes, and matplotlib.collections. import fiona import rasterio ...


9

I believe you are using the wrong bucket name. It should be s3://sentinel-s2-l1c, not s3://sentinel-pds Try: url = 's3://sentinel-s2-l1c/tiles/10/S/DG/2015/12/7/0/B01.jp2' See more info here: http://sentinel-pds.s3-website.eu-central-1.amazonaws.com/ and https://aws.amazon.com/public-datasets/sentinel-2/


9

Issues: The Anaconda default gdal may be built without BigTIFF support. If I create a non conda-forge env, i.e conda create -n testgdal gdal I can reproduce the md['DMD_CREATIONOPTIONLIST'].find('BigTIFF') == -1 no BigTIFF issue. There seems to be an incompatibility between the latest version of conda and vs2015_runtime and conda-forge. I updated my conda ...


8

If your GeoJSON geometry is unique (not a "type": "FeatureCollection") as in geoms = {'type': 'Polygon', 'coordinates': [[(250542.40328375285, 141691.07089614146), (250641.30366207045, 141400.7504307576), (250421.17056194422, 141512.41214821293), (250542.40328375285, 141691.07089614146)]]} and you try with rasterio.open("a.tif") as src: out_image, ...


8

This ended up being more straightforward than I thought, with all of the capabilities lying in the rasterio.open function. Here is an example using a proj4 string instead of wkt. import rasterio from rasterio.transform import from_origin arr = np.random.randint(5, size=(100,100)).astype(np.float) transform = from_origin(472137, 5015782, 0.5, 0.5) ...


7

import rasterio import numpy as np from rasterio.merge import merge from rasterio.plot import show src1 = rasterio.open('/path/to/your/raster1') src2 = rasterio.open('/path/to/your/raster2') # Taking a peek to make sure these are the rasters you want... show(src1) show(src2) srcs_to_mosaic = [src1, src2] # The merge function returns a single array and ...


7

The meta property contains the basic raster metadata. The profile property is a super set of meta and dataset creation options (i.e inc. tiling, block size, compression etc...). You would use src.profile when you want to ensure you create an exact (empty) clone of an existing dataset. And it is especially useful when you are reading in block windows from a ...


7

A rasterio way of doing this is pretty simple. Note this requires your raster be in the same projection as your coordinates. You can of course project your coordinates on the fly, but that's another question... import rasterio as rio infile = r"C:\Temp\test.tif" outfile = r'C:\Temp\test_{}.tif' coordinates = ( (130.5, -25.5) , # lon, lat of ~centre ...


6

First I would use bands 4(red) and 5(nir) for Landsat 8 according to the description of the OLI instrument, and 3(red) and 4(NIR) for the Landsat TM and ETM. Second, you define an output in dtype=rasterio.uint16, but NDVI should be a float (between -1 and 1). You should either initialize your raster as dtype=rasterio.float32 , or multiply your values by ...


6

Looking over your code, I don't see how it will be faster than pure Python. It's exclusively calling Python methods (rasterio methods, np.dstack is Python) and those aren't executed any faster just because the function is compiled with Cython. The key to speeding things up is in here: # Iterate over rows and columns in this block for x in xrange(...


6

That's because these bands come as unsigned integer 16 so the numpy division returns only positive integers. You can replace ndvi = (NIR-RED)/(NIR+RED) by ndvi = (NIR.astype(float) - RED.astype(float)) / (NIR+RED)


6

I'd recommend using gdal_contour . The results are likely to be much better than any attempt to re-implement it :-) Having said that, there should be a way to do this in rasterio. I've not tried this, so it may behave differently to how I'd expect it to work in QGIS. Step 1. Quantize the raster into contour bands using rasterio. Use rasterio as a raster ...


6

Are you looking to take the average for each grid cell in the stack, or the overall average? If it is the former, you could use the AverageOverlay tool in the WhiteboxTools library. This can be scripted in Python as follows: from whitebox_tools import WhiteboxTools wbt = WhiteboxTools() wbt.work_dir = "/path/to/data/" wbt.average_overlay(inputs='file1.tif;...


5

Try adding this to your script: from rasterio import logging log = logging.getLogger() log.setLevel(logging.ERROR) You can substitute ERROR with other levels: DEBUG, INFO, WARN, FATAL.


5

Answered at https://github.com/mapbox/rasterio/issues/710. Example reading a 30x30 window into a 3x3 array where overviews (if available) would kick in. arr = np.empty(shape=(3, 3)).astype(src.profile['dtype']) arr = src.read(1, out=arr, window=((0, 30), (0, 30))) arr array([[9195, 9116, 9134], [9158, 9144, 9085], [9010, 8935, ...


5

The source coordinate reference system (CRS) is WGS84 / UTM zone 18 North. The points are south of the equator. I'm assuming the destination CRS is WGS84, EPSG:4326. One of the programs is "reboxing" the raster while the other one isn't. When I unproject all 4 corners of the raster: 496950.0 -861870.0 496950.0 -1046760.0 599040.0 -1046760.0 599040.0 ...


5

The issue is resolved. The issue was I misread the documentation. On a second read, the rasterio.mask documentation clearly states that polygons should be a list of GeoJSON-like dicts. I found the following snippet of code to generate those lists from this answer: geoms = [feature["geometry"] for feature in shapefile] Here is the the full code that is ...


5

The previous answers are misleading or wrong. To modify the nodata value of a GeoTIFF with Rasterio, do this with rasterio.open(tiffname, 'r+') as dataset: dataset.nodata = -32767 The project has tests of this usage that you may see also: https://github.com/mapbox/rasterio/blob/master/tests/test_update.py#L59-L64. Note that you may have to close and ...


5

For me it works if I pass the matplotlib ax object explicitly to rasterio.plot.show: fig, ax = plt.subplots(figsize=(15, 15)) rasterio.plot.show(raster, ax=ax) countries.plot(ax=ax, facecolor='none', edgecolor='r'); Full example with raster data from naturalearth (link) (both datasets have the same CRS, so I don't do any reprojection): import geopandas ...


5

In rasterio <= 1.0.7 this is not possible. You have to write and close, then read. If you want to write again, just overwrite the memfile object. For example: from rasterio.io import MemoryFile import rasterio import numpy as np def create_memory_file(data, west_bound, north_bound, cellsize, driver='GTIFF'): #data is a numpy array if data.ndim ...


5

You could use shapely's affine_transform: import geopandas as gpd import rasterio as rio from shapely.affinity import affine_transform infile = "/tmp/test.geojson" outfile = "/tmp/out.geojson" raster = "/tmp/test.tif" df = gpd.read_file(infile) with rio.open(raster) as ds: t = ds.transform df.geometry = df.geometry.apply(lambda geom: ...


4

There is no such option as far as I know, therefore you need to use some tricks: option one: rasterize your polygon convert the polygons to line rasterize the line add the rasterized line to the existing rasterized polygon (1+1=2) option two apply a negative buffer of 1/2 pixel on each polygon rasterize the result --> countours and backgrounds will be ...


4

rasterio lets you define windows for read/write functions. Assuming your file has dimensions of 1000 by 1000 pixels you could call rasterio to only read the upper left quarter of the image: import rasterio with rasterio.open('your/data/2000/doy1.tif') as src: w = src.read(1, window=((0, 250), (0, 250))) Where window is given as ((row_start, row_stop), ...


4

Thanks to @LoïcDutrieux for pointing out in the comments that you can use pip to install the plugin. A quick search indicates the plugin does exist in PyPI: $ pip search rio-color rio-color (0.4.0) - Color correction plugin for rasterio To install, use: pip install rio-color If you are running Windows, pip is installed by default with ...


4

Rasterio can read and write GCPs and warp with them since version 1.0a3. src.crs returns nothing to make it clear that there is no coordinate reference system associated with the file's affine transformation matrix. If you want to see the file's ground control points and their CRS, do this (one of the project's test files shown for example). >>> ...


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