第一步:找到站点数据对应的格点数据,其本质是找到离站点经纬度最近的格点位置。
clc;
clear;
filename1 = 'D:\Stalocation_GC.xlsx';
filename2 ='D:\datagrid\wrfd03.dat';
sta = xlsread(filename1);
stnum = size(sta,1);
id = zeros(stnum,2);
fid = fopen(filename2,'r','b');
data2 = fread(fid,'real*4');
fclose(fid);
data2 = reshape(data2,228,150,84);
lon = squeeze(data2(:,:,83));
lat = squeeze(data2(:,:,82));
dd = zeros(228,150);
nx = 228;
ny = 150;
for ist=1:stnum
for i=1:nx
for j=1:ny
dd(i,j) = sqrt((sta(ist,1)-lon(i,j))^2+(sta(ist,2)-lat(i,j))^2);
end
end
mm = min(min(dd));
[row,col] = find(dd==mm);
id(ist,1) = row;
id(ist,2) = col;
disp([lon(row,col),lat(row,col)]);
end
第二步:从第一步知道idx = 146、idy=110,读取二进制数据并提取。不得不说matlab和python一样慢,有大佬知道提高读取速度的方法吗?求分享
matlab版
clc;
clear;
tic;
filename2='D:\datagrid\wrfd03.dat';
file = dir("E:\dry_data\*.grd"); ## 从2019年12月7日-12月31日
N = length(file);
bc = zeros(600,20);
rh = zeros(600,20);
t = zeros(600,20);
nx = 228;
ny = 150;
vars_vertical = 2905;
idx = 146;
idy = 110;
for i =1:N
temp = strcat(file(1).folder,"\",file(i).name);
disp(temp);
fid = fopen(temp,'r','b');
data1 = fread(fid,'float');
fclose(fid);
data1=reshape(data1,nx,ny,vars_vertical);
%bc = squeeze(sum(data1(:,:,:,1),3));
bc(i,:) = squeeze(data1(idx,idy,2286:2305));
rh(i,:) = squeeze(data1(idx,idy,101:120));
t(i,:) = squeeze(data1(idx,idy,81:100));
%clear data1;
end
toc;
%%% matlab读这么多文件需要58分钟,真慢。
python版
## 1.导入库
import numpy as np
import pandas as pd
import matplotlib as mpl
import gc
import os
import math
import glob
import datetime
import time
from matplotlib import pyplot as plt
from datetime import datetime,timedelta
from pathlib import Path
## 2.读文件
p = Path(r"E:\dry_data")
FileList = list(p.glob("*.grd"))
bc = np.zeros((600,20));
rh = np.zeros((600,20));
t = np.zeros((600,20));
i=-1
idx = 145 ## 一定要注意python是从0开始计数啊,朋友们。
idy = 109
nx = 228;
ny = 150;
vars_vertical = 2905
T1=time.time()
for file in FileList[:]:
print(i+1,file)
data_raw = np.fromfile(file,dtype=">f")
data = data_raw.reshape((vars_vertical,ny,nx),order="C")
i=i+1
rh[i,:]= data[100:120,idy,idx]
t[i,:] = data[80:100,idy,idx]
bc[i,:]= data[2285:2305,idy,idx]
T2=time.time()
print('程序运行时间:%s秒' % ((T2 - T1)))
## 运行时间大约也要一个小时~,气。
网友评论