-
-
Save jgomezdans/5488682 to your computer and use it in GitHub Desktop.
| #!/usr/bin/env python | |
| """ | |
| SYNOPSIS | |
| dn_2_rad.py [-h,--help] [-v,--verbose] [-i,--input] [-o,--output] | |
| DESCRIPTION | |
| This program is used to extract the gain parameters and to convert | |
| Landsat TM "digital numbers" (DNs) to top-of-atmosphere (TOA) radiances. | |
| The units of the output are W/(m2*ster*um). We assume that the user | |
| provides an input band and that the band names follow the standard USGS | |
| naming convention. | |
| EXAMPLES | |
| $ ./dn_2_rad.py --input LE71660672004161ASN01_B1.TIF --verbose | |
| Tue Apr 30 14:26:06 2013 | |
| Start reading LE71660672004161ASN01_B1.TIF | |
| Using blocksize 7831 x 1 | |
| Creating output LE71660672004161ASN01_B1_TOARAD.tif | |
| Start reading LE71660672004161ASN01_B2.TIF | |
| Using blocksize 7831 x 1 | |
| Creating output LE71660672004161ASN01_B2_TOARAD.tif | |
| Start reading LE71660672004161ASN01_B3.TIF | |
| Using blocksize 7831 x 1 | |
| Creating output LE71660672004161ASN01_B3_TOARAD.tif | |
| Start reading LE71660672004161ASN01_B4.TIF | |
| Using blocksize 7831 x 1 | |
| Creating output LE71660672004161ASN01_B4_TOARAD.tif | |
| Start reading LE71660672004161ASN01_B5.TIF | |
| Using blocksize 7831 x 1 | |
| Creating output LE71660672004161ASN01_B5_TOARAD.tif | |
| Start reading LE71660672004161ASN01_B7.TIF | |
| Using blocksize 7831 x 1 | |
| Creating output LE71660672004161ASN01_B7_TOARAD.tif | |
| Tue Apr 30 14:26:31 2013 | |
| TOTAL TIME IN MINUTES: 0.429338383675 | |
| EXIT STATUS | |
| -1 if numpy and/or GDAL aren't present | |
| AUTHOR | |
| J Gomez-Dans <j.gomez-dans@ucl.ac.uk> | |
| LICENSE | |
| This script is in the public domain, free from copyrights or restrictions. | |
| """ | |
| import sys | |
| import os | |
| import traceback | |
| import optparse | |
| import time | |
| try: | |
| import numpy as np | |
| except ImportError: | |
| print "You need numpy installed" | |
| sys.exit ( -1 ) | |
| try: | |
| from osgeo import gdal | |
| except ImportError: | |
| print "You need GDAL installed" | |
| sys.exit ( -1 ) | |
| GDAL_OPTS = ["COMPRESS=LZW", "INTERLEAVE=PIXEL", "TILED=YES",\ | |
| "SPARSE_OK=TRUE", "BIGTIFF=YES" ] | |
| def process_metadata ( fname ): | |
| """A function to extract the relelvant metadata from the | |
| USGS control file. Returns dicionaries with LMAX, LMIN, | |
| QCAL_LMIN and QCAL_LMAX for each of the bands of interest.""" | |
| fp = open( fname, 'r') # Open metadata file | |
| lmax = {} # Dicts to store constants | |
| lmin = {} | |
| qc_lmax = {} | |
| qc_lmin = {} | |
| gain = {} | |
| bias = {} | |
| for line in fp: # | |
| # Check for LMAX and LMIN strings | |
| # Note that parse logic is identical to the first case | |
| # This version of the code works, but is rather inelegant! | |
| if ( line.find ("RADIANCE_MULT_BAND") >= 0 ): | |
| s = line.split("=") # Split by equal sign | |
| the_band = int(s[0].split("_")[3]) # Band number as integer | |
| if the_band in [1,2,3,4,5,7]: # Is this one of the bands we want? | |
| gain[the_band] = float ( s[-1] ) # Get constant as float | |
| elif ( line.find ("RADIANCE_ADD_BAND") >= 0 ): | |
| s = line.split("=") # Split by equal sign | |
| the_band = int(s[0].split("_")[3]) # Band number as integer | |
| if the_band in [1,2,3,4,5,7]: # Is this one of the bands we want? | |
| bias[the_band] = float ( s[-1] ) # Get constant as float | |
| elif ( line.find ("QUANTIZE_CAL_MAX_BAND") >= 0 ): | |
| s = line.split("=") # Split by equal sign | |
| the_band = int(s[0].split("_")[4]) # Band number as integer | |
| if the_band in [1,2,3,4,5,7]: # Is this one of the bands we want? | |
| qc_lmax[the_band] = float ( s[-1] ) # Get constant as float | |
| elif ( line.find ("QUANTIZE_CAL_MIN_BAND") >= 0 ): | |
| s = line.split("=") # Split by equal sign | |
| the_band = int(s[0].split("_")[4]) # Band number as integer | |
| if the_band in [1,2,3,4,5,7]: # Is this one of the bands we want? | |
| qc_lmin[the_band] = float ( s[-1] ) # Get constant as float | |
| elif ( line.find ("RADIANCE_MAXIMUM_BAND") >= 0 ): | |
| s = line.split("=") # Split by equal sign | |
| the_band = int(s[0].split("_")[3]) # Band number as integer | |
| if the_band in [1,2,3,4,5,7]: # Is this one of the bands we want? | |
| lmax[the_band] = float ( s[-1] ) # Get constant as float | |
| elif ( line.find ("RADIANCE_MINIMUM_BAND") >= 0 ): | |
| s = line.split("=") # Split by equal sign | |
| the_band = int(s[0].split("_")[3]) # Band number as integer | |
| if the_band in [1,2,3,4,5,7]: # Is this one of the bands we want? | |
| lmin[the_band] = float ( s[-1] ) # Get constant as float | |
| return ( bias, gain, lmax, lmin, qc_lmax, qc_lmin ) | |
| def get_metadata ( fname ): | |
| """This function takes `fname`, a filename (opionally with a path), and | |
| and works out the associated metadata file""" | |
| original_fname = os.path.basename ( fname ) | |
| metadata_fname = original_fname.split("_")[0] + "_MTL.txt" | |
| metadata_fname = os.path.join ( os.path.dirname ( fname ), metadata_fname ) | |
| return metadata_fname | |
| def extract_chunk ( the_file, n_blocks=1 ): | |
| """A function that extracts a chunk from a datafile""" | |
| ds_config = {} | |
| g = gdal.Open ( the_file ) | |
| block_size = g.GetRasterBand(1).GetBlockSize() | |
| nx = g.RasterXSize | |
| ny = g.RasterYSize | |
| the_bands = np.arange ( g.RasterCount ) + 1 | |
| proj = g.GetProjectionRef() | |
| geoT = g.GetGeoTransform() | |
| ds_config['nx'] = nx | |
| ds_config['ny'] = ny | |
| ds_config['nb'] = g.RasterCount | |
| ds_config['geoT'] = geoT | |
| ds_config['proj'] = proj | |
| block_size = [ block_size[0]*n_blocks, block_size[1]*n_blocks ] | |
| # store these numbers in variables that may change later | |
| nx_valid = block_size[0] | |
| ny_valid = block_size[1] | |
| # find total x and y blocks to be read | |
| nx_blocks = (int)((nx + block_size[0] - 1) / block_size[0]); | |
| ny_blocks = (int)((ny + block_size[1] - 1) / block_size[1]); | |
| buf_size = block_size[0]*block_size[1] | |
| print("Using blocksize %s x %s" %(block_size[0], block_size[1])) | |
| sys.stdout.flush() | |
| ################################################################ | |
| # start looping through blocks of data | |
| ################################################################ | |
| # loop through X-lines | |
| for X in xrange( nx_blocks ): | |
| # change the block size of the final piece | |
| if X == nx_blocks-1: | |
| nx_valid = nx - X * block_size[0] | |
| buf_size = nx_valid*ny_valid | |
| # find X offset | |
| this_X = X*block_size[0] | |
| # reset buffer size for start of Y loop | |
| ny_valid = block_size[1] | |
| buf_size = nx_valid*ny_valid | |
| # loop through Y lines | |
| for Y in xrange( ny_blocks ): | |
| # change the block size of the final piece | |
| if Y == ny_blocks-1: | |
| ny_valid = ny - Y * block_size[1] | |
| buf_size = nx_valid*ny_valid | |
| # find Y offset | |
| this_Y = Y*block_size[1] | |
| buf = g.ReadRaster(this_X, this_Y, nx_valid, ny_valid, \ | |
| buf_xsize=nx_valid, buf_ysize=ny_valid, \ | |
| band_list= the_bands ) | |
| a = np.frombuffer(buf, dtype=np.uint8 ) | |
| if len(the_bands) > 1: | |
| data_in = a.reshape(( len(the_bands), ny_valid, nx_valid)) | |
| else: | |
| data_in = a.reshape(( ny_valid, nx_valid)) | |
| yield ( ds_config, this_X, this_Y, nx_valid, ny_valid, data_in ) | |
| def convert_to_radiance ( fname, band, gain, bias, lmax, lmin, qc_lmax, qc_lmin, verbose ): | |
| """ | |
| This function reads the input file chunk by chunk using ``extract_chunk``, | |
| and converts from DN to radiance using the methodology given by | |
| <http://landsat.usgs.gov/how_is_radiance_calculated.php>. | |
| """ | |
| # This variable is used to create the output dataset when the first | |
| # chunk of data is read. | |
| first_time = True | |
| # `extract_chunk` can be used to efficiently read chunks of data | |
| # from a GDAL dataset. It uses a "generator", which means that we | |
| # can iterate over it using e.g. a for loop | |
| if verbose: | |
| print "Start reading %s" % fname | |
| # The output filaname. | |
| output_fname = fname.replace(".TIF", "_TOARAD.tif" ) | |
| for ( ds_config, this_X, this_Y, nx_valid, ny_valid, data_in ) in \ | |
| extract_chunk ( fname ): | |
| if first_time: | |
| if verbose: | |
| print "Creating output %s" % output_fname | |
| # Create output dataset if `first_time` is true | |
| drv = gdal.GetDriverByName ( "GTiff" ) | |
| dst_ds = drv.Create ( output_fname, ds_config['nx'], \ | |
| ds_config['ny'], 1, gdal.GDT_Float32, options=GDAL_OPTS ) | |
| dst_ds.SetGeoTransform ( ds_config['geoT'] ) | |
| dst_ds.SetProjection ( ds_config['proj'] ) | |
| first_time = False | |
| # Create the output array. Not needed, but we can | |
| # conveniently set the output type here to float32 | |
| radiance = np.zeros_like ( data_in, dtype=np.float32 ) | |
| # DN to radiance conversion if we have a sensible DN | |
| passer = np.logical_and ( qc_lmin < data_in, data_in < qc_lmax ) | |
| radiance = np.where ( passer, \ | |
| (( lmax - lmin )/(qc_lmax-qc_lmin))*\ | |
| ( (data_in*gain + bias) - qc_lmin) + lmin, \ | |
| -999 ) | |
| # Write the output chunk | |
| dst_ds.GetRasterBand ( 1 ).WriteArray( radiance, \ | |
| xoff=this_X, yoff=this_Y) | |
| # Add some useful metadata to the dataset | |
| dst_ds.GetRasterBand( 1 ).SetNoDataValue ( -999 ) | |
| dst_ds.GetRasterBand( 1 ).SetMetadata ( {"Band": "%d" % band, \ | |
| "Units": "W/(m2*ster*um)", \ | |
| "Data": "TOA Radiance" } ) | |
| # Need to do this to flush the dataset to disk | |
| dst_ds = None | |
| def main (): | |
| """The main function""" | |
| global options | |
| global args | |
| metadata_file = get_metadata ( options.input_f ) | |
| ( gain, bias, lmax, lmin, qc_lmax, qc_lmin ) = process_metadata ( metadata_file ) | |
| original_fname = os.path.basename ( options.input_f) | |
| prefix = original_fname.split("_")[0] | |
| fname_prefix = os.path.join ( os.path.dirname ( options.input_f ), prefix ) | |
| for the_band in [1,2,3,4,5,7]: | |
| input_file = fname_prefix + "_B%d.TIF" % the_band | |
| convert_to_radiance ( input_file, the_band, \ | |
| gain[the_band], bias[the_band], \ | |
| lmax[the_band], lmin[the_band], qc_lmax[the_band], \ | |
| qc_lmin[the_band], options.verbose ) | |
| if __name__ == '__main__': | |
| start_time = time.time() | |
| parser = optparse.OptionParser(formatter=optparse.TitledHelpFormatter(), \ | |
| usage=globals()['__doc__'] ) | |
| parser.add_option ('-v', '--verbose', action='store_true', \ | |
| default=False, help='verbose output') | |
| parser.add_option ('-i', '--input', action='store', dest="input_f",\ | |
| type=str, help='Input filename') | |
| (options, args) = parser.parse_args() | |
| if options.verbose: print time.asctime() | |
| main() | |
| if options.verbose: print time.asctime() | |
| if options.verbose: print 'TOTAL TIME IN MINUTES:', | |
| if options.verbose: print (time.time() - start_time) / 60.0 |
Why is it necessary to put the fist_time flag? I tried to do some similar to your code but also implementing reflectance, DOS1 and DOS2, but only after the first calculations are made and image saved, the metadata file still exist, then is deleted i.e. I load metadata, then calculate radiance for band 1 and saved new image, then proceed to reflectance but at that moment the metadafile is delete from system, I debugged coded and find that in line
[227] dst_ds.SetGeoTransform ( ds_config['geoT'] ) or
[228]: dst_ds.SetProjection ( ds_config['proj'] )
is when metadata file is delete, do you know what is causing that?
Thanks.
from gdalconst import *
RED= gdal.Open("d:)
NIR= gdal.Open("d:)
print ('size ',RED.RasterXSize,'x',RED.RasterYSize,
'x',RED.RasterCount)
print ('proek ',RED.GetProjection())
geotransform = RED.GetGeoTransform()
if not geotransform is None:
print ('nachal cor (',geotransform[0], ',',geotransform[3],')')
print ('razmer pix = (',geotransform[1], ',',geotransform[5],')')
band_RED = RED.GetRasterBand(1)
band_NIR = NIR.GetRasterBand(1)
proj = RED.GetProjection()
print ('Тип данных',gdal.GetDataTypeName(band_RED.DataType))
print ('Тип данных',gdal.GetDataTypeName(band_NIR.DataType))
if band_RED.GetOverviewCount() > 0:
print ('Канал содержит ', band_RED.GetOverviewCount(),' обзорных изображений.')
if not band_RED.GetRasterColorTable() is None:
print ('Канал содержит таблицу цветов с ',
band_RED.GetRasterColorTable().GetCount(), ' записями.')
if band_NIR.GetOverviewCount() > 0:
print ('Канал содержит ', band_NIR.GetOverviewCount(),' обзорных изображений.')
if not band_NIR.GetRasterColorTable() is None:
print ('Канал содержит таблицу цветов с ',
band_NIR.GetRasterColorTable().GetCount(), ' записями.')
array_RED = band_RED.ReadAsArray()
cut_RED = array_RED[3000:4000,3000:4000]
geotransform_cut = RED.GetGeoTransform()
print (cut_RED)
array_NIR = band_NIR.ReadAsArray()
cut_NIR = array_NIR[3000:4000,3000:4000]
geotransform_cut = NIR.GetGeoTransform()
print (cut_NIR)
calib_RED = (cut_RED *mtl['REFLECTANCE_MULT_BAND_4'])+mtl['REFLECTANCE_ADD_BAND_4']
calib_NIR = (cut_NIR *mtl['REFLECTANCE_MULT_BAND_5'])+mtl['REFLECTANCE_ADD_BAND_5']
NDVI = (calib_NIR - calib_RED)/(calib_NIR + calib_RED)
#plt.imshow(NDVI)
#plt.show()
def Save (image, name):
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create("d:\Rul\Python\Landsat 8{}.TIF".format(name), 1000, 1000, 1, gdal. GDT_Float32)
dst_ds.SetGeoTransform([geotransform[0]+3000geotransform[1],geotransform[1], 0,
geotransform[3]+3000geotransform[5], 0, geotransform[5]])
dst_ds.SetProjection(proj)
dst_ds.GetRasterBand(1).WriteArray(image)
Save (calib_RED , 'RED')
Save (calib_NIR , 'NIR')
Save (NDVI , 'NDVI')
Is it possible you've switched around gain/bias?
https://gist.github.com/jgomezdans/5488682#file-dn_2_rad-py-L126
return ( bias, gain, lmax, lmin, qc_lmax, qc_lmin )
https://gist.github.com/jgomezdans/5488682#file-dn_2_rad-py-L260
( gain, bias, lmax, lmin, qc_lmax, qc_lmin ) = process_metadata ( metadata_file )
Thanks!
Is it possible you've switched around gain/bias?
https://gist.github.com/jgomezdans/5488682#file-dn_2_rad-py-L126
return ( bias, gain, lmax, lmin, qc_lmax, qc_lmin )https://gist.github.com/jgomezdans/5488682#file-dn_2_rad-py-L260
( gain, bias, lmax, lmin, qc_lmax, qc_lmin ) = process_metadata ( metadata_file )Thanks!
I think the order should be followed.
def test_return():
a = 1
b = 2
return (a , b)
and then
(b, a) = test_retun()
print('a: %s' % a)
print('b: %s' % b)
I think this piece of code inludes error:
"radiance = np.where ( passer,
(( lmax - lmin )/(qc_lmax-qc_lmin))_
( (data_in_gain + bias) - qc_lmin) + lmin,
-999 )"
I don't uderstand this "(data_in_gain + bias)",
I think that (( lmax - lmin )/(qc_lmax-qc_lmin))_( data_in - qc_lmin) + lmin is more true
article about calculation of radiance and reflectance
http://gis-lab.info/docs/revised_landsat-5_tm_radiometric_calibration_procedures_and_postcalibration_dynamic_ranges.pdf
I would like to use your code for my project