将 hdf5 gadget2 格式文件中的parttype属性从1更改为X（例如2）

问题描述 投票：0回答：1

我正在尝试适应我的需求的脚本：https://github.com/eltevo/StePS/blob/master/tools/Snapshot_Format_Converters/ascii2hdf5_snapshot.py

此脚本从 ascii 文件中读取所有粒子（该文件是由同一存储库中将 hdf5 转换为 ascii 的其他脚本创建的）并将它们放入 hdf5 文件中。

我需要通过粒子的质量符号来区分粒子。如果为正，则进入 PartType1。如果为负数，则进入 PartType2（此处不存在）。这样我就可以使用 gadgetviewer 或 glnemo（或其他）在同一视图中可视化这两个组。

遗憾的是，我对 python 的了解不够，无法将这样的功能添加到这个脚本中。


#*******************************************************************************#
#  ascii2hdf5_snapshot.py - An ASCII to hdf5 file converter for StePS           #
#     (STEreographically Projected cosmological Simulations) snapshots.         #
#    Copyright (C) 2017-2022 Gabor Racz                                         #
#                                                                               #
#    This program is free software; you can redistribute it and/or modify       #
#    it under the terms of the GNU General Public License as published by       #
#    the Free Software Foundation; either version 2 of the License, or          #
#    (at your option) any later version.                                        #
#                                                                               #
#    This program is distributed in the hope that it will be useful,            #
#    but WITHOUT ANY WARRANTY; without even the implied warranty of             #
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              #
#    GNU General Public License for more details.                               #
#*******************************************************************************#

import numpy as np
import h5py
import sys
import time
# %matplotlib inline

#Beginning of the script
if len(sys.argv) != 4:
    print("Error:")
    print("usage: ./ascii2hdf5_snapshot.py <input ASCII snapshot> <output HDF5 snapshot> <precision 0: 32bit 1: 64bit>\nExiting.")
    sys.exit(2)

if int(sys.argv[3]) != 0 and int(sys.argv[3]) != 1:
    print("Error:")
    print("Unkown output precision.\nExiting.")
    sys.exit(2)

print("Reading the input ASCII file...")
start = time.time()
ASCII_snapshot=np.fromfile(str(sys.argv[1]), count=-1, sep='\t', dtype=np.float64)
ASCII_snapshot = ASCII_snapshot.reshape(int(len(ASCII_snapshot)/7),7)
N=len(ASCII_snapshot)
M_min = np.min(ASCII_snapshot[:,6])
R_max = np.max(np.abs(ASCII_snapshot[:,0:3]))
print("Number of particles:\t%i\nMinimal mass:\t%f*10e11M_sol\nMaximal radius:\t%fMpc" % (N,M_min,R_max))
end = time.time()
print("..done in %fs. \n\n" % (end-start))

print("Saving the snapshot in HDF5 format...")
start = time.time()
HDF5_snapshot = h5py.File(str(sys.argv[2]), "w")
#Creating the header
header_group = HDF5_snapshot.create_group("/Header")
#Writing the header attributes
header_group.attrs['NumPart_ThisFile'] = np.array([0,N,0,0,0,0],dtype=np.uint32)
header_group.attrs['NumPart_Total'] = np.array([0,N,0,0,0,0],dtype=np.uint32)
header_group.attrs['NumPart_Total_HighWord'] = np.array([0,0,0,0,0,0],dtype=np.uint32)
header_group.attrs['MassTable'] = np.array([0,0,0,0,0,0],dtype=np.float64)
header_group.attrs['Time'] = np.double(1.0)
header_group.attrs['Redshift'] = np.double(0.0)
header_group.attrs['Lbox'] = np.double(R_max*2.01)
header_group.attrs['NumFilesPerSnapshot'] = int(1)
header_group.attrs['Omega0'] = np.double(1.0)
header_group.attrs['OmegaLambda'] = np.double(0.0)
header_group.attrs['HubbleParam'] = np.double(1.0)
header_group.attrs['Flag_Sfr'] = int(0)
header_group.attrs['Flag_Cooling'] = int(0)
header_group.attrs['Flag_StellarAge'] = int(0)
header_group.attrs['Flag_Metals'] = int(0)
header_group.attrs['Flag_Feedback'] = int(0)
header_group.attrs['Flag_Entropy_ICs'] = int(0)
#Header created.
#Creating datasets for the particle data
particle_group = HDF5_snapshot.create_group("/PartType1")
if int(sys.argv[3]) == 0:
    HDF5datatype = 'float32'
    npdatatype = np.float32
