Создание netCDF файла из кучи бинарных при помощи NCL
Задача: взять данные, которые распространяются в виде жуткой кучи бинарников и сделать из неё один красивый, легко обрабатываемый файл netCDF
?нструмент: NCL
В незапамятные времена, когда интернет у всех был медленным и дисковое пространство было на вес золота, королём форматов для распространения данных был, так называемый, бинарный формат. Это обычная последовательность нулей и единиц, которая не содержит в себе никакой информации о данных находящихся внутри. Чтобы «расшифровать» его вам нужно добыть информацию о структуре файла, о порядке байтов (endianness), о типе данных (float, integer), о том не умножены ли данные на какое-нибудь значение, что за единицы измерения используются, что за координаты, за какую дату и время эти данные и так далее и тому подобное. Геморрой. Единственное достоинство бинарного формата это размер получающихся файлов, он сравнительно мал.
Прогресс не стоял на месте и были изобретены, так называемые, self-describing форматы, одним из которых является netCDF. Преимущество этих форматов в том, что вам не надо пытать создателей файлов для того, чтобы получить информацию о его содержании, в идеале вся необходимая информация должна уже содержаться в самом файле. Это позволяет создавать приложения однообразно работающие с различными данными представленными в self-describing форматах, а не заново изобретать велосипед для каждого нового бинарника.
Размер файлов netCDF больше чем бинарных, но если это было очень важным в 90е, то в 2000х большинству нет разницы качать один или два мегабайта, телефонные модемы постепенно вымирают и слава богу. Тем не менее, некоторые несознательные граждане продолжают предоставлять данные в бинарном формате, то-ли по привычке, то-ли ленясь изменять схемы обработки, работающие с 90х (работает, не трогай). Ещё одной причиной может быть то, что создание правильного netCDF файла это не самая простая задача (в отличии от его чтения).
Далее я приведу код для довольно простого создания netCDF файла из бинарных при помощи языка программирования для наук о Земле NCL о котором я рассказывал в одном из предыдущих постов.
Обрабатывать будем данные по концентрации морского льда в Северном Ледовитом океане, описание которых расположено тут http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
Файлы, с которыми мы будем работать, это бинарники, тип данных byte, размер 304 columns x 448 rows. Это не совсем клинический случай бинарников, поскольку некое описание у них всё таки имеется в первых 300 байтах. Нам оно не нужно, поэтому придётся его вырезать. Сами значения нужно поделить на 250 для того, чтобы получить концентрацию льда (от нуля до единицы). Всё что больше единицы — missing values.
Сетка, на которой сидят наши данные, нерегулярная, поэтому нужны массивы с широтами и долготами для каждой точки. Они также упрятаны в бинарники, которые можно скачать отсюда ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/ Файлы psn25lats_v2.dat и psn25lons_v2.dat. Значения хранятся в 4-byte integers (little endian) и ко всему умножены на 100000 (то есть нам надо поделить 🙂 ). Размерность, естественно, такая же как у самих данных 304 columns x 448 rows.
Вот эту кучу нам предстоит разгрести и создать из неё нормальный netCDF файл.
Диспозиция должна быть следующей. В папке с .ncl скриптом должны находиться файлы psn25lats_v2.dat и psn25lons_v2.dat . ?сходные бинарники должны лежать по отношению к папке со скриптом в ./data/year . Где year — папка в которой лежат файлы за определённый год.
Скачивать файлы по одному слава богу нет необходимости, поскольку сервер NSIDC поддерживает tar-on-the-fly и схема закачки, например 1998 года будет следующей
(login as anonymous, and use your e-mail address as password)
ftp>cd /pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/north/daily/
ftp>get 1998.tar
Большое спасибо Маше Цукерник из NCAR, которая предоставила мне основу для данного скрипта. Вот собственно и он
;Data description http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
;Data location ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/
;Created by Kolya Koldunov koldunovn@gmail.com, IMPRS-ESM, http://kodunov.net
;Based on script by Maria (Masha) Tsukernik, NCAR
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;name for lat nad lon files
latfilename = "psn25lats_v2.dat"
lonfilename = "psn25lons_v2.dat"
;open lat lon files
lat_raw = cbinread(latfilename,-1,"integer")
lon_raw = cbinread(lonfilename,-1,"integer")
lat = lat_raw/100000.
lon = lon_raw/100000.
lat2d = onedtond(lat,(/448,304/))
lon2d = onedtond(lon,(/448,304/))
lat2d@units = "degrees_north"
lon2d@units = "degrees_east"
lat2d@long_name = "latitude"
lat2d@lstandard_name = "grid_latitude"
lon2d@long_name = "longitude"
lon2d@standard_name = "grid_longitude" ;
yyyy = 1998
dirin="./data/"+yyyy+"/"
if (isleapyear(yyyy)) then ntimes=366
else ntimes=365
end if
fili=systemfunc("cd "+dirin+"; ls nt_"+yyyy+"*.bin")
nfil = dimsizes(fili)
;print("Number of files in yr "+yyyy+" : "+nfil+", number of actual timesteps: "+ntimes)
dimi = dimsizes(lat2d)
nlat = dimi(0)
nlon = dimi(1)
ice3d = new((/ntimes,nlat,nlon/),typeof(lat2d))
x = fspan(60,90,448)
y = fspan(-180,180,304)
time = new((/ntimes/),"double")
x@units = "degrees_north"
y@units = "degrees_east"
tunits= "days since 1978-01-01"
time@units=tunits
x!0 = "x"
y!0 = "y"
time!0 = "time"
x&x = x
y&y = y
time&time = time
lat2d!0 = "x"
lat2d!1 = "y"
lon2d!0 = "x"
lon2d!1 = "y"
lat2d&x = x
lat2d&y = y
lon2d&x = x
lon2d&y = y
yyyymmdd=new((/ntimes/),"integer")
yyyymmdd!0="time"
yyyymmdd@long_name="yyyymmdd: year, month, day"
do j=0,ntimes-1
;==============================
;j = 120
day=j+1
mmdd=monthday(yyyy,day)
mon=mmdd/100
dd=mod(mmdd,100)
;yyyymmdd(j)=yyyy+10000+mon*100+dd
yyyymmdd(j)=yyyy*10000+mmdd
time(j)=ut_inv_calendar(yyyy,mon,dd,0,0,0,time@units,0)
stringday=(sprinti("%0.4i",mmdd))
everyfil=dirin+"nt_"+yyyy+stringday+"_f13_v01_n.bin"
; if (isfilepresent(everyfil)) then
print("opening file: "+everyfil)
;open the data file to inn variable
inn = cbinread(everyfil,-1,"byte")
;convert our data from byte to float and drop first 300 records (which is the header we don't need)
ice_raw = byte2flt(inn(300:))
;to get a real values you have to divide values by 250
ice_raw= ice_raw/250
;in out data set we have several values that are not the data
;251/250 = 1.004 Circular mask used in the Arctic to cover the irregularly-shaped data gap around the pole (caused by the orbit inclination and instrument swath)
;252 Unused
;253/250 = 1.012 Coastlines
;254/250 = 1.016 Superimposed land mask
;255/250 = 1.02 Missing data
ice_raw04 = ice_raw
ice_raw04@_FillValue = 1.004
ice_raw12 = ice_raw
ice_raw12@_FillValue = 1.012
ice_raw16 = ice_raw
ice_raw16@_FillValue = 1.016
ice_raw02 = ice_raw
ice_raw02@_FillValue = 1.02
; here is a strange part
; we summ up all variables to "collect" _FillValue from 4 different variables. As result we have a variable with one _FillValue for everything that is not tsea ice concentration. Then we divide it to return to original values.
ice_raw_sum = ice_raw12+ice_raw16+ice_raw04+ice_raw02
ice_fill = ice_raw_sum/4
;reshape our 1D arrays to 2D arrays
ice2d = onedtond(ice_fill,(/448,304/))
;set up atributes for lat and lon files, and then make this files atributes of our data (ice2d). This is done to be able to plot data on curveliniat coordinates. If you plan laiter to put data to netCDF, you better remove this attributes from ice2d.
;ice2d@lat2d = lat2d
;ice2d@lon2d = lon2d
ice3d(j,:,:) = ice2d
end do
ice3d!0 = "time"
ice3d!1 = "x"
ice3d!2 = "y"
ice3d&time = time
ice3d&x = x
ice3d&y = y
ice3d@long_name = "ice compactness"
ice3d@units = "frac."
ice3d@code = 15
ice3d@coordinates = "lon2d lat2d"
ice3d@missinf_value = ice3d@_FillValue
ntim = dimsizes(time) ; get dimension sizes
nx = dimsizes(x)
ny = dimsizes(y)
delete_VarAtts(ice2d,(/"lon2d","lat2d"/))
system("/bin/rm -f out.nc") ; remove if exists
fout = addfile ("out.nc", "c") ; open output file
setfileoption(fout,"DefineMode",True)
fAtt = True ; assign file attributes
fAtt@title = "Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I Passive Microwave Data"
fAtt@source_file = "bynary files from ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/"
fAtt@Conventions = "None"
fAtt@creation_date = systemfunc ("date")
fileattdef( fout, fAtt ) ; copy file attributes
dimNames = (/"time", "x", "y"/)
dimSizes = (/ -1 , nx, ny /)
dimUnlim = (/ True , False, False/)
filedimdef(fout,dimNames,dimSizes,dimUnlim)
filevardef(fout, "time" ,typeof(time),getvardims(time))
filevardef(fout, "x" ,typeof(x),getvardims(x) )
filevardef(fout, "y" ,typeof(y),getvardims(y))
filevardef(fout, "lon2d" ,typeof(lon2d),getvardims(lon2d))
filevardef(fout, "lat2d" ,typeof(lat2d) ,getvardims(lat2d))
filevardef(fout, "ICE" ,typeof(ice3d) ,getvardims(ice3d))
filevarattdef(fout,"time" ,time) ; copy time attributes
;filevarattdef(fout,"x" ,x) ; copy lev attributes
;filevarattdef(fout,"y" ,y) ; copy lat attributes
filevarattdef(fout,"lon2d" ,lon2d) ; copy lon attributes
filevarattdef(fout,"lat2d" ,lat2d) ; copy PS attributes
filevarattdef(fout,"ICE",ice3d) ; copy TOPOG attributes
setfileoption(fout,"DefineMode",False)
;
fout->time = (/time/)
;fout->x = (/x/)
;fout->y = (/y/)
fout->lon2d = (/lon2d/)
fout->lat2d = (/lat2d/)
fout->ICE = (/ice3d/)
Построчно разбирать скрипт у меня нет ни сил ни желания 😉 Если будут какие-то непонятные моменты, спрашивайте в комментариях, попробую пояснить. В принципе важные вещи откомментированы и не должны вызвать больших затруднений, плюс у NCL отличная документация, к которой всегда можно обратиться и выяснить, что делает та или иная функция.
В итоге, мы получаем файл out.nc содержащий в себе ежедневную концентрацию морского льда в Арктике за 1998 год. Метаданные в нём минимальные, их можно посмотреть сделав ncdump -h out.nc
dimensions:
time = UNLIMITED ; // (365 currently)
x = 448 ;
y = 304 ;
variables:
double time(time) ;
time:units = "days since 1978-01-01" ;
time:_FillValue = -9999. ;
float x(x) ;
float y(y) ;
float lon2d(x, y) ;
lon2d:standard_name = "grid_longitude" ;
lon2d:long_name = "longitude" ;
lon2d:units = "degrees_east" ;
float lat2d(x, y) ;
lat2d:lstandard_name = "grid_latitude" ;
lat2d:long_name = "latitude" ;
lat2d:units = "degrees_north" ;
float ICE(time, x, y) ;
ICE:missinf_value = 1.012f ;
ICE:coordinates = "lon2d lat2d" ;
ICE:code = 15 ;
ICE:units = "frac." ;
ICE:long_name = "ice compactness" ;
ICE:_FillValue = 1.012f ;
// global attributes:
:creation_date = "Сбт Фев 7 18:05:27 CET 2009" ;
:Conventions = "None" ;
:source_file = "bynary files from ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/" ;
:title = "Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I Passive Microwave Data" ;
}
Но, этой информации достаточно для того, чтобы просматривать файл в самом распространённом просмотрщике netCDF — ncview.
?ли посчитать среднемесячные значения при помощи cdo
Надеюсь, что вы сможете модифицировать этот скрипт для ваших данных и сделать свой вклад в изживание бинарников со свету и приумножение netCDF файлов! )