1

After following the steps in the tutorial of the Planet website (https://developers.planet.com/tutorials/convert-planetscope-imagery-from-radiance-to-reflectance/), I want to save each individual band of a RapidEye image Level 3B as a separate GeoTIFF layer file to be able to open it in any other GIS-software such as QGIS, ArcGIS other. I could not succeed in saving each band separately as GeoTIFF layer file. Here is the code:

#!/usr/bin/env python
# coding: utf-8

import rasterio
import numpy as np

filename = "20180308_133037_1024_3B_AnalyticMS.tif"

# Load red and NIR bands - note all PlanetScope 4-band images have band order BGRN
with rasterio.open(filename) as src:
    band_blue_radiance = src.read(1)

with rasterio.open(filename) as src:
    band_green_radiance = src.read(2)

with rasterio.open(filename) as src:
    band_red_radiance = src.read(3)

with rasterio.open(filename) as src:
    band_nir_radiance = src.read(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]

import numpy as np
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))

Here I tried unsuccessfully to save each individual band as GeoTIFF file.

I have tried, any idea how to do the save of each band correctly to further open in QGIS orArGIS for calculating NDVI there?

dst_crs='EPSG:32720'

##  saving bands individualy  
new_dataset = rasterio.open(
       'ReflectanceB1.tif',
       'w',
       driver='GTiff',
       height=band_blue_reflectance.shape[0],
       width=band_blue_reflectance.shape[1],
       count=1,
       dtype=band_blue_reflectance.dtype,
       crs=dst_crs,
)    

new_dataset.write(band_blue_reflectance, 1)
new_dataset.close()


new_dataset = rasterio.open(
       'ReflectanceB2.tif',
       'w',
       driver='GTiff',
       height=band_red_reflectance.shape[0],
       width=band_red_reflectance.shape[1],
       count=1,
       dtype=band_red_reflectance.dtype,
       crs=dst_crs,
)    

new_dataset.write(band_red_reflectance, 1)
new_dataset.close()

new_dataset = rasterio.open(
       'ReflectanceB3.tif',
       'w',
       driver='GTiff',
       height=band_red_reflectance.shape[0],
       width=band_red_reflectance.shape[1],
       count=1,
       dtype=band_red_reflectance.dtype,
       crs=dst_crs,
)    

new_dataset.write(band_red_reflectance, 1)
new_dataset.close()    


new_dataset = rasterio.open(
       'ReflectanceB4.tif',
       'w',
       driver='GTiff',
       height=band_red_reflectance.shape[0],
       width=band_red_reflectance.shape[1],
       count=1,
       dtype=band_red_reflectance.dtype,
       crs=dst_crs,
)    

new_dataset.write(band_red_reflectance, 1)
new_dataset.close()

I'm using as input image "20180308_133037_1024_3B_AnalyticMS.tif" instead of "20180308_133037_1024_3B_AnalyticMS_SR.tif" don't know if is the problem

I'm getting resulting GeoTIFF file but when I open it in ArcMap I get this error

enter image description here

The GeoTIFF band are there and can be seen at ArcMap screen... but the error appear anyway... might be related to coordinate systems to be defined possibly while saving in rasterio? And need to have it correctly projected for clipping with shapefiles.

I have done the change in data frame properties before loading the raster layer as

enter image description here

But the same worning appears

In Jupyter notebook I got the following warning...

enter image description here

  • That is not an error, but a warning. You need to report the canvas coordinate reference -- it could be the data added previously that caused the problem. – Vince May 18 '19 at 17:19
  • how can I do that? thakns – Gab May 18 '19 at 17:38
  • Vince can you give me a hint on how to do that? – Gab May 18 '19 at 22:27
  • Not in two of the three products you've got listed. An ArcMap canvas has Data Frame properties through a right click on the LayerList. You need to provide individual spatial reference properties for all the objects as well. This is the core of your problem, but it's not described anywhere in the question. – Vince May 18 '19 at 22:38
  • Hi Vince, see above I tried to define the data frame you mention but the warning still appears.... the warning say something about the extent, may it be possible that when I use rasterio package saving function in the arguments there is some input argument not defined related to either extent, projected coordinates systems or maybe transform .... ?? – Gab May 19 '19 at 0:34
0

I noticed that you are not writing out the transform. I think that is the main problem.

Here is how you can retrieve and write the information:


with rasterio.open(filename) as src:
    dst_transform = src.transform
    dst_crs = src.crs

##  saving bands individualy  
new_dataset = rasterio.open(
       'ReflectanceB1.tif',
       'w',
       driver='GTiff',
       height=band_blue_reflectance.shape[0],
       width=band_blue_reflectance.shape[1],
       count=1,
       dtype=band_blue_reflectance.dtype,
       crs=dst_crs,
       transform=dst_transform,
)   

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.