I am working with LandSat 8 data, calculating a vegetation index from a single scene. I can read the necessary bands and do the calculation. My issue occurs when I try to save the information back to a GeoTiff.
None of the spatial characteristics of the data have changed - same height, width, count (which is 1 for LandSat data files - each band is in its own file), crs, and transform.
And yet, when I save this back to a GeoTiff, the resulting code file only contains zeros.
To help folks reproduce this, below are some example cloud-optimized files:
cogs = ['s3://usgs-landsat/collection02/level-2/standard/oli-tirs/2019/026/030/LC08_L2SP_026030_20190505_20200828_02_T1/LC08_L2SP_026030_20190505_20200828_02_T1_SR_B4.TIF',
's3://usgs-landsat/collection02/level-2/standard/oli-tirs/2019/026/030/LC08_L2SP_026030_20190505_20200828_02_T1/LC08_L2SP_026030_20190505_20200828_02_T1_SR_B5.TIF',
's3://usgs-landsat/collection02/level-2/standard/oli-tirs/2019/026/030/LC08_L2SP_026030_20190505_20200828_02_T1/LC08_L2SP_026030_20190505_20200828_02_T1_QA_PIXEL.TIF']
Import needed packages
import matplotlib.pyplot as plt
import numpy as np
import rasterio as rio
from rasterio.session import AWSSession
import boto3
Read in the data as variables
aws_session = AWSSession(boto3.Session(region_name='us-west-2'), requester_pays=True)
with rio.Env(aws_session,
AWS_NO_SIGN_REQUEST='NO',
GDAL_DISABLE_READDIR_ON_OPEN='TRUE'):
with rio.open(cogs[0]) as src:
profile = src.profile
red_raw = src.read(1)
with rio.open(cogs[1]) as src_nir:
nir_raw = src_nir.read(1)
with rio.open(cogs[2]) as src_qa:
qa = src_qa.read(1)
Apply the data and cloud mask, scale the good quality data
cloud_free_vals = [21824, 21952, 22080, 22208, 23888, 24144, 30048, 30304,
54596, 54724, 54852, 54980, 56660, 56916, 62820, 63076]
cloud_mask = np.isin(qa, cloud_free_vals)
red = np.ma.masked_equal(red_raw,0.0)
red = np.ma.array(red, mask=cloud_mask)
red = red * 0.0000275 + -0.2
red[red>1] = 1
red[red<0] = 0
nir = nir_raw * 0.0000275 + -0.2
nir = np.ma.masked_less(nir, 0.0)
nir[nir>1] = 1
nirv = (nir-red)/(nir+red)*nir
When I plot the results, I get something reasonable
masked LandSat data.
But when I try and save the result as a geotiff to be accessed later (so that we preserve the geographical information - I have a lot of files to process!) I only end up saving a file full of zeros.
Here's the code I used to save the geotiff, using the profile saved from the red band file above.
save_example_file = rio.open(
'test.tif',
'w',
**profile)
save_example_file.write(nirv.data, 1)
save_example_file.close()
Opening it again like so yields a field of zeros:
test2 = rio.open('test.tif', 'r')
plot_example = test2.read(1)
plt.figure()
plt.imshow(plot2)
plt.colorbar()
plt.show()
Which results in Bunch of zeros!
