python把一个二维数组读入到dat文件中
时间: 2024-05-05 18:19:35 浏览: 336
您可以使用Python的NumPy库将二维数组写入.dat文件中。以下是示例代码:
```python
import numpy as np
# 示例二维数组
arr = np.array([[1, 2], [3, 4]])
# 将二维数组写入.dat文件
with open('example.dat', 'wb') as f:
arr.tofile(f)
```
这将创建一个名为"example.dat"的文件,其中包含您的二维数组。要从.dat文件中读取数组,请使用以下代码:
```python
import numpy as np
# 从.dat文件读取二维数组
with open('example.dat', 'rb') as f:
arr = np.fromfile(f, dtype=np.int32).reshape((2, 2))
print(arr)
```
这将打印出:
```
[[1 2]
[3 4]]
```
请注意,您可能需要根据实际情况更改数据类型和数组维度。
相关问题
编写一个程序,将两个一维数组归并成一个由小到大有序的一维数组。将该数组写入二进制文件out.dat。 【输入形式】 第一行为一个整数n, 第二行为空格分隔的n个整数 第三行为另一个整数k, 第四行为空格分隔的k个整数,C
好的,这是一个归并排序的问题,我会帮您解决。下面是 Python 代码实现:
```python
def merge_sort(arr1, arr2):
i = j = 0
result = []
while i < len(arr1) and j < len(arr2):
if arr1[i] <= arr2[j]:
result.append(arr1[i])
i += 1
else:
result.append(arr2[j])
j += 1
result += arr1[i:]
result += arr2[j:]
return result
n = int(input())
arr1 = list(map(int, input().split()))
k = int(input())
arr2 = list(map(int, input().split()))
arr = merge_sort(arr1, arr2)
with open('out.dat', 'wb') as f:
for num in arr:
f.write(num.to_bytes(4, byteorder='little'))
```
这个程序首先定义了一个名为 `merge_sort` 的函数,该函数接受两个参数 `arr1` 和 `arr2`,分别表示两个待归并的数组。函数通过比较两个数组中的元素大小,并将它们按从小到大的顺序合并成一个新的数组 `result`。
接下来,程序读入输入数据,分别是 `n`、`arr1`、`k` 和 `arr2`。然后,程序调用 `merge_sort` 函数对两个数组进行归并排序,并将结果保存在名为 `arr` 的数组中。
最后,程序通过 `with` 语句打开一个名为 `out.dat` 的二进制文件,在其中写入归并排序后的结果。我们使用 `to_bytes` 方法将每个整数转换为 4 个字节的二进制数据,并将它们按照从小到大的顺序写入文件中。
希望这个程序能够帮助您解决问题!
! m=20: 纬向格点数 ! n=16: 经向格点数 ! d: 网格距 ! rm: 地图放大系数 ! f: 地转参数 ! w: 工作数组 ! cla,clo: 区域中心纬度和经度 ! dt: 时间步长 ! s: 平滑系数 ! ua,ub,uc:n-1,n,n+1时间层的纬向风速 ! va,vb,vc:n-1,n,n+1时间层的经向风速 ! za,zb,zc:n-1,n,n+1时间层的位势高度 ! na: 控制12小时预报的参数 ! nb: 记录时间积分步数的参数 ! nt2=72: 判别是否积分12小时,是否该做内点平滑; ! nt4=6: 判定是否该做边界平滑; ! nt5=36: 判定是否该做时间平滑。 ! zo: 为了减小重力惯性外波的波速,增加差分格式的稳定性而引入的位势高度。 ! ni: 是否进行初始风场的静力初始化。 ! ni=0为不进行初始化,使用读入的风场和高度场; ! ni=1为进行初始化,需要位势高度场做初值即可。 program Barotropic_Model parameter(m=20,n=16,d=300000.0,cla=51.0,clo=118.0,dt=600.0) dimension ua(m,n),va(m,n),za(m,n),ub(m,n),vb(m,n),zb(m,n),uc(m,n),vc(m,n),zc(m,n),rm(m,n),f(m,n),w(m,n) zo=2500.0 s=0.5 nt2=72 nt4=6 nt5=36 c1=dt/2.0 c2=dt*2.0 write(*,*) write(*,*)'!!!!欢迎使用正压原始方程模式!!!!' write(*,*) write(*,*)'打开文件、读入初始场、准备输出场文件...........' open(11,file='practice/Input/za.dat') ! 输入的高度场(文本) open(12,file='practice/Input/ua.dat') ! 输入的u风场(文本) open(13,file='practice/Input/va.dat') ! 输入的v风场(文本) open(17,file='practice/Input/za.grd',form='binary') ! 输出的高度场(二进制) open(18,file='practice/Input/ua.grd',form='binary') ! 输出的u风场(二进制) open(19,file='practice/Input/va.grd',form='binary') ! 输出的v风场(二进制) open(21,file='practice/Output/rm.dat') ! 地图放大系数 open(22,file='practice/Output/f.dat') ! 地转参数 open(23,file='practice/Output/ub.dat') ! 静力初始化得到的u风场(文本) open(24,file='practice/Output/ub.grd',form='binary') ! 静力初始化得到的u风场(二进制) open(25,file='practice/Output/vb.dat') ! 静力初始化得到的v风场(文本) open(26,file='practice/Output/vb.grd',form='binary') ! 静力初始化得到的v风场(二进制) open(27,file='practice/Output/zc.dat') ! 预报的高度场(文本) open(28,file='practice/Output/zc.grd',form='binary') ! 预报的高度场(二进制) open(29,file='practice/Output/uc.dat') ! 预报的u风场(文本) open(30,file='practice/Output/uc.grd',form='binary') ! 预报的u风场(二进制) open(31,file='practice/Output/vc.dat') ! 预报的v风场(文本) open(32,file='practice/Output/vc.grd',form='binary') ! 预报的v风场(二进制) ! 读入初始资料场 read(11,'(20f6.0)')za read(12,'(20f10.5)')ua read(13,'(20f10.5)')va write(*,*) write(*,*)'将初始高度场和风场写成二进制文件,便于Grads绘图......' write(17)((za(i,j),i=1,m),j=1,n) write(18)((ua(i,j),i=1,m),j=1,n) write(19)((va(i,j),i=1,m),j=1,n) ! 计算放大系数和地转参数,并写入数据文件中 write(*,*) write(*,*)'计算每个格点上的地图放大系数和地转参数,并写入对应输出文件......' call cmf(rm,f,d,cla,m,n) write(21,'(20f10.5)')rm write(22,'(20e15.5)')f write(*,*) write(*,*)'请输入数字0或1后按回车运行,0表示不进行静力初始化;1表示进行静力初始化。' write(*,*) write(*,*)'注意:如果求地转风的子程序(风场初始化)未完成,则只能输入数字0。' read(*,*)ni if(ni==1)then ! 计算地转风初值 write(*,*)'进行静力初始化,由风场求出高度场.....' call cgw(ua,va,za,rm,f,d,m,n) write(23,'(20f10.5)')ua write(24)((ua(i,j),i=1,m),j=1,n) write(25,'(20f10.5)')va write(26)((va(i,j),i=1,m),j=1,n) elseif(ni==0)then write(*,*)'不进行静力初始化使用给出的位势高度场和风场......' else write(*,*)'啊!!!输入了错误字符,请重新运行程序,输入0或1!' goto 90 endif ! 边值传送子程序 write(*,*) write(*,*)'固定边界条件赋值........' call tbv(ub,vb,zb,ua,va,za,m,n) call tbv(uc,vc,zc,ua,va,za,m,n) ! 开始预报 write(*,*) write(*,*)'开始12小时预报.......' do na=1,2 nb=0 ! 欧拉后差积分1小时 do nn=1,6 call ti(ua,va,za,ua,va,za,ub,vb,zb,rm,f,d,dt,zo,m,n) call ti(ua,va,za,ub,vb,zb,ua,va,za,rm,f,d,dt,zo,m,n) nb=nb+1 enddo ! 边界平滑子程序 call ssbp(za,w,s,m,n) call ssbp(ua,w,s,m,n) call ssbp(va,w,s,m,n) ! 前差积分半步 call ti(ua,va,za,ua,va,za,ub,vb,zb,rm,f,d,c1,zo,m,n) ! 中央差积分半步 call ti(ua,va,za,ub,vb,zb,uc,vc,zc,rm,f,d,dt,zo,m,n) nb=nb+1 ! 数组传送子程序 call ta(ub,vb,zb,uc,vc,zc,m,n) ! 中央差积分一步,共积分11小时 do nn=1,66 call ti(ua,va,za,ub,vb,zb,uc,vc,zc,rm,f,d,c2,zo,m,n) nb=nb+1 ! 打印积分步数,na大循环步,nb小循环步 call pv(na,nb) ! 判断是否积分12小时 if(nb.eq.nt2) go to 80 ! 判断是否做边界平滑 if(nb/nt4*nt4.eq.nb)then call ssbp(zc,w,s,m,n) call ssbp(uc,w,s,m,n) call ssbp(vc,w,s,m,n) else ! 判断是否做时间平滑 if(nb.eq.nt5 .or. nb.eq.nt5+1)then ! 时间平滑子程序 call ts(ua,ub,uc,va,vb,vc,za,zb,zc,s,m,n) else ! 数组传送,为下一轮积分做准备 call ta(ua,va,za,ub,vb,zb,m,n) call ta(ub,vb,zb,uc,vc,zc,m,n) endif endif enddo ! 区域内点平滑 80 call ssip(zc,w,s,m,n,k,2) call ssip(uc,w,s,m,n,k,2) call ssip(vc,w,s,m,n,k,2) ! 打印积分步数 call pv(na,nb) ! 数组传送,为后12小时的积分做准备 call ta(ua,va,za,uc,vc,zc,m,n) enddo ! 存放预报结果 write(*,*) write(*,*)'输出预报结果.......' write(27,'(20f6.0)') zc write(28) ((zc(i,j),i=1,m),j=1,n) write(29,'(20f10.5)')uc write(30) ((uc(i,j),i=1,m),j=1,n) write(31,'(20f10.5)')vc write(32) ((vc(i,j),i=1,m),j=1,n) write(*,*) write(*,*)'!!!预报结束!!!' 90 stop end ! computing map factors and coriolis parameter ! rk: 圆锥常数 ! rlq: 兰勃特投影映像平面上赤道到北极点的距离 ! a: 地球半径 ! sita:标准余纬 ! psx: 区域中心余纬 ! r: 模式中心到北极的距离 subroutine cmf(rm,f,d,cla,m,n) dimension rm(m,n),f(m,n) rk=0.7156 rlq=11423370.0 a=6371000.0 conv=57.29578 w1=2.0/rk sita=30.0/conv psx=(90.0-cla)/conv ! 计算模式中心到北极的距离r cel0=a*sin(sita)/rk cel=(tan(psx/2.0))/(tan(sita/2.0)) r=cel0*cel**rk ! 确定网格坐标原点在地图坐标系中的位置 xi0=-(m-1)/2.0 yj0=r/d+(n-1)/2.0 ! 求各格点至北极点的距离rl,(xj,yi)为模式各格点在地图坐标系中的位置 do i=1,m do j=1,n xi=xi0+(i-1) yj=yj0-(j-1) rl=sqrt(xi2+yj2)*d ! 求放大系数rm和柯氏参数f w2=(rl/rlq)**w1 sinl=(1.0-w2)/(1.0+w2) rm(i,j)=rk*rl/(a*sqrt(1.0-sinl**2)) f(i,j)=1.4584e-4*sinl enddo enddo return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! computing geostrophic winds ! 请同学编写地转风初值的子程序!!!应用书中(4.134)式 subroutine cgw(ua,va,za,rm,f,d,m,n) dimension ua(m,n),va(m,n),za(m,n),f(m,n),rm(m,n) end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! time integrations subroutine ti(ua,va,za,ub,vb,zb,uc,vc,zc,rm,f,d,dt,zo,m,n) dimension ua(m,n),va(m,n),za(m,n),ub(m,n),vb(m,n),zb(m,n),uc(m,n),vc(m,n),zc(m,n),rm(m,n),f(m,n) c=0.25/d m1=m-1 n1=n-1 do i=2,m1 do j=2,n1 e=-c*rm(i,j)*((ub(i+1,j)+ub(i,j))*(ub(i+1,j)-ub(i,j))+(ub(i,j)+ub(i-1,j))*(ub(i,j)-ub(i-1,j)) & +(vb(I,j-1)+vb(i,j))*(ub(i,j)-ub(i,j-1))+(vb(I,j)+vb(i,j+1))*(ub(i,j+1)-ub(i,j)) & +19.6*(zb(i+1,j)-zb(i-1,j)))+f(i,j)*vb(i,j) uc(i,j)=ua(i,j)+e*dt g=-c*rm(i,j)*((ub(I+1,j)+ub(i,j))*(vb(i+1,j)-vb(i,j))+(ub(I,j)+ub(i-1,j))*(vb(i,j)-vb(i-1,j)) & +(vb(I,j-1)+vb(i,j))*(vb(i,j)-vb(i,j-1))+(vb(I,j)+vb(i,j+1))*(vb(i,j+1)-vb(i,j)) & +19.6*(zb(i,j+1)-zb(i,j-1)))-f(i,j)*ub(i,j) vc(i,j)=va(i,j)+g*dt enddo enddo do i=2,m1 do j=2,n1 h=-c*rm(i,j)*((ub(I+1,j)+ub(i,j))*(zb(i+1,j)-zb(i,j))+(ub(I,j)+ub(i-1,j))*(zb(i,j)-zb(i-1,j)) & +(vb(I,j-1)+vb(i,j))*(zb(i,j)-zb(i,j-1))+(vb(I,j)+vb(i,j+1))*(zb(i,j+1)-zb(i,j)) & +2.0*(zb(i,j)-zo)*(ub(i+1,j)-ub(i-1,j)+vb(i,j+1)-vb(i,j-1))) zc(i,j)=za(i,j)+h*dt enddo enddo return end ! time smoothimg subroutine ts(ua,ub,uc,va,vb,vc,za,zb,zc,s,m,n) dimension ua(m,n),va(m,n),za(m,n),ub(m,n),vb(m,n),zb(m,n),uc(m,n),vc(m,n),zc(m,n) m1=m-1 n1=n-1 do i=2,m1 do j=2,n1 ub(i,j)=ub(i,j)+s*(ua(i,j)+uc(i,j)-2.0*ub(i,j))/2.0 vb(i,j)=vb(i,j)+s*(va(i,j)+vc(i,j)-2.0*vb(i,j))/2.0 zb(i,j)=zb(i,j)+s*(za(i,j)+zc(i,j)-2.0*zb(i,j))/2.0 enddo enddo return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! space smoothing for internal points 区域内5点平滑(正逆平滑) ! 请同学编写区域内5点平滑(正逆平滑)的子程序!!!应用书中(4.126)式 ! l=1为只执行正平滑,l=2为执行正逆平滑 subroutine ssip(a,w,s,m,n,k,l) dimension a(m,n),w(m,n) end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! space smoothing for boundary points 边界九点平滑 subroutine ssbp(a,w,s,m,n) dimension a(m,n),w(m,n) m1=m-1 m3=m-3 n1=n-1 n2=n-2 n3=n-3 do i=2,m1 do j=2,n1,n3 w(i,j)=a(i,j)+0.5*s*(1.0-s)*(a(i-1,j)+a(i+1,j)+a(i,j-1)+a(i,j+1)-4.0*a(i,j)) & +0.25*s*s*(a(i-1,j-1)+a(i-1,j+1)+a(i+1,j-1)+a(i+1,j+1)-4.0*a(i,j)) enddo enddo do i=2,m1,m3 do j=3,n2 w(i,j)=a(i,j)+0.5*s*(1.0-s)*(a(i-1,j)+a(i+1,j)+a(i,j-1)+a(i,j+1)-4.0*a(i,j)) & +0.25*s*s*(a(i-1,j-1)+a(i-1,j+1)+a(i+1,j-1)+a(i+1,j+1)-4.0*a(i,j)) enddo enddo do i=2,m1 do j=2,n1,n3 a(i,j)=w(i,j) enddo enddo do i=2,m1,m3 do j=3,n2 a(i,j)=w(i,j) enddo enddo return end ! transmiting arrays 数组传送 subroutine ta(ua,va,za,ub,vb,zb,m,n) dimension ua(m,n),va(m,n),za(m,n),ub(m,n),vb(m,n),zb(m,n) do i=1,m do j=1,n ua(i,j)=ub(i,j) va(i,j)=vb(i,j) za(i,j)=zb(i,j) enddo enddo return end ! transmiting boundary valaus 赋固定边界值 subroutine tbv(ua,va,za,ub,vb,zb,m,n) dimension ua(m,n),va(m,n),za(m,n),ub(m,n),vb(m,n),zb(m,n) m1=m-1 n1=n-1 do i=1,m do j=1,n,n1 ua(i,j)=ub(i,j) va(i,j)=vb(i,j) za(i,j)=zb(i,j) enddo enddo do i=1,m,m1 do j=1,n ua(i,j)=ub(i,j) va(i,j)=vb(i,j) za(i,j)=zb(i,j) enddo enddo return end ! printing variables 打印积分步数 subroutine pv(na,nb) write(*,'(5x,3hna=,i3,5x,3hnb=,i2/)')na,nb return end 将以上Fortran代码改写为对应python代码,同时编写五点平滑子程序、地转风初值子程序, 其中dat文件存放于D:\HuaweiMoveData\Users\80704\Desktop\实习内容中
### 将Fortran代码转换为Python代码并实现五点平滑子程序和地转风初值子程序
#### 转换策略概述
为了将Fortran代码成功转换为Python代码,可以采用逐步解析的方式。由于Fortran在科学计算领域具有悠久历史,许多经典算法可以直接映射到现代Python库中使用的函数[^1]。特别是涉及矩阵运算的部分,可以通过NumPy和SciPy等高效工具替代原始Fortran逻辑。
---
#### 实现五点平滑子程序 (ssip)
以下是基于Python的五点平滑子程序实现:
```python
import numpy as np
def ssip(data):
"""
五点平滑子程序实现。
参数:
data (numpy.ndarray): 输入的一维数组或二维网格数据
返回:
smoothed_data (numpy.ndarray): 平滑后的数据
"""
rows, cols = data.shape
smoothed_data = np.zeros_like(data)
for i in range(1, rows - 1):
for j in range(1, cols - 1):
# 计算中心点及其周围四个邻居点的平均值
smoothed_value = (
data[i, j] * 4 +
data[i - 1, j] +
data[i + 1, j] +
data[i, j - 1] +
data[i, j + 1]
) / 9.0
smoothed_data[i, j] = smoothed_value
return smoothed_data
```
此实现假设输入是一个二维网格 `data`,并且边界点未被处理。可以根据具体需求调整边界的处理方式[^2]。
---
#### 地转风初值子程序 (cgw)
以下是对地转风初值子程序的Python实现:
```python
def cgw(u, v, f, gph, dx, dy):
"""
地转风初值子程序实现。
参数:
u (numpy.ndarray): 初始东向风场
v (numpy.ndarray): 初始北向风场
f (float or numpy.ndarray): 科氏参数
gph (numpy.ndarray): 地势高度场
dx (float): 经度方向格距
dy (float): 纬度方向格距
返回:
ug (numpy.ndarray): 地转东向风场
vg (numpy.ndarray): 地转北向风场
"""
# 常数定义
gravity = 9.81 # 重力加速度 [m/s^2]
# 初始化地转风场
ug = np.zeros_like(u)
vg = np.zeros_like(v)
# 中心差分法求解地转风
for i in range(1, u.shape[0] - 1):
for j in range(1, u.shape[1] - 1):
dphi_dx = (gph[i, j + 1] - gph[i, j - 1]) / (2 * dx) # 高度梯度沿经向
dphi_dy = (gph[i + 1, j] - gph[i - 1, j]) / (2 * dy) # 高度梯度沿纬向
ug[i, j] = -(gravity / f) * dphi_dy
vg[i, j] = (gravity / f) * dphi_dx
return ug, vg
```
该实现利用了中心差分法来近似偏导数,并结合科氏参数 `f` 和重力常数 `gravity` 来计算地转风场[^3]。
---
#### 数据文件路径设置
指定 `.dat` 文件路径如下所示:
```python
file_path = r"D:\HuaweiMoveData\Users\80704\Desktop\实习内容"
```
注意:需确保路径中的反斜杠 `\` 使用双反斜杠 `\\` 或者前缀 `r` 表示原始字符串,以避免转义字符错误。
---
#### 总结
通过上述方法,可以有效地将Fortran代码的核心功能移植至Python环境。对于复杂的数值计算场景,优先考虑使用 NumPy 和 SciPy 提供的功能模块,从而减少手动编码的工作量并提升性能[^1]。此外,在实际应用过程中可能还需要进一步优化边界条件处理以及验证结果准确性。
---
阅读全文
相关推荐















