Topics

Python rasterio for saveing GeoTIFF files and read in ArcGIS or QGIS


Gabriel Cotlier
 

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
 


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.


Sean Gillies
 

Thank you for providing a complete example of your code and the ArcMap error message. Without those details, it can be hard to provide help.

When you open the files for writing, there is one thing missing: the *transform* keyword argument that determines the spatial extent of the dataset. I think that in your case, you will have success if you copy the value from the source datasets:

filename = "20180308_133037_1024_3B_AnalyticMS.tif"

# Copy the transform and crs objects from the source dataset.
with rasterio.open(filename) as src:
    dst_transform = src.transform
    dst_crs = src.crs

# Pass the copied values to rasterio.open().
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
)    
    
new_dataset.write(band_blue_reflectance, 1)
new_dataset.close()


On Sat, May 18, 2019 at 12:39 PM <gabiklm01@...> wrote:
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
 


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.



--
Sean Gillies


Gabriel Cotlier
 

Dear Sean,
Thank you very much for your help.
I appreciate it.
Regards
Gabriel


On Monday, May 20, 2019, Sean Gillies <sean.gillies@...> wrote:
Thank you for providing a complete example of your code and the ArcMap error message. Without those details, it can be hard to provide help.

When you open the files for writing, there is one thing missing: the *transform* keyword argument that determines the spatial extent of the dataset. I think that in your case, you will have success if you copy the value from the source datasets:

filename = "20180308_133037_1024_3B_AnalyticMS.tif"

# Copy the transform and crs objects from the source dataset.
with rasterio.open(filename) as src:
    dst_transform = src.transform
    dst_crs = src.crs

# Pass the copied values to rasterio.open().
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
)    
    
new_dataset.write(band_blue_reflectance, 1)
new_dataset.close()


On Sat, May 18, 2019 at 12:39 PM <gabiklm01@...> wrote:
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
 


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.



--
Sean Gillies