04 Work with different data formats

This is part of Python for Geosciences notes.

================

Binary data

Open binary

Data from Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data.

In [ ]:
!wget ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/north/monthly/nt_200709_f17_v01_n.bin

Create file id:

In [14]:
ice = fromfile('nt_200709_f17_v01_n.bin', dtype='uint8')

We use uint8 data type. List of numpy data types

The file format consists of a 300-byte descriptive header followed by a two-dimensional array.

In [15]:
ice = ice[300:]

Reshape

In [16]:
ice = ice.reshape(448,304)

Simple visualisation of array with imshow (Matplotlib function):

In [17]:
imshow(ice)
colorbar()
Out[17]:

To convert to the fractional parameter range of 0.0 to 1.0, divide the scaled data in the file by 250.

In [18]:
ice = ice/250.
imshow(ice)
colorbar()
Out[18]:

Let's mask all land and missing values:

In [19]:
ice_masked = ma.masked_greater(ice, 1.0)
imshow(ice_masked)
colorbar()
Out[19]:

Masking in this case is similar to using NaN in Matlab. More about NumPy masked arrays

Save binary

In [20]:
fid = open('My_ice_2007.bin', 'wb')
ice.tofile(fid)
fid.close()

In order to work with other data formats we need to use one of the SciPy submodules:

SciPy

General purpose scientific library (that consist of bunch of sublibraries) and builds on NumPy arrays.

We are going to use only scipy.io library.

scipy.io

Open .mat files

First we have to load function that works with Matlab files:

In [21]:
from scipy.io import loadmat

We are going to download Polar science center Hydrographic Climatology (PHC) for January in Matlab format.

In [ ]:
!wget https://www.dropbox.com/s/0kuzvz03gw6d393/PHC_jan.mat

Open file:

In [22]:
all_variables = loadmat('PHC_jan.mat')

We can look at the names of variables stored in the file:

In [23]:
all_variables.keys()
Out[23]:
['PTEMP1', 'LON', '__header__', '__globals__', 'DEPTH', 'LAT', '__version__']

We need only PTEMP1 (3d potential temperature).

In [24]:
temp = numpy.array(all_variables['PTEMP1'])

Check variable's shape:

In [25]:
temp.shape
Out[25]:
(33, 180, 360)

Show surface level:

In [26]:
imshow(temp[0,:,:])
colorbar()
Out[26]:

Open netCDF files

Import nessesary function:

In [29]:
from scipy.io import netcdf

I am going to download NCEP reanalysis data. Surface 4 daily air temperature for 2012.

In [27]:
!wget ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface/air.sig995.2012.nc

#Alternative for the times of US goverment shutdowns:
#!wget http://database.rish.kyoto-u.ac.jp/arch/ncep/data/ncep.reanalysis/surface/air.sig995.2012.nc

--2013-10-27 00:19:39--  ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface/air.sig995.2012.nc
           => ‘air.sig995.2012.nc’
Resolving ftp.cdc.noaa.gov (ftp.cdc.noaa.gov)... 140.172.38.117
Connecting to ftp.cdc.noaa.gov (ftp.cdc.noaa.gov)|140.172.38.117|:21... connected.
Logging in as anonymous ... Logged in!
==> SYST ... done.    ==> PWD ... done.
==> TYPE I ... done.  ==> CWD (1) /Datasets/ncep.reanalysis/surface ... done.
==> SIZE air.sig995.2012.nc ... 30793412
==> PASV ... done.    ==> RETR air.sig995.2012.nc ... done.
Length: 30793412 (29M) (unauthoritative)

100%[======================================>] 30.793.412  1,22MB/s   in 27s    

2013-10-27 00:20:09 (1,10 MB/s) - ‘air.sig995.2012.nc’ saved [30793412]

Create file id:

In [30]:
fnc = netcdf.netcdf_file('air.sig995.2012.nc')

It's not really file id, it's netcdf_file object, that have some methods and attributes:

In [32]:
fnc.description
Out[32]:
'Data is from NMC initialized reanalysis\n(4x/day).  These are the 0.9950 sigma level values.'
In [33]:
fnc.history
Out[33]:
'created 2011/12 by Hoop (netCDF2.3)'

list variables

In [34]:
fnc.variables
Out[34]:
{'air': ,
 'lat': ,
 'lon': ,
 'time': }

Access information about variables

In [35]:
air = fnc.variables['air']

This time we create netcdf_variable object, that contain among other things attributes of the netCDF variable as well as data themselves.

In [37]:
air.actual_range
Out[37]:
array([ 191.1000061,  323.       ], dtype=float32)
In [38]:
air.long_name
Out[38]:
'4xDaily Air temperature at sigma level 995'
In [39]:
air.units
Out[39]:
'degK'
In [40]:
print(air.add_offset)
print(air.scale_factor)

512.81
0.01
In [41]:
air.shape
Out[41]:
(1464, 73, 144)

We can access the data by simply using array syntax. Here we show first time step of our data set:

In [44]:
imshow(air[0,:,:])
colorbar()
Out[44]:

Data are packed. In order to unpack we have to use scale_factor and add_offset values from attributes of the netCDF air variable. We also convert to $^{\circ}$C:

In [45]:
air_c = ((air[:] * air.scale_factor) + air.add_offset) - 273.15

Check the values of our new variable:

In [46]:
imshow(air_c[0,:,:])
colorbar()
Out[46]:

Save netCDF file

Minimalistic variant :)

In [ ]:
!rm test_netcdf.nc
fw = netcdf.netcdf_file('test_netcdf.nc', 'w')

fw.createDimension('t', 1464)
fw.createDimension('y', 73)
fw.createDimension('x', 144)

air_var = fw.createVariable( 'air','float32', ('t', 'y', 'x'))
air_var[:] = air_c[:]
fw.close()

More descriptive variant:

In [ ]:
!rm test_netcdf.nc
fw = netcdf.netcdf_file('test_netcdf.nc', 'w')

fw.createDimension('TIME', 1464)
fw.createDimension('LATITUDE', 73)
fw.createDimension('LONGITUDE', 144)

time = fw.createVariable('TIME', 'f', ('TIME',))
time[:] = fnc.variables['time'][:]
time.units = 'hours since 1-1-1 00:00:0.0' 

lat  = fw.createVariable('LATITUDE', 'f', ('LATITUDE',))
lat[:] = fnc.variables['lat'][:]

lon = fw.createVariable('LONGITUDE', 'f', ('LONGITUDE',))
lon[:] = fnc.variables['lon'][:]

ha = fw.createVariable('New_air','f', ('TIME', 'LATITUDE', 'LONGITUDE'))
ha[:] = air_c[:]
ha.missing_value = -9999.

fw.close()

Comments !

links

social