Hi,

We are currently using rasterio version 1.1.8 and it was installed following this procedure : https://rasterio.readthedocs.io/en/latest/installation.html#installing-with-anaconda

Thank you for the hint, I wasn´t aware of this ArcGIS behaviour concerning binary raster files. A first test was successful!

Hello everybody,

I´m quite new to rasterio and Python, but my task is, to produce binary raster files based on the information of sentinel2 satellite images (GTiff). Therefore I use Spyder IDE and the numpy.where() method to calculate new arrays and save them as GTiffs.

Although the output arrays contain 0 and 1 values, the output tif only shows black (in the shape of the original gtiff), if opened with standard image viewer or ArcGIS Desktop. Only if opened via Python script in Spyder, I am able to display the output tifs as expected (using rasterio.plot.show() method).

If anybody has an idea, what the problem might be, I would appreciate any help.

Thanks!

It may be possible that in-memory HDF4 is not supported by GDAL.

Hi Sean, 

Thanks for your quick answer and for noting this, i had indeed a typo. However, this did not solve my problem. Am I missing something else ? There is not much info in the documentation and I feel like I included everything needed.

Replacing your FTP code with direct access using the /vsicurl/ virtual filesystem:

This is a very elegant solution, thanks for pointing it out and thanks for the sample code! The only drawback I see is that it is not possible to access multiple products (files within the archive) while at the same time downloading the archive only once.

Thanks again,
F

Hi,

I noticed some odd behavior when using WarpedVRT in combination with a VRT with mixed Data Types. 

When trying to open the mixed data type raster using rasterio, it throws an error due to no support for mixed data types. However, when trying to open the same dataset as WarpedVRT, I do not get such error. Instead it returns an all-zero array.

I would either expect to get the same error message, receive an array with all data cast to the largest data type or receive an array with mixed data types.

import rasterio
import subprocess
from affine import Affine

from rasterio.vrt import WarpedVRT

# First raster saved as uint16
a = np.array([[0,1,2,3,4]]).astype("uint16")

with rasterio.open("uint16.tif", "w", width=5, height=1, count=1, dtype=rasterio.uint16, nodata=0) as dst:

    dst.write(a, 1)

# Second raster saved as int32
b = np.array([[1,2,3,4,0]]).astype("int32")

with rasterio.open("int32.tif", "w", width=5, height=1, count=1, dtype=rasterio.int32, nodata=0) as dst:

    dst.write(b, 1)

# create a VRT combining both datasets into a two band raster using gdalbuildvrt

subprocess.run(["gdalbuildvrt", "-separate", "mixed_datatype.tif", "uint16.tif", "int32.tif"], capture_output=True)

# Try to open the two band raster using rasterio

with rasterio.open("mixed_datatype.tif") as src:

    print(src.profile)

    print(src.read())  # This throws an error due to mixed data types

{'driver': 'VRT', 'dtype': 'uint16', 'nodata': 0.0, 'width': 5, 'height': 1, 'count': 2, 'crs': None, 'transform': Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0), 'tiled': False}

Traceback (most recent call last):

  File "<input>", line 3, in <module>

  File "rasterio/_io.pyx", line 278, in rasterio._io.DatasetReaderBase.read

ValueError: more than one 'dtype' found

# Try to open the two band raster using WarpedVRT

with rasterio.open("mixed_datatype.tif") as src:

    with WarpedVRT(src, crs="EPSG:4326", transform=transform, width=5, height=1) as vrt:

        print(vrt.profile)

        print(vrt.read())  # This does not throw an error but returns an array with only zeros

    

{'driver': 'VRT', 'dtype': 'int32', 'nodata': 0.0, 'width': 5, 'height': 1, 'count': 2, 'crs': CRS.from_epsg(4326), 'transform': Affine(1.0, 0.0, 0.0,

       0.0, -1.0, 0.0), 'tiled': False}

[[[0 0 0 0 0]]

 [[0 0 0 0 0]]]

Hi Sean,

Ok thanks for the clarifications. I think I found GDAL's logic around here: https://github.com/OSGeo/gdal/blob/d01976cd4fc11144662f44ec4e658ce1871f1713/gdal/apps/gdalwarp_lib.cpp#L3105

It would indeed be useful to have a function that returns the best transform and shape given any bounds, resolution and shape. That's what I tried to do here in a Python implementation: https://github.com/GlacioHack/GeoUtils/blob/12bd607fe8024720ff264545abf0cc5ba0156ea6/geoutils/georaster.py#L879 (until L 932), if that can be of any help. In the meantime, I think I'll just keep my current implementation, which seems to provide similar results to what gdalwarp would do (I'll check the code in detail later).

Amaury

Hello,
I'm sorry if this is question is an easy one, I'm a beginner with rasterio. I have a python script that reads in a dataset, calculates an overlay layer and then outputs the calculated layer. This part works as expected. I then need to take this output geotiff and read it in for display inside a web application. Inside my web application, I need to extract data such as the geotiff's coordinate system , resolution and origin. Currently, I have some other geotiffs that were not generated using rasterio  and I'm able to fetch the needed information by extracting the geoKeys with npm's geotiff package : https://www.npmjs.com/package/geotiff. However, when reading my geotifffs generated with rasterio, no geoKeys are found. Below you can find my code used to open the existing dataset and the open and create options I'm using. I'm wondering if there is something that I'm missing to be able to correctly generate these geoKeys. Any help will be greatly appreciated!

def layers(directory, gis_data, maptype, mapname):
    h = gis_data.height
    w = gis_data.width
    c = gis_data.crs
    t = gis_data.transform
    output_folder = directory + 'output'
    
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)
    
    mapname = rio.open(
        output_folder + '/' + mapname + '.tiff',
        'w',
        driver = 'GTiff',
        height = h,
        width = w,
        count = 1,
        dtype = np.uint8,
        crs = c,
        transform = t,
        profile='GeoTIFF',
        geotiff_keys_flavour='STANDARD',
        georef_sources='INTERNAL',
    )
    mapname.write(maptype, 1)
    mapname.close()