如何根据站点数据提取模式数据

第一步：找到站点数据对应的格点数据，其本质是找到离站点经纬度最近的格点位置。

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)))
## 运行时间大约也要一个小时~，气。