0

I have Landsat-8 OLI images in directories named after image WRS-2 scenes. I tried to compose script that would iterate these directories and stack rasters within each directory - using earthpy python package.

earthpy.spatial.stack() -function returns array and a rasterio object, and writes a new stacked .tif file to the folder. Currently the function returns ValueError: the array's dtype 'uint8' does not match the files dtype 'uint16'. Is it possible to work around this problem without scaling the uint16 data to uint8 format?

my script:


# Import needed packages

import glob2 as glob
import rasterio
import os
import earthpy
import matplotlib.pyplot as plt
import rasterio as rio
from rasterio.plot import plotting_extent
import geopandas as gpd
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep

# WRS-2 naming convention: path 210, row 4 = 210004)
wrs = ["210004", "211004", "212004", "213004", "214004","215004","216004","217003","217004","218003","219003"]

#Local path for reference images 
refPath = "F:/ref_Images/"

for sceneId in wrs:
    ''' Iteratively create datacubes by stacking bands within each acquistion scene folder''' 
    #check sub-directories (acquistions) under every WRS-scene
    subDirs = os.listdir(refPath + sceneId + "/")
    for directory in subDirs:
        #iterate subdirectories and compile datacubes within them 
        bandList = glob.glob(refPath + sceneId + "/" + directory + "/*.tif")
        output = refPath + sceneId + "/" + directory + "/"
        filename = directory + "_STACK.tif"
        array, raster = es.stack(bandList, out_path=output + filename, nodata = 0)

Value error returned: 

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-82fa289db5a0> in <module>
     10         output = refPath + sceneId + "/" + directory + "/"
     11         filename = directory + "_STACK.tif"
---> 12         array, raster = es.stack(bandList, out_path=output + filename, nodata = 0)

C:\ProgramData\Anaconda3\envs\lopputyo\lib\site-packages\earthpy\spatial.py in stack(band_paths, out_path, nodata)
    252             # Write stacked gtif file
    253             with rio.open(out_path, "w", **dest_kwargs) as dest:
--> 254                 _stack_bands(sources, write_raster, dest)
    255 
    256             # Read and return array

C:\ProgramData\Anaconda3\envs\lopputyo\lib\site-packages\earthpy\spatial.py in _stack_bands(sources, write_raster, dest)
    313                 bands = bands[np.newaxis, ...]
    314             for band in bands:
--> 315                 dest.write(band, ii + 1)
    316 
    317     else:

rasterio\_io.pyx in rasterio._io.DatasetWriterBase.write()

ValueError: the array's dtype 'uint8' does not match the file's dtype 'uint16'
PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Mitro Müller
  • 71
  • 8
  • Would you be willing to use other alternatives to earthpy? I did not find anything about this error in the documentation so you'd probably need to modify the function in the source code itself. A similar function can be easily written using numpy and either gdal or rasterio. – Marcelo Villa Jan 14 '20 at 19:39
  • Alternative suggestions are welcome. Problem just surprised me, as some earthpy examples use Landsat-8 data. Like this one: https://earthpy.readthedocs.io/en/latest/gallery_vignettes/plot_raster_stack_crop.html#id1 – Mitro Müller Jan 15 '20 at 05:22
  • I have not found how to change the datatype directly in earthpy.spatial.stack() -function, or why the function works at earthpy docs examples. As this function utlizes rasterio, answer how to force written data to uint16 data-format has been presented in: https://gis.stackexchange.com/questions/284179/add-bands-name-and-description-to-the-metadata-when-stacking-using-rasterio – Mitro Müller Jan 15 '20 at 08:11

0 Answers0