20

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) ...


18

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....


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)


12

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, ...


11

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 ...


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

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 ...


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) ...


8

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 ...


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 ...


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

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, ...


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

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 ...


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;...


6

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

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

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

You can read specific bands in a single call using rasterio by passing a list/tuple of band numbers (Following the GDAL convention, bands are indexed from 1): import rasterio rasterio.__version__ '1.0a8' dataset = rasterio.open('multiband.tif') dataset.count 4 dataset.read((1,2)) #read 1st two bands into an array. array([[[ 85, 98, 75, ..., 53, 55, ...


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

You can directly write the result blocks to the big result tiff file via windowed writing. import rasterio src_RP1 = rasterio.open(r'RP1_filename.tif') src_RP2 = rasterio.open(r'RP2_filename.tif') out_profile = src_RP1.profile.copy() dst = rasterio.open(r'result.tif', 'w', **out_profile) for block_index, window in src_RP1.block_windows(1): RP1_block =...


5

The Rasterio Plotting documentation describes how to visualize multiband imagery. For example, using 4-band NAIP imagery: import rasterio from rasterio.plot import show src = rasterio.open("path/to/your/image/m_3511642_sw_11_1_20140704.tif") show(src) To visualize specific band combination use the following approach (source). In this case, I am creating a ...


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: ...


5

Yes, you can write a one bit raster with rasterio*. You need to: write to a format that supports a 1bit dataset, such as GeoTIFF; ensure your numpy array is np.uint8/ubyte so rasterio doesnt raise the TypeError: invalid dtype: 'bool' exception; and pass the NBITS=1 creation option to tell the underlying GDAL GeoTIFF driver to create a one bit file. import ...


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