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:

8 comments:

  1. Hi, I was wondering what version of Python you're using?

    ReplyDelete
    Replies
    1. Hi, I am using Python 2.7, but this post here is purely gdal command line without python (installation I used http://geoinformaticstutorial.blogspot.no/search/label/Installation )

      Delete
    2. I copied the file "nt_197810_n07_v01_n2.bin" onto my desktop and changed the icefile assignment in the code accordingly. Both pieces of code ran in Spyder without any errors, but the resulting .tif file on my desktop is solid black. Any idea what the problem may be?

      Delete
    3. It only contains 0s? Cant say just now, but maybe the 1978 files are somewhat different? Can you try with the same file I used just to see if that is working?

      Delete
    4. I think it was just a simple mistake I was making, but I was able to convert and open the files in ArcGIS. Thanks so much!

      Does the code also take care of projecting the data or is that something I will still have to do?

      Delete
    5. The data is in EPSG:3411 but you will have to apply the projection info. How to do this, look at the end of the Bin2GeoTiff function here https://github.com/npolar/RemoteSensing/blob/master/IcePersistency/IcePersistency.py

      Delete
    6. Would I still use the code you wrote above on this page, or would I just use the entire function Bin2GeoTiff by itself in Python?

      Delete
    7. You use lines 103-113 (lines may change at a later update) but you of course have to adjust these commands -- you would need some basic knowledge of Python and gdal. But in these lines you see the projection information, which you of course can apply in ArcGIS or QGIS more simply if you only have a single file

      Delete