利用griddata进行插值

利用griddata进行插值因为最近在做算法优化,所以对数据统一性有一定要求,在最近的研究中主要用一个简单的最近邻插值对数据集进行降尺度处理。主要运用到的函数时scipy里面的griddata第一步:导入相关库importxarrayasxrfromscipy.interpolateimportgriddata#插值函数importnumpyasnp第二步:给出插值到的经纬度信息(目标经纬度)mask_tmp=xr.open_dataset(‘G:/China/temperatu

大家好,又见面了,我是你们的朋友全栈君。

因为最近在做算法优化,所以对数据统一性有一定要求,在最近的研究中主要用一个简单的最近邻插值对数据集进行降尺度处理。

主要运用到的函数时scipy里面的
griddata

griddata函数讲解

`scipy.interpolate.griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)`
  • points:2-D ndarray of floats with shape (n, D), 或 length D tuple of 1-D ndarrays with shape (n,).
  • values:ndarray of float 或 complex, shape (n,)
  • xi:2-D ndarray of floats with shape (m, D), 或 length D tuple of ndarrays broadcastable to the same shape.
  • method:{‘linear’, ‘nearest’, ‘cubic’}, 可选参数

nearest:返回最接近插值点的数据点的值
linear:线性插值
cubic:三次样条插值

第一步:导入相关库

import xarray as xr    
from scipy.interpolate import griddata    # 插值函数
import numpy as np

第二步:给出插值到的经纬度信息(目标经纬度)

mask_tmp = xr.open_dataset('G:/China/temperature_max/nc/2000/tmx_2000_1.nc')
# 待插值的标准经纬度
mask_tmp_lon = mask_tmp.lon.values  # (7680,) 一维
mask_tmp_lat = mask_tmp.lat.values  # (4717,) 一维
mask_LON,mask_LAT = np.meshgrid(mask_tmp_lon,mask_tmp_lat)  # (4717, 7680) 生成经纬度网格

第三步:待插值数据

rad = xr.open_dataset('D:/Users/62692/Desktop/rad.nc')
# 辐射数据经纬度
rad_lon = rad.lon.values    # 辐射数据经度 (641,)
rad_lat = rad.lat.values    # 辐射数据纬度 (394,)
rad_LON, rad_LAT = np.meshgrid(rad_lon,rad_lat)   # (394, 641)
rad_LON = rad_LON.ravel().reshape(-1,1)  # 展平 (252554, 1)
rad_LAT = rad_LAT.ravel().reshape(-1,1)  # 展平 (252554, 1)
rad_values = rad['rad'].values     # 需要插值的辐射数据 (394, 641)
points = np.concatenate([rad_LON,rad_LAT],axis = 1)   # (252554, 2)

第四步:插值

# 插值
data = griddata(points, rad_values.ravel(),(mask_LON,mask_LAT),method='nearest')  # 用最近邻插值即可
# rad_values.ravel() (252554,)
# 这里用最邻近主要考虑到辐射数据的连续性变化,对于线性插值或者三次插值并没有多大影响

汇总成函数

''' Created on 1 23, 2022 @author: GongHaixing 将一个文件夹里面所有的nc文件进行插值 '''
def interp2D(maskpath,mask_lon='lon',mask_lat='lat',inputpath='', outputpath='',data_lon='lon',data_lat='lat',variable='',interp_method='nearest',save=True):
    """输入插值目标的相关信息以及需要插值的数据 :maskpath: 需要插值到对应数据的数据路径 :mask_lon: 标准数据的经度名称,比如:x,lon :mask_lat: 标准数据的纬度名称,比如:y,lat :inputpath: 需要做插值处理的nc文件所在的目录 :outputpath: 插值完nc文件保存的路径,注意要是'/' :data_lon: 需要做插值数据经度名称,比如:'x','lon' :data_lat: 需要做插值数据经度名称,比如:'y','lat' :variable:需要做插值数据变量的名称,比如:'tmp','ndvi' :interp_method: griddata的插值方法,比如:'nearest','linear','cubic' :save:是否对文件进行存储 """
    #导入相关库
    import xarray as xr
    import os
    from scipy.interpolate import griddata    # 插值函数
    import numpy as np

    ### 目标插值
    mask_data = xr.open_dataset(maskpath)
    mask_data_lon = mask_data[mask_lon].values
    mask_data_lat = mask_data[mask_lat].values
    mask_LON,mask_LAT = np.meshgrid(mask_data_lon,mask_data_lat)
    mask_LON1 = np.array(mask_LON)
    mask_LAT1 = np.array(mask_LAT)

    ### 插值对象
    os.chdir(inputpath)                            # 给出nc文件所在的目录(路径)
    files = os.listdir()
    savepath = outputpath
    for file in files:
        inputfile = xr.open_dataset(file)
        inputfile_lon = inputfile[data_lon].values                                 # 数据的经纬度
        inputfile_lat = inputfile[data_lat].values
        inputfile_LON, inputfile_LAT = np.meshgrid(inputfile_lon,inputfile_lat)
        inputfile_LON = inputfile_LON.ravel().reshape(-1,1)
        inputfile_LAT = inputfile_LAT.ravel().reshape(-1,1)
        inputfile_values = np.array(inputfile[variable].values,dtype=np.float32)     # 需要插值的数据
        points = np.concatenate([inputfile_LON,inputfile_LAT],axis = 1)
        # 插值
        print('开始对'+file+'进行插值')
        inputfile_interp = griddata(points, inputfile_values.ravel(),(mask_LON1,mask_LAT1),method=interp_method).astype(np.float32)  # 用最近邻插值即可,不然数据会nan
        outfile= xr.Dataset(
                      data_vars = { 
   variable:(('lat','lon'),inputfile_interp)},
                      coords = { 
   
                                'lon':('lon',np.array(mask_data_lon)),
                                'lat':('lat',np.array(mask_data_lat))}
                                )
        if save == True:
            outfile.to_netcdf(outputpath+'/'+file)
            print(file+'已经插值成功,且已经保存到'+outputpath+'路径下')
        else:
            print(file+'已经插值成功,但是我没有保存文件')
