Hi,

My team and I have been facing some problems when reading netcdf files in
`LUSTRE` filesystem with `ncdump` tool regarding an HDF5 1.10.1 issue in
`ubuntu:bionic` (see for example: https://github.com/ALPSCore/ALPSCore/issues/410)


We made a workaround by installing `netcdf-bin` using repositories of
`ubuntu:xenial` in `bionic`. Although `ncdump` have successfully read the
netcdf file, we couldn't fix it (sort of...) for `gdalinfo` and `rasterio`
for example:

this works:

ncdump /LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc|head
-n 20

gdalinfo /LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc
|head -n 20

gdalinfo -sd 1 /LUSTRE/MADMEX/dir_test_mount/
madmex_003_37_-32_1996-01-01.nc |head -n 20

but this neither of next lines work:

gdalinfo
netcdf:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:
blue_mean

gdalinfo
NetCDF:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:
blue_mean

```
ERROR 4:
NetCDF:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:blue_mean:
No such file or directory
gdalinfo failed - unable to open
'NetCDF:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:
blue_mean'.
```

This ultimately leads to error with `rio insp` in both next lines :

rio insp
netcdf:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:
blue_mean

rio insp
NetCDF:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:
blue_mean

```
ERROR:root:Exception caught during processing
Traceback (most recent call last):
File "rasterio/_base.pyx", line 214, in
rasterio._base.DatasetBase.__init__
File "rasterio/_shim.pyx", line 64, in rasterio._shim.open_dataset
File "rasterio/_err.pyx", line 205, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_OpenFailedError:
netcdf:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:blue_mean:
No such file or directory

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/usr/local/lib/python3.6/dist-packages/rasterio/rio/insp.py", line
77, in insp
with rasterio.open(input, mode) as src:
File "/usr/local/lib/python3.6/dist-packages/rasterio/env.py", line 423,
in wrapper
return f(*args, **kwds)
File "/usr/local/lib/python3.6/dist-packages/rasterio/__init__.py", line
216, in open
s = DatasetReader(path, driver=driver, **kwargs)
File "rasterio/_base.pyx", line 216, in
rasterio._base.DatasetBase.__init__
rasterio.errors.RasterioIOError:
netcdf:/LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc:blue_mean:
No such file or directory
Aborted!

```

But this works for both `gdalinfo` and `rio insp` tools:

```
gdalinfo hdf5:/LUSTRE/MADMEX/dir_test_mount/
madmex_003_37_-32_1996-01-01.nc://blue_mean

Driver: HDF5Image/HDF5 Dataset
Files: /LUSTRE/MADMEX/dir_test_mount/madmex_003_37_-32_1996-01-01.nc
Size is 1667, 1667
Coordinate System is `'
Metadata:
Conventions=CF-1.6, ACDD-1.3
date_created=2019-04-26T21:07:47.595416
geospatial_bounds=POLYGON ((-98.8618335626467
19.8093029971201,-98.8722811078025 19.3587577175966,-98.3949028300917
19.3481841431337,-98.382861955109 19.7986884105649,-98.8618335626467
19.8093029971201))
geospatial_bounds_crs=EPSG:4326
...
```

```
rio insp hdf5:/LUSTRE/MADMEX/dir_test_mount/
madmex_003_37_-32_1996-01-01.nc://blue_mean

/usr/local/lib/python3.6/dist-packages/rasterio/__init__.py:216:
NotGeoreferencedWarning: Dataset has no geotransform set. The identity
matrix may be returned.
s = DatasetReader(path, driver=driver, **kwargs)
Rasterio 1.0.23 Interactive Inspector (Python 3.6.7)
Type "src.meta", "src.read(1)", or "help(src)" for more information.

src.meta

{'driver': 'HDF5Image', 'dtype': 'int16', 'nodata': None, 'width': 1667,
'height': 1667, 'count': 1, 'crs': None, 'transform': Affine(1.0, 0.0, 0.0,
0.0, 1.0, 0.0)}

a=src.read()
a

array([[[632, 603, 586, ..., 505, 483, 471],
[567, 538, 536, ..., 468, 478, 523],
[614, 580, 537, ..., 540, 633, 675],
...,
[828, 810, 804, ..., 275, 268, 299],
[857, 844, 823, ..., 310, 290, 307],
[840, 854, 836, ..., 320, 288, 294]]], dtype=int16)

```

This is related with LUSTRE FS because if we copy our NetCDF file to /root
location then:

```
this works:

gdalinfo netcdf:/root/madmex_003_37_-32_1996-01-01.nc:blue_mean

Driver: netCDF/Network Common Data Format
Files: /root/madmex_003_37_-32_1996-01-01.nc
Size is 1667, 1667
Coordinate System is:
PROJCS["unnamed",
GEOGCS["WGS 84",
...
```

```
rio insp netcdf:/root/madmex_003_37_-32_1996-01-01.nc:blue_mean
Rasterio 1.0.23 Interactive Inspector (Python 3.6.7)
Type "src.meta", "src.read(1)", or "help(src)" for more information.

src.meta

{'driver': 'netCDF', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1667,
'height': 1667, 'count': 1, 'crs':
CRS.from_wkt('PROJCS["unnamed",GEOGCS["WGS
84",DATUM["unknown",SPHEROID["WGS84",6378137,6556752.3141]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",17.5],PARAMETER["standard_parallel_2",29.5],PARAMETER["latitude_of_origin",12],PARAMETER["central_meridian",-102],PARAMETER["false_easting",2500000],PARAMETER["false_northing",0]]'),
'transform': Affine(30.0, 0.0, 2827530.0,
0.0, -30.0, 876410.0)}

a=src.read()
a

array([[[632, 603, 586, ..., 505, 483, 471],
[567, 538, 536, ..., 468, 478, 523],
[614, 580, 537, ..., 540, 633, 675],
...,
[828, 810, 804, ..., 275, 268, 299],
[857, 844, 823, ..., 310, 290, 307],
[840, 854, 836, ..., 320, 288, 294]]], dtype=int16)
```

We are using `ubuntu:bionic` docker image so you can reproduce this error
using next lines:

```
sudo docker run -v /LUSTRE/:/LUSTRE/ \
-v /LUSTRE/MADMEX/dir_test_mount:/temporal \
--name mounting_volume_test_bionic \
--hostname antares3-datacube \
-dit ubuntu:bionic /bin/bash
```

enter to docker container:

```
sudo docker exec -it mounting_volume_test_bionic bash
```

and execute:

```
apt-get update

export DEBIAN_FRONTEND=noninteractive && echo "America/Mexico_City" >
/etc/timezone && apt-get install -y tzdata
apt-get update && apt-get install -y \
wget curl \
openssh-server \
openssl \
sudo \
nano \
software-properties-common \
git \
vim \
vim-gtk \
htop \
build-essential \
libssl-dev \
libffi-dev \
cmake \
python3-dev \
python3-pip \
python3-setuptools \
ca-certificates \
postgresql-client \
libudunits2-dev \
nodejs && pip3 install --upgrade pip

#install netcdf:

echo -e 'deb http://security.ubuntu.com/ubuntu/ xenial-security
universe\ndeb http://archive.ubuntu.com/ubuntu/ xenial universe' >>
/etc/apt/sources.list
apt update
apt install -t xenial-security libcurl3-gnutls
apt install -y -t xenial netcdf-bin
ncdump /LUSTRE/MADMEX/madmex_003_37_-32_1996-01-01.nc | head -n 30


#Install spatial libraries
add-apt-repository -y ppa:ubuntugis/ubuntugis-unstable && apt-get -qq
update

apt-get install -y \
ncview \
libproj-dev \
libgeos-dev \
gdal-bin \
libgdal-dev
pip install numpy && pip install GDAL==$(gdal-config --version)
--global-option=build_ext --global-option='-I/usr/include/gdal' && pip
install rasterio
export LC_ALL=C.UTF-8
export LANG=C.UTF-8
export GDAL_DATA=/usr/share/gdal/
```

Hope you can help us

Cheers,

Erick

Hi Dion,

That error comes from these lines in GDALOpenEx when the function fails to return successfully: https://github.com/OSGeo/gdal/blob/fe0b5dd644abb4ac0c9869db28e9cf977181fbce/gdal/gcore/gdaldataset.cpp#L3455. The way the errors are re-raised obscures the problem.

First, are the /vsis3/ identifiers for the datasets correct? Is Rasterio mangling them? Can you access the data successfully any time at all?

Are you using any advanced HTTP features? HTTP/2? The Rasterio wheels use a slightly old version of curl and that project is ever fixing bugs.

Is multithreading involved? I recommend that you consult https://gdal.org/gdal_vrttut.html#gdal_vrttut_mt as soon as it is back online. Multithreaded access to VRTs is complicated.

If the error occurs every time 

On Fri, May 24, 2019 at 8:13 AM Dion Häfner <dion.haefner@...> wrote:

We've only tried rasterio wheels.

On 23/05/2019 21.23, vincent.sarago via Groups.Io wrote:

This is not something I've been seen for now. does this happens with Rasterio wheels and rasterio build on gdal source ? 

--

Sean Gillies

This is not something I've been seen for now. does this happens with Rasterio wheels and rasterio build on gdal source ? 

Thank you for providing a complete example of your code and the ArcMap error message. Without those details, it can be hard to provide help.

When you open the files for writing, there is one thing missing: the *transform* keyword argument that determines the spatial extent of the dataset. I think that in your case, you will have success if you copy the value from the source datasets:

filename = "20180308_133037_1024_3B_AnalyticMS.tif"

# Copy the transform and crs objects from the source dataset.

with rasterio.open(filename) as src:

    dst_transform = src.transform

    dst_crs = src.crs

# Pass the copied values to rasterio.open().

new_dataset = rasterio.open(

       'ReflectanceB1.tif',

       'w',

       driver='GTiff',

       height=band_blue_reflectance.shape[0],

       width=band_blue_reflectance.shape[1],

       count=1,

       dtype=band_blue_reflectance.dtype,

       crs=dst_crs,

       transform=dst_transform

)    

    

new_dataset.write(band_blue_reflectance, 1)

new_dataset.close()

On Sat, May 18, 2019 at 12:39 PM <gabiklm01@...> wrote:

After following the steps in the tutorial of the Planet website (https://developers.planet.com/tutorials/convert-planetscope-imagery-from-radiance-to-reflectance/), I want to save each individual band of a RapidEye image Level 3B as a separate GeoTIFF layer file to be able to open it in any other GIS-software such as QGIS, ArcGIS other. I could not succeed in saving each band separately as GeoTIFF layer file. Here is the code:

 

```

#!/usr/bin/env python

# coding: utf-8

 

import rasterio

import numpy as np

 

filename = "20180308_133037_1024_3B_AnalyticMS.tif"

 

# Load red and NIR bands - note all PlanetScope 4-band images have band order BGRN

with rasterio.open(filename) as src:

    band_blue_radiance = src.read(1)

    

with rasterio.open(filename) as src:

    band_green_radiance = src.read(2)

 

with rasterio.open(filename) as src:

    band_red_radiance = src.read(3)

 

with rasterio.open(filename) as src:

    band_nir_radiance = src.read(4)

 

from xml.dom import minidom

 

xmldoc = minidom.parse("20180308_133037_1024_3B_AnalyticMS_metadata.xml")

nodes = xmldoc.getElementsByTagName("ps:bandSpecificMetadata")

 

# XML parser refers to bands by numbers 1-4

coeffs = {}

for node in nodes:

    bn = node.getElementsByTagName("ps:bandNumber")[0].firstChild.data

    if bn in ['1', '2', '3', '4']:

        i = int(bn)

        value = node.getElementsByTagName("ps:reflectanceCoefficient")[0].firstChild.data

        coeffs[i] = float(value)

 

print "Conversion coefficients:", coeffs

 

# Multiply the Digital Number (DN) values in each band by the TOA reflectance coefficients

band_blue_reflectance = band_blue_radiance * coeffs[1]

band_green_reflectance = band_green_radiance * coeffs[2]

band_red_reflectance = band_red_radiance * coeffs[3]

band_nir_reflectance = band_nir_radiance * coeffs[4]

 

import numpy as np

print "Red band radiance is from {} to {}".format(np.amin(band_red_radiance), np.amax(band_red_radiance))

print "Red band reflectance is from {} to {}".format(np.amin(band_red_reflectance), np.amax(band_red_reflectance))

 

print "Before Scaling, red band reflectance is from {} to {}".format(np.amin(band_red_reflectance), np.amax(band_red_reflectance))

 

# Here we include a fixed scaling factor. This is common practice.

scale = 10000

blue_ref_scaled = scale * band_blue_reflectance

green_ref_scaled = scale * band_green_reflectance

red_ref_scaled = scale * band_red_reflectance

nir_ref_scaled = scale * band_nir_reflectance

 

print "After Scaling, red band reflectance is from {} to {}".format(np.amin(red_ref_scaled), np.amax(red_ref_scaled))

```

Here I tried unsuccessfully to save each individual band as GeotiFF file.

 

I have tried, any idea how to do the save of each band correctly to further open in QGIS orArGIS for calculating NDVI there?

 

```

dst_crs='EPSG:32720'

 

##  saving bands individualy  

new_dataset = rasterio.open(

       'ReflectanceB1.tif',

       'w',

       driver='GTiff',

       height=band_blue_reflectance.shape[0],

       width=band_blue_reflectance.shape[1],

       count=1,

       dtype=band_blue_reflectance.dtype,

       crs=dst_crs,

)    

    

new_dataset.write(band_blue_reflectance, 1)

new_dataset.close()

   

   

new_dataset = rasterio.open(

       'ReflectanceB2.tif',

       'w',

       driver='GTiff',

       height=band_red_reflectance.shape[0],

       width=band_red_reflectance.shape[1],

       count=1,

       dtype=band_red_reflectance.dtype,

       crs=dst_crs,

)    

    

new_dataset.write(band_red_reflectance, 1)

new_dataset.close()

 

new_dataset = rasterio.open(

       'ReflectanceB3.tif',

       'w',

       driver='GTiff',

       height=band_red_reflectance.shape[0],

       width=band_red_reflectance.shape[1],

       count=1,

       dtype=band_red_reflectance.dtype,

       crs=dst_crs,

)    

    

new_dataset.write(band_red_reflectance, 1)

new_dataset.close()    

 

 

new_dataset = rasterio.open(

       'ReflectanceB4.tif',

       'w',

       driver='GTiff',

       height=band_red_reflectance.shape[0],

       width=band_red_reflectance.shape[1],

       count=1,

       dtype=band_red_reflectance.dtype,

       crs=dst_crs,

)    

    

new_dataset.write(band_red_reflectance, 1)

new_dataset.close()

```

I'm using as input image "20180308_133037_1024_3B_AnalyticMS.tif" instead of "20180308_133037_1024_3B_AnalyticMS_SR.tif" don't know if is the problem

 

I'm getting resulting GeoTIFF file but when I open it in ArcMap I get this error

 

The GeoTIFF band are there and can be seen at ArcMap screen...   but the error appear anyway... might be related to coordinate systems to be defined possibly while saving in rasterio? And need to have it correctly projected for clipping with shapefiles.

--

Sean Gillies