Date   

Re: masking in rio calc

Gregory, Matthew
 

Great, thanks Sean.

 

As I’m sure you figured out from my original example, but I wasn’t explicit about, was that filling with 255 when the numpy dtype was bool resulted in True, which then gets set to 1 on the cast to uint8.  Therefore, all input mask pixels become 1 in the output.

 

matt

 

-----------------------

 

Update: https://github.com/mapbox/rasterio/issues/2041. Thanks for bringing this up, Matt.

 

On Wed, Nov 18, 2020 at 1:54 PM Sean Gillies <sean.gillies@...> wrote:

Hi Matt,

 

On Tue, Nov 17, 2020 at 5:57 PM Gregory, Matthew <matt.gregory@...> wrote:

Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do.  I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked.  I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output.  My command is:

  rio calc
    -f GTiff
    -t uint8
    --overwrite
    --masked
    --profile nodata=255
    "(>= (read 1) 0.5)"
    raw.tif
    threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype.  When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result.  But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

  res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt

 

You may get better results using numpy.where, like "(where (>= (read 1) 0.5) 1 0)".

 

I'll dig into the filling code, that looks like a bug to me.

 

--

Sean Gillies


Re: masking in rio calc

Sean Gillies
 

Update: https://github.com/mapbox/rasterio/issues/2041. Thanks for bringing this up, Matt.

On Wed, Nov 18, 2020 at 1:54 PM Sean Gillies <sean.gillies@...> wrote:
Hi Matt,

On Tue, Nov 17, 2020 at 5:57 PM Gregory, Matthew <matt.gregory@...> wrote:
Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do.  I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked.  I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output.  My command is:

  rio calc
    -f GTiff
    -t uint8
    --overwrite
    --masked
    --profile nodata=255
    "(>= (read 1) 0.5)"
    raw.tif
    threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype.  When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result.  But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

  res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt

You may get better results using numpy.where, like "(where (>= (read 1) 0.5) 1 0)".

I'll dig into the filling code, that looks like a bug to me.

--
Sean Gillies


Re: masking in rio calc

Sean Gillies
 

Hi Matt,

On Tue, Nov 17, 2020 at 5:57 PM Gregory, Matthew <matt.gregory@...> wrote:
Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do.  I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked.  I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output.  My command is:

  rio calc
    -f GTiff
    -t uint8
    --overwrite
    --masked
    --profile nodata=255
    "(>= (read 1) 0.5)"
    raw.tif
    threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype.  When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result.  But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

  res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt

You may get better results using numpy.where, like "(where (>= (read 1) 0.5) 1 0)".

I'll dig into the filling code, that looks like a bug to me.

--
Sean Gillies


masking in rio calc

Gregory, Matthew
 

Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do. I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked. I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output. My command is:

rio calc
-f GTiff
-t uint8
--overwrite
--masked
--profile nodata=255
"(>= (read 1) 0.5)"
raw.tif
threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype. When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result. But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt


Re: Sampling from one single point location using sample()

Alexander Jüstel
 

Hey Sean, 

I made it work using your tests! Thanks for that ;)

Cheers Alex

On Mon, Nov 16, 2020 at 4:03 PM Sean Gillies via groups.io <sean=mapbox.com@groups.io> wrote:
Hi Alex,

On Sun, Nov 15, 2020 at 9:06 AM Alexander Jüstel <alexander.juestel@...> wrote:
Hello, 
 
I would like to sample the value of my raster at one single point location. I have passed the coordinates as an iterable to the sample() function. However, I am not able to extract the value from the generator. The sampling works fine for more than one point. However, I need it to work for one point only. I have tried using list(gen) or next(gen) to get the value but nothing works so far.
 
Thanks for your support
Cheers
Alex

Is your approach any different than the examples in https://github.com/mapbox/rasterio/blob/master/tests/test_sampling.py?

--
Sean Gillies


Re: Sampling from one single point location using sample()

