; Converts the currently loaded data to a format that Amira can read ; ; Author: John Wise ( December 2007 ) ; ; Modified: tweeked the format a bit (TA) ;------------------------------------------------------------------------ pro hdf5_write, dataset_id, name, data, ATTRIBUTE=ATTRIBUTE, DATASET=DATASET if (NOT KEYWORD_SET(ATTRIBUTE) AND KEYWORD_SET(DATASET)) then $ DATASET = 1 ATTRIBUTE = KEYWORD_SET(ATTRIBUTE) DATASET = KEYWORD_SET(DATASET) datatype_id = H5T_IDL_CREATE(data) ndim = size(data, /DIMENSIONS) if (ndim gt 0) then $ dataspace_id = H5S_CREATE_SIMPLE(size( data, /DIMENSIONS )) $ else $ dataspace_id = H5S_CREATE_SCALAR() if (KEYWORD_SET(ATTRIBUTE)) then begin a_id = H5A_CREATE(dataset_id, name, datatype_id, dataspace_id) H5A_WRITE, a_id, data H5A_CLOSE, a_id endif if (KEYWORD_SET(DATASET)) then begin d_id = H5D_CREATE(dataset_id, name, datatype_id, dataspace_id) H5D_WRITE, d_id, data H5D_CLOSE, d_id endif H5S_CLOSE, dataspace_id H5T_CLOSE, datatype_id end function return_filename_parts, full_path parts = STRSPLIT(full_path, ":::", /EXTRACT) last_slash = STRPOS(parts[0], "/", /REVERSE_SEARCH)+1 filename = STRMID(parts[0], last_slash) return, [filename, parts[1]] end ;------------------------------------------------------------------------ pro construct_amira_export, only_current=ocurrent @common_blocks.inc if N_elements(ocurrent) eq 0 then begin ocurrent = 0 these_grids = all_grid_info endif else these_grids=grid_info total_ngrids = N_ELEMENTS(these_grids) test_grid = these_grids[total_ngrids-1] ;; Check what kind of file we're reading data_format = determine_data_format_from_file_name(test_grid.baryon_file) use_references = (data_format eq 5 or data_format eq 6) if (not use_references) then begin allsize = TOTAL( PRODUCT(these_grids.dim, 1) ) warning_message = $ ["Original data are HDF4. Data duplication is necessary for export to Amira.",$ "Total exported data = " + $ STRING(4*allsize/1048576, format='(G0.3)') + " MB. Continue?"] result = DIALOG_MESSAGE(warning_message, /QUESTION) if (result eq "No") then return endif if (data_format eq 4 or data_format eq 5 and $ MAX(these_grids.level) le 0) then $ these_grids = set_level_info(these_grids) use_times = times outputtimes = N_elements(use_times) if ocurrent then BEGIN outputtimes=1 these_grids.timestep = 1 these_grids.time = max(these_grids.time) times = (these_grids.time) times = times[UNIQ(times)] endif timesteps = these_grids.timestep timesteps = timesteps[UNIQ(timesteps)] field_name = list_str[var_index] default_filename = "Movie" + field_name + "_" + $ STRING(timesteps[0], timesteps[N_ELEMENTS(timesteps)-1], $ format='(I0,"-",I0)') + ".a5" filename = DIALOG_PICKFILE(FILTER="*.h5;*.a5", PATH=data_dir, $ FILE=default_filename, DIALOG_PARENT=BASE2, $ TITLE="Select Amira Export Filename", $ GET_PATH=data_dir, /NOCONFIRM, /WRITE, $ /OVERWRITE_PROMPT) if (filename eq "") then return spawn, 'rm '+filename fid = H5F_CREATE(filename) ; Global metadata rg1 = these_grids[(where(these_grids.level eq 0))[0]] root_dx = (rg1.Right_edge - rg1.Left_edge) / rg1.dim stag = 1L if (ocurrent and data_format eq 6) then stag = 0L gid = H5G_CREATE(fid, "globalMetaData") hdf5_write, gid, "datatype", 0L, /ATTRIBUTE hdf5_write, gid, "staggering", stag, /ATTRIBUTE hdf5_write, gid, "fieldtype", 1L, /ATTRIBUTE hdf5_write, gid, "rootDelta", root_dx, /ATTRIBUTE hdf5_write, gid, "minTimeStep", LONG(MIN(timesteps)), /ATTRIBUTE ; older versions of amira need this: hdf5_write, gid, "maxTimeStep", LONG(MAX(timesteps)), /ATTRIBUTE hdf5_write, gid, "minTime", MIN(times), /ATTRIBUTE hdf5_write, gid, "maxTime", MAX(times), /ATTRIBUTE hdf5_write, gid, "numTimeSteps", LONG(N_ELEMENTS(timesteps)), /ATTRIBUTE H5G_CLOSE, gid ; Sorted times and integer timesteps sorted_times = times[SORT(times)] sorted_timesteps = timesteps[SORT(timesteps)] hdf5_write, fid, "sorted_times", sorted_times, /DATASET hdf5_write, fid, "sorted_timesteps", sorted_timesteps, /DATASET ; Grids sorted by times, then levels if N_elements(these_grids) gt 1 then $ hh_times = HISTOGRAM(these_grids.timestep, REVERSE_INDICES=rtime) $ else begin hh_times = these_grids.timestep rtime = [0,1] endelse refine_by = [2L,2L,2L] ghostzoneFlags = LONARR(6) numGhostzones = LONARR(3) ;breaki for i = 0L, outputtimes-1 do begin if (rtime[i] eq rtime[i+1]) then continue if ocurrent then tgi = these_grids else $ tgi = these_grids[rtime[(rtime[i]):(rtime[i+1]-1)]] if (N_ELEMENTS(tgi) eq 1) then levels = [tgi.level] else levels = tgi.level hh_level = HISTOGRAM(levels, REVERSE_INDICES=rlevel, $ OMIN=min_level, OMAX=max_level) time_id = H5G_CREATE(fid, STRING(sorted_timesteps[i], FORMAT='("time-",I0)')) hdf5_write, time_id, "time", times[i], /ATTRIBUTE hdf5_write, time_id, "numLevels", (max_level-min_level+1L), /ATTRIBUTE print, sorted_timesteps[i], format='("Constructing amira file for timestep ",I0)' for level = min_level, max_level do begin bin = level-min_level dx = root_dx / refine_by^level ngrids = hh_level[bin] if (ngrids eq 0) then continue tlgi = tgi[rlevel[(rlevel[bin]):(rlevel[bin+1]-1)]] ;; Read data if HDF4 if (not use_references or ocurrent) then begin all_data = get_data(tlgi, field_name, const_sub) endif level_id = H5G_CREATE(time_id, STRING(level, FORMAT='("level-",I0)')) hdf5_write, level_id, "delta", dx, /ATTRIBUTE hdf5_write, level_id, "relativeRefinementFactor", refine_by, /ATTRIBUTE hdf5_write, level_id, "numGrids", ngrids, /ATTRIBUTE reverse_dims = 1 ;reverse_dims = 0 for grid = 0L, ngrids-1 do begin integer_origin = LONG64(tlgi[grid].left_edge / dx) if (use_references and not ocurrent) then begin grid_id = H5G_CREATE(level_id, STRING(grid, FORMAT='("grid-",I0)')) endif else begin ;; first data field = *all_data[grid].data[0] datatype_id = H5T_IDL_CREATE(*all_data[grid].data[0]) if reverse_dims then dataspace_id = H5S_CREATE_SIMPLE(reverse(all_data[grid].dims)) $ else dataspace_id = H5S_CREATE_SIMPLE((all_data[grid].dims)) ; grid_id = H5D_CREATE(level_id, STRING(grid, FORMAT='("grid-",I0)'), $ ; datatype_id, dataspace_id) grid_id = H5G_CREATE(level_id, STRING(grid, FORMAT='("grid-",I0)')) ; datatype_id, dataspace_id) ds_id = H5D_CREATE(grid_id, 'grid-data', $ datatype_id, dataspace_id) H5D_WRITE, ds_id, *all_data[grid].data[0] endelse hdf5_write, grid_id, "integerOrigin", integer_origin, /ATTRIBUTE hdf5_write, grid_id, "origin", tlgi[grid].left_edge, /ATTRIBUTE hdf5_write, grid_id, "ghostzoneFlags", ghostzoneFlags, /ATTRIBUTE hdf5_write, grid_id, "numGhostzones", numGhostzones, /ATTRIBUTE if reverse_dims then hdf5_write, grid_id, "dims", reverse(LONG(tlgi[grid].dim)), /ATTRIBUTE $ else hdf5_write, grid_id, "dims", (LONG(tlgi[grid].dim)), /ATTRIBUTE if (use_references and not ocurrent) then begin file_parts = return_filename_parts(tlgi[grid].baryon_file) file_parts[1] += field_name hdf5_write, grid_id, "referenceFileName", file_parts[0], /ATTRIBUTE hdf5_write, grid_id, "referenceDataPath", file_parts[1], /ATTRIBUTE endif ; older version: ; if (use_references) then H5G_CLOSE, grid_id else ; H5D_CLOSE, grid_id H5G_CLOSE, grid_id if (use_references) then H5D_CLOSE, ds_id endfor ;; grids H5G_CLOSE, level_id if (not use_references) then undefine, all_data endfor ;; level H5G_CLOSE, time_id heap_gc ; dirty but useful ... endfor ;; times H5F_CLOSE, fid end