3

Is there a way to preserve the mask of a numpy array when reprojecting it with rasterio? See the reproducible example with real data below:

import os
import tempfile

import rasterio
from rasterio.crs import CRS
from rasterio.warp import reproject, calculate_default_transform
import matplotlib.pyplot as plt
import numpy as np
import requests

print(rasterio.__version__) #1.0a10

# Download function, only downloads the file once
def download(url):
    filename = os.path.join(tempfile.gettempdir(), os.path.basename(url))
    if not os.path.isfile(filename):
        r = requests.get(url, stream=True)
        with open(filename, 'wb') as dst:
            for chunk in r.iter_content(chunk_size=1024):
                dst.write(chunk)
    return filename

# Download test data (14MB)
scene_id = 'LC08_L1TP_026047_20170616_20170629_01_T1'
url = 'https://earthexplorer.usgs.gov/browse/gisready/landsat_8/%s.zip' % scene_id
archive = download(url)

# Read data as numpy array and apply mask
rgb = '/vsizip/%s/%s.tif' % (archive, scene_id)
mask = '/vsizip/%s/%s_QB.tif' % (archive, scene_id)
with rasterio.open(mask) as src:
    m = src.read(1)
with rasterio.open(rgb) as src:
    src_crs = src.crs
    src_transform = src.transform
    src_bounds = src.bounds
    src_width = src.width
    src_height = src.height
    red_array = np.ma.array(src.read(1), mask=m!=0)

# VIsualize masked array (clouds and cloud shadows are properly masked)
# plt.imshow(red_array)
# plt.show()

# Prepare destination array for warping
dst_crs = CRS({u'proj': u'laea', u'lon_0': -106, u'lat_0': 24}) # Equal area crs more or less centered on mexico center
dst_transform, dst_width, dst_height = calculate_default_transform(src_crs,
                                                                dst_crs,
                                                                src_width,
                                                                src_height,
                                                                *src_bounds,
                                                                resolution=(30, 30))
red_laea = np.ma.zeros((dst_height, dst_width), dtype=np.uint8)

# reproject array
reproject(source=red_array,
          destination=red_laea,
          src_transform=src_transform,
          src_crs=src_crs,
          dst_transform=dst_transform,
          dst_crs=dst_crs,
          dst_nodata=0)

# Visualize (clouds are no longer masked)
plt.imshow(red_laea)
plt.show()

What I would like/expect is that the mask (cloud and shadows in the example above) is transferred to the destination array. I was sort of expecting that behaviour from the description of the src_nodata argument (rasterio.warp.reproject function), but I may have misinterpreted it.

2

Very nice example. I put it in a notebook here. I can reproduce your issue with rasterio 1.0a11.

Problems

Rasterio's reproject does intend to support masked arrays, but it seems to fail in this case.

  1. One problem is that the .fill_value of your masked array is used as src_nodata (see these lines). That is wrong, because the actual .data in the masked array underneath the mask might have different values than .fill_value. But even filling the data underneath does not help.

  2. Also, your red_array.fill_value is 999999, which is actually too large for the red_array.dtype, which is uint8. That is a big problem with NumPy's masked arrays problem. Even worse: taking a .filled() copy, you do not get an error but 63...

  3. The reproject function will not fill your destination with np.ma.masked. Especially not when you set dst_nodata=0. dst_nodata=np.ma.masked would be better, but fails due to dtype checking. So you will need to create a masked array from the returned data with something like

    red_laea[red_laea == dst_nodata] = np.ma.masked
    

In summary

I can understand that you think you can use masked arrays with rasterio.warp.reproject, but it is not yet implemented well and masked arrays can have issues that are not easy to see (invalid .fill_value). So I recommend that you use regular arrays and keep track of the nodata value yourself. And consider raising an issue on rasterio.

masked

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.