Sean Gillies
 

Hi Alex,

On Sun, Nov 15, 2020 at 9:06 AM Alexander Jüstel <alexander.juestel@...> wrote:
Hello, 
 
I would like to sample the value of my raster at one single point location. I have passed the coordinates as an iterable to the sample() function. However, I am not able to extract the value from the generator. The sampling works fine for more than one point. However, I need it to work for one point only. I have tried using list(gen) or next(gen) to get the value but nothing works so far.
 
Thanks for your support
Cheers
Alex

Is your approach any different than the examples in https://github.com/mapbox/rasterio/blob/master/tests/test_sampling.py?

--
Sean Gillies


Re: Is it possible to ignore existing overview when performing decimated read?

Sean Gillies
 

Hi Loïc,

A test that requires 3.3 would be welcome. 

The next time we have a rasterio release I will patch GDAL the GDAL library in the wheels we publish on PyPI.


On Mon, Nov 16, 2020, 3:26 AM Loïc Dutrieux <loic.dutrieux@...> wrote:
Many thanks Sean and Even,

So with rasterio compiled against gdal 3.3 I'll be able to force ignore overviews?
If that deserves a new rasterio unit test I'm happy to send a pull request in the coming days.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Sean Gillies via groups.io <sean=mapbox.com@groups.io>
Sent: 13 November 2020 20:57
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

Hi Even,

On Fri, Nov 13, 2020 at 10:38 AM Even Rouault <even.rouault@...<mailto:even.rouault@...>> wrote:
On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io<https://urldefense.com/v3/__http://groups.io__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTeRRzVSk$> wrote:
> HI Loïc,
>
> I don't think we currently have any option other than to remove the
> overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
> the implementation isn't there. For example, if you used gdalwarp with
> `-ovr NONE` it will set an internal overview level to the value of -1. But
> that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
> `OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.

Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181<https://urldefense.com/v3/__https://github.com/OSGeo/gdal/pull/3181__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTK8LUNCg$>

Even

Thanks


Re: Is it possible to ignore existing overview when performing decimated read?

Loïc Dutrieux
 

Many thanks Sean and Even,

So with rasterio compiled against gdal 3.3 I'll be able to force ignore overviews?
If that deserves a new rasterio unit test I'm happy to send a pull request in the coming days.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Sean Gillies via groups.io <sean=mapbox.com@groups.io>
Sent: 13 November 2020 20:57
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

Hi Even,

On Fri, Nov 13, 2020 at 10:38 AM Even Rouault <even.rouault@...<mailto:even.rouault@...>> wrote:
On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io<https://urldefense.com/v3/__http://groups.io__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTeRRzVSk$> wrote:
HI Loïc,

I don't think we currently have any option other than to remove the
overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
the implementation isn't there. For example, if you used gdalwarp with
`-ovr NONE` it will set an internal overview level to the value of -1. But
that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
`OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.
Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181<https://urldefense.com/v3/__https://github.com/OSGeo/gdal/pull/3181__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTK8LUNCg$>

Even

Thanks!

--
Sean Gillies


Sampling from one single point location using sample()

Alexander Jüstel
 

Hello, 
 
I would like to sample the value of my raster at one single point location. I have passed the coordinates as an iterable to the sample() function. However, I am not able to extract the value from the generator. The sampling works fine for more than one point. However, I need it to work for one point only. I have tried using list(gen) or next(gen) to get the value but nothing works so far.
 
Thanks for your support
Cheers
Alex


Re: Is it possible to ignore existing overview when performing decimated read?

Sean Gillies
 

Hi Even,

