File: hdfwrite.pro

package info (click to toggle)
coyote 2019.01.29-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster
  • size: 6,316 kB
  • sloc: python: 184; makefile: 14; sh: 13
file content (146 lines) | stat: -rw-r--r-- 3,823 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
PRO HDFWRITE, filename

  ; Create randomly-distributed data.
  
seed = -1L
x = RANDOMU(seed, 40)
y = RANDOMU(seed, 40)
distribution = SHIFT(DIST(40,40), 25, 15)
distribution = EXP(-(distribution/15)^2)

lat = x * (24./1.0) + 24
lon = y * 50.0/1.0 - 122
temp = distribution(x*40, y*40) * 273

   ; Select the name of a new file and open it.

IF N_ELEMENTS(filename) EQ 0 THEN filename = Pickfile(/Write)
IF filename EQ '' THEN RETURN
PRINT, ''
PRINT, 'Opening HDF file "' + filename + '" ...'

   ; If there is a problem opening the file, catch the error.
   
CATCH, error
IF error NE 0 THEN BEGIN
   PRINT, ''
   PRINT, 'Unable to obtain an SDS file ID.'
   PRINT, 'This file may already be open by an HDF routine.'
   PRINT, 'Returning...'
   PRINT, ''
   RETURN
ENDIF
fileID = HDF_SD_START(filename, /CREATE)
CATCH, /CANCEL

   ; Write some attributes into the file.

PRINT, 'Writing file attributes ...'   
version = 'MacOS 4.0.1b'
date = "Jan 1, 1997"
experiment = "Experiment 25A36M"
name = 'David Fanning'
email = 'david@idlcoyote.com'

HDF_SD_ATTRSET, fileID, 'VERSION', version
HDF_SD_ATTRSET, fileID, 'DATE', date
HDF_SD_ATTRSET, fileID, 'EXPERIMENT', experiment
HDF_SD_ATTRSET, fileID, 'NAME', name
HDF_SD_ATTRSET, fileID, 'EMAIL ADDRESS', email

   ; Create the SDS data sets for the raw data. 
   
PRINT, 'Writing raw data ...'
latsdsID = HDF_SD_CREATE(fileID, "Raw Latitude", [40L], /FLOAT)
lonsdsID = HDF_SD_CREATE(fileID, "Raw Longitude", [40L], /FLOAT)
tempsdsID = HDF_SD_CREATE(fileID, "Raw Temperature", [40L], /FLOAT)

   ; Write the raw data.
   
HDF_SD_ADDDATA, latsdsID, lat
HDF_SD_ADDDATA, lonsdsID, lon
HDF_SD_ADDDATA, tempsdsID, temp

   ; Terminate access to the raw data SDSs.

HDF_SD_ENDACCESS, latsdsID
HDF_SD_ENDACCESS, lonsdsID
HDF_SD_ENDACCESS, tempsdsID

   ; Grid the irregularly spaced, raw data.
   
latMax = 49.0
latMin = 24.0
lonMax = -67.0
lonMin = -125.0
mapBounds = [lonMin, latMin, lonMax, latMax]
mapSpacing = [0.5, 0.25]
 
TRIANGULATE, lon, lat, FVALUE=temp, SPHERE=triangles, /DEGREES
gridData = TRIGRID(temp, SPHERE=triangles, /DEGREES, /EXTRAPOLATE, $
  mapSpacing, mapBounds)

   ; Calculate xlon and ylat vectors corresponding to gridded data.
   
s = SIZE(gridData)
gridlon = FINDGEN(s(1))*((lonMax - lonMin)/(s(1)-1)) + lonMin
gridlat = FINDGEN(s(2))*((latMax - latMin)/(s(2)-1)) + latMin

   ; Now store the gridded data.

PRINT, 'Writing gridded data ...'
gridID = HDF_SD_CREATE(fileID, "Gridded Data", [s(1), s(2)], /FLOAT)
HDF_SD_ADDDATA, gridID, gridData

   ; Store local SDS attributes with the gridded data.
   
method = 'Delauney Triangulation'
routines = 'TRIANGULATE, TRIGRID'

PRINT, 'Writing gridded data set attributes ...'
HDF_SD_AttrSet, gridID, 'METHOD', method
HDF_SD_AttrSet, gridID, 'IDL ROUTINES', routines
HDF_SD_AttrSet, gridID, 'GRID SPACING', mapSpacing
HDF_SD_AttrSet, gridID, 'MAP BOUNDRIES', mapBounds

   ; Store the scales for the two dimensions of the gridded data.

PRINT, 'Writing dimension data ...'
longitudeDimID = HDF_SD_DIMGETID(gridID, 0)
HDF_SD_DIMSET, longitudeDimID, $
   LABEL='Longitude', $
   NAME='Longitude Dimension', $
   SCALE=gridlon, $
   UNIT='Degrees'
   
latitudeDimID = HDF_SD_DIMGETID(gridID, 1)
HDF_SD_DIMSET, latitudeDimID, $
   LABEL='Latitude', $
   NAME='Latitude Dimension', $
   SCALE=gridlat, $
   UNIT='Degrees'
   
   ; Load a color palette you will use to display the data.
   
thisDevice = !D.NAME
SET_PLOT, 'Z'
LOADCT, 5, /SILENT
TVLCT, r, g, b, /GET
palette = BYTARR(3,256)
palette(0,*) = r
palette(1,*) = g
palette(2,*) = b
SET_PLOT, thisDevice

   ; Store the color palette along with the data.

PRINT, 'Writing color palette ...'
HDF_DFP_ADDPAL, filename, palette

   ; Close everything up properly.
   
HDF_SD_ENDACCESS, gridID
HDF_SD_END, fileID
PRINT, 'Write Operation Completed'
PRINT, ''
END