from interp2D import *
import xarray as xr
import os
import pandas as pd
from scipy.interpolate import griddata    # 插值函数


maskpath = 'H:/China/temperature_max/nc/2000/tmx_2000_1.nc'
mask_lon='lon'
mask_lat='lat'
inputpath='H:/China/LAI/nc'
outputpath='H:/China/LAI/nc_1km'
data_lon='lon'
data_lat='lat'
variable='LAI'
interp_method='nearest'
save=True

interp2D(maskpath=maskpath,mask_lon='lon',mask_lat='lat',inputpath=inputpath, outputpath=outputpath,data_lon='lon',data_lat='lat',variable=variable,interp_method='nearest',save=True)

结果对比

插值前(10km)

在这里插入图片描述

插值后(1km)

在这里插入图片描述

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请联系我们举报,一经查实,本站将立刻删除。

发布者:全栈程序员-站长,转载请注明出处:https://javaforall.net/141272.html原文链接:https://javaforall.net

(0)
全栈程序员-站长的头像全栈程序员-站长


相关推荐

  • 计算机科学丛书豆瓣_有关计算机科学与技术专业的书籍

    计算机科学丛书豆瓣_有关计算机科学与技术专业的书籍1.TheArtofComputerProgramming  Author:Donald.E.Knuth  Website:http://www-cs-faculty.stanford.edu/~knuth/taocp.html  BookInfo:这部书被誉为20世纪最重要的20部著作之一,与Einstein的>并列,是计算机科学领域的权威著作.全书共分7卷,目前已经出版了3卷

    2025年6月16日
    0
  • 彻底搞懂0-1背包问题(动态规划)

    彻底搞懂0-1背包问题(动态规划)看了很多网上的博客,发现对于0-1背包问题很多讲的都很专业,初学者学起来还是比较吃力,今天我就用最简单最形象的语言来描述一下0-1背包问题,为什么不能用贪婪算法,而要选择使用动态规划。首先对于0-1背包问题,我们需要知道的是:每一个物品只有1个,要么全拿,要么不拿,最后使得拿到的物品的总价值最大。假如一个小偷有一个可以容纳4千克的背包,但是发现面前只有有3样物品可以偷:台灯(30元,4千克)、音响(20元,3千克)、充电宝(15元,1千克)(价格和重量可能有点奇怪????)。问,小偷能够偷到的物品的

    2022年7月26日
    8
  • javaweb-spring-52

    javaweb-spring-52

    2021年5月17日
    89
  • Pycharm自动导入模块小技巧「建议收藏」

    Pycharm自动导入模块小技巧「建议收藏」周末愉快!不知道大家周末写不写代码,哈哈,反正我已经加完班回来了,今天分享一个能提高编码效率的小技巧,可能你早就在用了,也可能像我一样刚学会,还是趁热跟大家分享一下。如果能把工具熟练运用,往往能达到事半功倍的效果,Pycharm是很多Python开发者的首选IDE,提供各种快捷键、重构功能、调试技巧等,Python是动态语言,对于自动导入模块没有静态语言那么方便,但有了Pycharm,还是…

    2022年8月27日
    0
  • 中英文字典树_字典树详解

    中英文字典树_字典树详解英文字典树英文字典树的结构图是这样的。按照树型结构存储字符串,每个结点存一个字符,自顶向下做标记的就是词的词尾,比如,app,apple,application,abstract,absorb,block,black,blake…等等介绍一下英文字典树的结点数据结构:1.词频int型变量记录词频2.结点型数组,长度26下标对应0-25(也…

    2022年9月6日
    2
  • JS数组的定义

    JS数组的定义

    2021年8月25日
    53

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注全栈程序员社区公众号