I would recommend using a wrapper around rasterio
for this task called rioxarray.
Here is how you would do the calculations:
#!/usr/bin/env python
# coding: utf-8
import numpy as np
import rioxarray
import xarray
filename = "20180308_133037_1024_3B_AnalyticMS.tif"
# Load red and NIR bands - note all PlanetScope 4-band images have band order BGRN
xds = xarray.open_rasterio(filename)
band_blue_radiance = xds.sel(band=1)
band_green_radiance = xds.sel(band=2)
band_red_radiance = xds.sel(band=3)
band_nir_radiance = xds.sel(band=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]
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)))
And here is how you would write it out to a raster:
def to_raster(in_xds, template_xds, out_file):
in_xds = in_xds.rio.write_crs(template_xds.rio.crs)
if template_xds.rio.nodata is not None:
in_xds.attrs["_FillValue"] = template_xds.rio.nodata
in_xds.rio.to_raster(out_file)
to_raster(blue_ref_scaled, xds, "ReflectanceB1.tif")
to_raster(green_ref_scaled, xds, "ReflectanceB2.tif")
to_raster(red_ref_scaled, xds, "ReflectanceB3.tif")
to_raster(nir_ref_scaled, xds, "ReflectanceB4.tif")
Best Answer
The previous answers are misleading or wrong. To modify the nodata value of a GeoTIFF with Rasterio, do this
The project has tests of this usage that you may see also: https://github.com/mapbox/rasterio/blob/master/tests/test_update.py#L59-L64. Note that you may have to close and reopen the file (as done in the test) to see the nodata value take effect in your program.
The
meta
property of a dataset is a copy of some of its important metadata. Modifying that object has no effect on the dataset.