if int(sys.argv[3]) == 1:
    HDF5datatype = 'double'
    npdatatype = np.float64
X = particle_group.create_dataset("Coordinates", (N,3),dtype=HDF5datatype)
V = particle_group.create_dataset("Velocities", (N,3),dtype=HDF5datatype)
IDs = particle_group.create_dataset("ParticleIDs", (N,),dtype='uint64')
M = particle_group.create_dataset("Masses", (N,),dtype=HDF5datatype)
#Saving the particle data
X[:,:] = ASCII_snapshot[:,0:3]
V[:,:] = ASCII_snapshot[:,3:6]
M[:] = ASCII_snapshot[:,6]
IDs[:] = np.arange(N, dtype=np.uint64)
HDF5_snapshot.close()
end = time.time()
print("..done in %fs. \n\n" % (end-start))

我询问了聊天 GPT 但我们无法找到解决方案呵呵。

python-3.x hdf5
1个回答
0
投票

我更新了这样的代码：

#!/usr/bin/env python3

#*******************************************************************************#
#  ascii2hdf5_snapshot.py - An ASCII to hdf5 file converter for StePS           #
#     (STEreographically Projected cosmological Simulations) snapshots.         #
#    Copyright (C) 2017-2022 Gabor Racz                                         #
#                                                                               #
#    This program is free software; you can redistribute it and/or modify       #
#    it under the terms of the GNU General Public License as published by       #
#    the Free Software Foundation; either version 2 of the License, or          #
#    (at your option) any later version.                                        #
#                                                                               #
#    This program is distributed in the hope that it will be useful,            #
#    but WITHOUT ANY WARRANTY; without even the implied warranty of             #
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              #
#    GNU General Public License for more details.                               #
#*******************************************************************************#

import numpy as np
import h5py
import sys
import time
# %matplotlib inline

#Beginning of the script
if len(sys.argv) != 4:
    print("Error:")
    print("usage: ./ascii2hdf5_snapshot.py <input ASCII snapshot> <output HDF5 snapshot> <precision 0: 32bit 1: 64bit>\nExiting.")
    sys.exit(2)

if int(sys.argv[3]) != 0 and int(sys.argv[3]) != 1:
    print("Error:")
    print("Unkown output precision.\nExiting.")
    sys.exit(2)

print("Reading the input ASCII file...")
start = time.time()


ASCII_snapshot=np.fromfile(str(sys.argv[1]), count=-1, sep='\t', dtype=np.float64)
ASCII_snapshot = ASCII_snapshot.reshape(int(len(ASCII_snapshot)/7),7)

N=len(ASCII_snapshot)

posmass=np.where(ASCII_snapshot[:, 6, None] > 0, ASCII_snapshot, 0)
negmass=np.where(ASCII_snapshot[:, 6, None] < 0, ASCII_snapshot, 0)
posmass=posmass[~np.all(posmass==0, axis=1)]
negmass=negmass[~np.all(negmass==0, axis=1)]
N_pos=len(posmass)
N_neg=len(negmass)
M_min = np.min(negmass[:,6])
R_max = np.max(np.abs(ASCII_snapshot[:,0:3]))

print("Number of particles:\t%i\nNumber of positive masses:\t%i\nNumber of Negative Masses\t%i\nMinimal mass:\t%f*10e11M_sol\nMaximal radius:\t%fMpc" % (N,N_neg,N_pos,M_min,R_max))

