美文网首页
如何根据站点数据提取模式数据

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

作者: Aerosols | 来源:发表于2022-04-17 21:44 被阅读0次

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

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

    相关文章

      网友评论

          本文标题:如何根据站点数据提取模式数据

          本文链接:https://www.haomeiwen.com/subject/qbkisrtx.html