CODEWITHSUNDEEP
pythongrads
iyui
iyui

Reputation: 31

How to write in GrADS-readable binary format using python?

I have been trying to export data stored as numpy array to GrADS-flat binary. It seems GrADS does not recognize Z dimension given in the .ctl file. Whatever value I use for 'set z integer', GrADS only shows the first level.

Here is the minimal reproduction of my problem

my python code:

import numpy as np
from array import array

data = np.linspace(1, 60, num=60, endpoint=True)
data = np.reshape(data, [5, 4, 3])
print(data)

with open('temp.dat', 'ab') as wf:
  float_array = array('f', data.flatten())
  float_array.tofile(wf)
wf.close()

executing this writes numbers [1, 2, 3, ..., 60] as single-precision float to a binary file

my .ctl file:

DSET  ^temp.dat
TITLE title
UNDEF -9.99E33
XDEF 3 LINEAR 0.0  1 
YDEF 4 LINEAR 0.0  1 
ZDEF 5 LEVELS 0 1 2 3 4
TDEF 1 LINEAR 0Z10apr1991 12hr 
VARS      1
var     0  99  some var
ENDVARS

This set of .dat and .ctl files shows first 12 numbers as first level of the field as expected.

ga-> open temp.ctl 
Scanning description file:  temp.ctl
Data file temp.dat is open as file 1
LON set to 0 2 
LAT set to 0 3 
LEV set to 0 0 
Time values set: 1991:4:10:0 1991:4:10:0 
E set to 1 1 
ga-> set digsize 0.6
digsiz = 0.6 
ga-> set lon -1 3
LON set to -1 3 
ga-> set lat -1 4
LAT set to -1 4 
ga-> set gxout grid
ga-> d var

enter image description here

However, if I try to 'set z 2' it still shows the first level. Orelse, var(z=3) and var(z=1) being

var(z=1)=
[[ 1.  2.  3.]
  [ 4.  5.  6.]
  [ 7.  8.  9.]
  [10. 11. 12.]]
  
var(z=3)=
 [[25. 26. 27.]
  [28. 29. 30.]
  [31. 32. 33.]
  [34. 35. 36.]]

this should show the constant field of "24" ... but GrADS is showing "0" as if z=3 is the same with z=1.

ga-> c
ga-> d var(z=3) - var(z=1)

enter image description here

What is even more ununderstantable is that if I give var2 in .ctl file, grads recognize what it should be var(z=2) as var2(z=1)!

I know there are lots of better visualization tools other than GrADS, but I want to use legacy GrADS codes so it is unavoidable.
Did I write the binary in the wrong order? or is the binary file missing header or separator or something? I am grad if anyone knows why it is happening. Thanks in advance.

Upvotes: 0

Views: 578

Answers (1)

iyui
iyui

Reputation: 31

My colleague has just pointed out that the problem is in .ctl file. I have set the layer number to zero, which GrADS recognize as "surface layer", which is a single layer variable field and can be overlayed upon any layer above it.

I am able to show any layer of the field with the corrected .ctl file.

my corrected .ctl.

DSET  ^temp.dat
TITLE title
UNDEF -9.99E33
XDEF 3 LINEAR 0.0  1 
YDEF 4 LINEAR 0.0  1 
ZDEF 5 LEVELS 0 1 2 3 4
TDEF 1 LINEAR 0Z10apr1991 12hr 
VARS      1
// this should be 5 not zero because there are five layers!
// var     0  99  some var
var     5  99  some var
ENDVARS

Upvotes: 0

Related Questions