end = time.time()
print("..done in %fs. \n\n" % (end-start))

print("Saving the snapshot in HDF5 format...")
start = time.time()
HDF5_snapshot = h5py.File(str(sys.argv[2]), "w")
#Creating the header
header_group = HDF5_snapshot.create_group("/Header")
#Writing the header attributes
header_group.attrs['NumPart_ThisFile'] = np.array([0,N,0,0,0,0],dtype=np.uint32)
header_group.attrs['NumPart_Total'] = np.array([0,N,0,0,0,0],dtype=np.uint32)
header_group.attrs['NumPart_Total_HighWord'] = np.array([0,0,0,0,0,0],dtype=np.uint32)
header_group.attrs['MassTable'] = np.array([0,M_min,0,0,0,0],dtype=np.float64)
header_group.attrs['Time'] = np.double(1.0)
header_group.attrs['Redshift'] = np.double(0.0)
header_group.attrs['Lbox'] = np.double(R_max*2.01)
header_group.attrs['NumFilesPerSnapshot'] = int(1)
header_group.attrs['Omega0'] = np.double(1.0)
header_group.attrs['OmegaLambda'] = np.double(0.0)
header_group.attrs['HubbleParam'] = np.double(1.0)
header_group.attrs['Flag_Sfr'] = int(0)
header_group.attrs['Flag_Cooling'] = int(0)
header_group.attrs['Flag_StellarAge'] = int(0)
header_group.attrs['Flag_Metals'] = int(0)
header_group.attrs['Flag_Feedback'] = int(0)
header_group.attrs['Flag_Entropy_ICs'] = int(0)
#Header created.
#Creating datasets for the POSITIVE particle data
particle_group = HDF5_snapshot.create_group("/PartType1")
if int(sys.argv[3]) == 0:
    HDF5datatype = 'float32'
    npdatatype = np.float32
if int(sys.argv[3]) == 1:
    HDF5datatype = 'double'
    npdatatype = np.float64
X = particle_group.create_dataset("Coordinates", (N_pos,3),dtype=HDF5datatype)
V = particle_group.create_dataset("Velocities", (N_pos,3),dtype=HDF5datatype)
IDs = particle_group.create_dataset("ParticleIDs", (N_pos,),dtype='uint64')
M = particle_group.create_dataset("Masses", (N_pos,),dtype=HDF5datatype)
#Saving the particle data
X[:,:] = posmass[:,0:3]
V[:,:] = posmass[:,3:6]
M[:] = posmass[:,6]
IDs[:] = np.arange(N_pos, dtype=np.uint64)
###################################################################################
#Creating datasets for the NEGATIVE particle data
particle_groupN = HDF5_snapshot.create_group("/PartType2")
if int(sys.argv[3]) == 0:
    HDF5datatype = 'float32'
    npdatatype = np.float32
if int(sys.argv[3]) == 1:
    HDF5datatype = 'double'
    npdatatype = np.float64
Xn = particle_groupN.create_dataset("Coordinates", (N_neg,3),dtype=HDF5datatype)
Vn = particle_groupN.create_dataset("Velocities", (N_neg,3),dtype=HDF5datatype)
IDsn = particle_groupN.create_dataset("ParticleIDs", (N_neg,),dtype='uint64')
Mn = particle_groupN.create_dataset("Masses", (N_neg,),dtype=HDF5datatype)
#Saving the particle data
Xn[:,:] = negmass[:,0:3]
Vn[:,:] = negmass[:,3:6]
Mn[:] = negmass[:,6]
IDsn[:] = np.arange(N_neg, dtype=np.uint64)
HDF5_snapshot.close()
end = time.time()
print("..done in %fs. \n\n" % (end-start))

这段代码执行得很好。也在合理的时间内。

但是当我尝试在 gadgetviewer 中查看结果时，它无法正常工作，它崩溃了。在 GLnemo2 中工作也很奇怪。我得到 3D 绘图，但没有组。由于某种原因，Gadgetviewer 拒绝绘制任何不是parttype1 组的内容。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.