Thursday, February 27, 2014

Reading binary data (NSIDC Sea Ice Concentrations) into GeoTIFF raster with Python

I need to read NSIDC sea ice concentration charts from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data into a raster file. These files are according to the documentation "scaled, unsigned flat binary with one byte per pixel, and therefore have no byte order, or endianness. Data are stored as one-byte integers representing scaled sea ice concentration values.[...] The file format consists of a 300-byte descriptive header followed by a two-dimensional array of one-byte values containing the data."

The blog entry, which helped med finding out how to do this is this one but it had to be modified for my files.


Part One uses the Python struct command to unpack he binary file into a tuple
import struct, numpy, gdal

#Dimensions from https://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
height = 448
width = 304

#for this code, inspiration found at https://stevendkay.wordpress.com/category/python/
icefile = open(r"U:\SSMI\nt_20120101_f17_v01_n.bin", "rb")
contents = icefile.read()
icefile.close()

# unpack binary data into a flat tuple z
s="%dB" % (int(width*height),)
z=struct.unpack_from(s, contents, offset = 300)

nsidc = numpy.array(z).reshape((448,304))
nsidc = numpy.rot90(nsidc, 1)
The "unpack_from" considers the offset of the header bytes and the tuple "z" contains now all the values in one long row of values -- this one-dimensional string of values is converted into an array and the reshaped according to the dimensions of the raster. Now I have my values in a two-dimensional array, the next part of the code writes it into a GeoTIFF file:
#write the data to a Geotiff
driver = gdal.GetDriverByName("GTiff")
outraster = driver.Create('C:\Users\max\Desktop\\test4.tif', height, width, 1, gdal.GDT_Int16 )
raster = numpy.zeros((width, height), numpy.float) 
outraster.GetRasterBand(1).WriteArray( raster )

outband = outraster.GetRasterBand(1)

#Write to file     
outband.WriteArray(nsidc)
outband.FlushCache()

#Clear arrays and close files
outband = None
iceraster = None
outraster = None
outarray = None

This then is the resulting and imported image:

Friday, February 21, 2014

Reading ERS / ASAR Data and Creating Quicklooks

In this previous post it was described how gdal reads Radarsat SAR data. I want to read ERS /ASAR data and create quicklooks for each file.

ERS /ASAR data has different formats, files ending with *.001 , *.E1, *.E2 or *. N1  Gdal can read all these formats, also the *.001 files which not all software reads. For the latter, the "DAT_01.001" is the file to be passed to gdal.

As first step, I create a list containing all files with these endings. I use os.walk to recursively search through all folders:

filelist_CEOS = []
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, 'DAT_01.001'):
      filelist_CEOS.append(os.path.join(root, filename))

filelist_E1E2 = []
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, '*.E1'):
      filelist_E1E2.append(os.path.join(root, filename))      
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, '*.E2'):
      filelist_E1E2.append(os.path.join(root, filename))  
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, '*.N1'):
      filelist_E1E2.append(os.path.join(root, filename)) 

gdalwarp is able to read all of this formats and I can proceed to process all the files from this filelist, calling gdalwarp from within a Python script. This map-projects the raw SAR file and creates a GeoTIFF:


os.system('gdalwarp -tps  -t_srs EPSG:32633 ' + file_from_filelist + ' ' + outputfilename )  


For a small quicklook a convert it to jpeg and make it smaller:
os.system('gdal_translate -of JPEG -ot byte -outsize 20% 20% -scale 0 1000 0 255 ' + outputfilename + ' ' + browseimage )

The whole script reading in all ERS /ASAR images and creating a quicklook can be found here and is hopefully documented well enough.