On Fri, Nov 13, 2020 at 10:38 AM Even Rouault <even.rouault@...> wrote:
On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io wrote:
> HI Loïc,
>
> I don't think we currently have any option other than to remove the
> overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
> the implementation isn't there. For example, if you used gdalwarp with
> `-ovr NONE` it will set an internal overview level to the value of -1. But
> that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
> `OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.

Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181

Even

Thanks!

--
Sean Gillies


Re: Is it possible to ignore existing overview when performing decimated read?

Even Rouault
 

On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io wrote:
HI Loïc,

I don't think we currently have any option other than to remove the
overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
the implementation isn't there. For example, if you used gdalwarp with
`-ovr NONE` it will set an internal overview level to the value of -1. But
that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
`OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.
Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com


Re: Is it possible to ignore existing overview when performing decimated read?

Sean Gillies
 

HI Loïc,

I don't think we currently have any option other than to remove the overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but the implementation isn't there. For example, if you used gdalwarp with `-ovr NONE` it will set an internal overview level to the value of -1. But that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass `OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.

To illustrate

>>> rasterio.open("tests/data/green.tif").shape
(64, 64)
>>> rasterio.open("tests/data/green.tif", OVERVIEW_LEVEL=0).shape
(32, 32)
>>> rasterio.open("tests/data/green.tif", OVERVIEW_LEVEL="NONE").shape
(32, 32)
>>> rasterio.open("tests/data/green.tif", OVERVIEW_LEVEL=-1).shape
Traceback (most recent call last):
  File "rasterio/_base.pyx", line 261, in rasterio._base.DatasetBase.__init__
    self._hds = open_dataset(filename, flags, driver, kwargs, None)
  File "rasterio/_shim.pyx", line 78, in rasterio._shim.open_dataset
    return exc_wrap_pointer(hds)
  File "rasterio/_err.pyx", line 213, in rasterio._err.exc_wrap_pointer
    raise exc
rasterio._err.CPLE_OpenFailedError: Cannot open overview level -1 of tests/data/green.tif

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sean/projects/rasterio/rasterio/env.py", line 433, in wrapper
    return f(*args, **kwds)
  File "/home/sean/projects/rasterio/rasterio/__init__.py", line 221, in open
    s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
  File "rasterio/_base.pyx", line 263, in rasterio._base.DatasetBase.__init__
    raise RasterioIOError(str(err))
rasterio.errors.RasterioIOError: Cannot open overview level -1 of tests/data/green.tif


On Fri, Nov 13, 2020 at 4:33 AM Loïc Dutrieux <loic.dutrieux@...> wrote:
Hi Denis,

Thanks for pointing to that issue. I may have misunderstood it but I don't think it helps answering my question.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Denis Rykov <rykovd@...>
Sent: 10 November 2020 21:36:13
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

You can find an example how to do that here: https://github.com/mapbox/rasterio/issues/1929<https://urldefense.com/v3/__https://github.com/mapbox/rasterio/issues/1929__;!!DOxrgLBm!XaInc3_J0VCLXQo6gYXo0rBWHs4b0p4ZkMuimU2g9HtJy0eRb55BN-J7JDPt-MB3Dhg5K3E$>.

On Tue, Nov 10, 2020, 8:12 PM Loïc Dutrieux <loic.dutrieux@...<mailto:loic.dutrieux@...>> wrote:
Hi everyone,

I'm trying to perform a decimated read of a dataset that already contains overviews. Regardless of which value I pass to the resampling= argument of the read method, it seems that the already existing overview is used. Is there any way to ignore it? See the reproducible example below.

Cheers,
Loïc

---

import tempfile
import os

import numpy as np
import rasterio
from rasterio.enums import Resampling


filename = os.path.join(tempfile.gettempdir(), 'overview_test.tif')

# Create random int array
shape = (100, 100)
arr = np.random.randint(0, 9999, size=shape, dtype=np.uint16)

meta = {'height': 100,
        'width': 100,
        'driver': 'GTiff',
        'dtype': np.uint16,
        'count': 1}

# Write random array to file and compute first overview
with rasterio.open(filename, 'w', **meta) as dst:
    dst.write(arr, 1)
    dst.build_overviews([2], Resampling.nearest)

# Read with downsample
with rasterio.open(filename) as src:
    arr_avg = src.read(1, out_shape=(1,50,50), resampling=Resampling.average,
                                  out_dtype=np.float)
    arr_nrt = src.read(1, out_shape=(1,50,50), resampling=Resampling.nearest)

print(np.max(arr_avg - arr_nrt))
# when the source file does not contain overviews, the max of the difference array
# is > 0


--
Sean Gillies


Re: Is it possible to ignore existing overview when performing decimated read?

Loïc Dutrieux
 

Hi Denis,

Thanks for pointing to that issue. I may have misunderstood it but I don't think it helps answering my question.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Denis Rykov <rykovd@...>
Sent: 10 November 2020 21:36:13
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

You can find an example how to do that here: https://github.com/mapbox/rasterio/issues/1929<https://urldefense.com/v3/__https://github.com/mapbox/rasterio/issues/1929__;!!DOxrgLBm!XaInc3_J0VCLXQo6gYXo0rBWHs4b0p4ZkMuimU2g9HtJy0eRb55BN-J7JDPt-MB3Dhg5K3E$>.

On Tue, Nov 10, 2020, 8:12 PM Loïc Dutrieux <@loicdtx<mailto:@loicdtx>> wrote:
Hi everyone,

I'm trying to perform a decimated read of a dataset that already contains overviews. Regardless of which value I pass to the resampling= argument of the read method, it seems that the already existing overview is used. Is there any way to ignore it? See the reproducible example below.

Cheers,
Loïc

---

import tempfile
import os

import numpy as np
import rasterio
from rasterio.enums import Resampling


filename = os.path.join(tempfile.gettempdir(), 'overview_test.tif')

# Create random int array
shape = (100, 100)
arr = np.random.randint(0, 9999, size=shape, dtype=np.uint16)

meta = {'height': 100,
'width': 100,
'driver': 'GTiff',
'dtype': np.uint16,
'count': 1}

# Write random array to file and compute first overview
with rasterio.open(filename, 'w', **meta) as dst:
dst.write(arr, 1)
dst.build_overviews([2], Resampling.nearest)

# Read with downsample
with rasterio.open(filename) as src:
arr_avg = src.read(1, out_shape=(1,50,50), resampling=Resampling.average,
out_dtype=np.float)
arr_nrt = src.read(1, out_shape=(1,50,50), resampling=Resampling.nearest)

print(np.max(arr_avg - arr_nrt))
# when the source file does not contain overviews, the max of the difference array
# is > 0


Re: Rasterizing polygon from numpy arrays

alexisshakas@...
 

Works like a charm, thank you!


Re: Rasterizing polygon from numpy arrays

Brendan Ward
 

Yes, you should be able to pass a list / tuple (these are iterables) of Shapely geometries as input to geometry_mask.


Re: Rasterizing polygon from numpy arrays

alexisshakas@...
 

Hi Brendan and thnk you for the quick response!

The coordinates should not be an issue, since the polygon coordinates are in the same frame as the image.

Can i directly input a Shapely polygon to the geometry_mask class? I was a bit lost on whether this needs to be in some other format, like a sort of iterable dictionary.


Re: Rasterizing polygon from numpy arrays

Brendan Ward
 

You need to convert each polygon in your numpy array of coordinates into a GeoJSON-like polygon.

You could use Shapely to do so.

You should then be able to use those as input to rasterize or geometry_mask.  Note that your polygons must be in the coordinate system of your images in order for this to work.  Since it sounds like your x,y dimensions of your images are in integer degrees, this means your x,y coordinates must also be in degrees (vertical axis shouldn't matter for rasterize / geometry_mask).

Hope that helps!


Rasterizing polygon from numpy arrays

alexisshakas@...
 

Hello,
I've been struggling with this for quite some time now, so its time to ask for help.

My problem is set up as follows:
I want to create a binary mask for an image that is quite large (~10k rows and 360 columns). I have several polygons, all closed, and I want to assign all values that fall in the polygons as True.
The polygons are in the form of a numpy array, with (x,y) coordinates. Each polygon has roughly 1000 values. I have about 200 such polygons for a given image.
Ideally, since the polygons are quite small, I should not query all points of the image (as suggested for example here https://stackoverflow.com/questions/36399381/whats-the-fastest-way-of-checking-if-a-point-is-inside-a-polygon-in-python)



I found documentation about how to rasterize shapefiles here https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html
using for example
rasterio.features.geometry_mask(geometries, out_shape, transform, all_touched=False, invert=False)
but I could not figure out how to do this with a simple numpy array.

Is it possible?

I should also mention that my image axes are mixed. The vertical axis is a float (depths, from 0m to ~100m at 1 mm spacing) while the horizontal axis contains integers (angles, from 0 to 360 in 1 degree increments).

But I could also change these to be simply indices (both integers).

Thank you for any feedback!



Re: Is it possible to ignore existing overview when performing decimated read?

Denis Rykov
 

You can find an example how to do that here: https://github.com/mapbox/rasterio/issues/1929.


On Tue, Nov 10, 2020, 8:12 PM Loïc Dutrieux <loic.dutrieux@...> wrote:
Hi everyone,

I'm trying to perform a decimated read of a dataset that already contains overviews. Regardless of which value I pass to the resampling= argument of the read method, it seems that the already existing overview is used. Is there any way to ignore it? See the reproducible example below.

Cheers,
Loïc

---

import tempfile
import os

import numpy as np
import rasterio
from rasterio.enums import Resampling


filename = os.path.join(tempfile.gettempdir(), 'overview_test.tif')

# Create random int array
shape = (100, 100)
arr = np.random.randint(0, 9999, size=shape, dtype=np.uint16)

meta = {'height': 100,
        'width': 100,
        'driver': 'GTiff',
        'dtype': np.uint16,
        'count': 1}

# Write random array to file and compute first overview
with rasterio.open(filename, 'w', **meta) as dst:
    dst.write(arr, 1)
    dst.build_overviews([2], Resampling.nearest)

# Read with downsample
with rasterio.open(filename) as src:
    arr_avg = src.read(1, out_shape=(1,50,50), resampling=Resampling.average,
                                  out_dtype=np.float)
    arr_nrt = src.read(1, out_shape=(1,50,50), resampling=Resampling.nearest)

print(np.max(arr_avg - arr_nrt))
# when the source file does not contain overviews, the max of the difference array
# is > 0









Is it possible to ignore existing overview when performing decimated read?

Loïc Dutrieux
 

Hi everyone,

I'm trying to perform a decimated read of a dataset that already contains overviews. Regardless of which value I pass to the resampling= argument of the read method, it seems that the already existing overview is used. Is there any way to ignore it? See the reproducible example below.

Cheers,
Loïc

---

import tempfile
import os

import numpy as np
import rasterio
from rasterio.enums import Resampling


filename = os.path.join(tempfile.gettempdir(), 'overview_test.tif')

# Create random int array
shape = (100, 100)
arr = np.random.randint(0, 9999, size=shape, dtype=np.uint16)

meta = {'height': 100,
        'width': 100,
        'driver': 'GTiff',
        'dtype': np.uint16,
        'count': 1}

# Write random array to file and compute first overview
with rasterio.open(filename, 'w', **meta) as dst:
    dst.write(arr, 1)
    dst.build_overviews([2], Resampling.nearest)

# Read with downsample
with rasterio.open(filename) as src:
    arr_avg = src.read(1, out_shape=(1,50,50), resampling=Resampling.average,
                        out_dtype=np.float)
    arr_nrt = src.read(1, out_shape=(1,50,50), resampling=Resampling.nearest)

print(np.max(arr_avg - arr_nrt))
# when the source file does not contain overviews, the max of the difference array
# is > 0