利用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)
全栈程序员-站长的头像全栈程序员-站长


相关推荐

  • 谈谈网站防盗链

    谈谈网站防盗链#开启重写引擎RewriteEngineon#图片防盗链RewriteCond%{HTTP_REFERER}!http://rewrite/.*RewriteRule.*\.(jpg|jpeg|png|gif)daotu.jpg引子:明明引用了一个正确的图片地址,但显示出来的却是一个红叉或写有“此图片仅限于***网站用户交流沟

    2022年7月23日
    11
  • phpstome2021.5.1 激活码(最新序列号破解)[通俗易懂]

    phpstome2021.5.1 激活码(最新序列号破解),https://javaforall.net/100143.html。详细ieda激活码不妨到全栈程序员必看教程网一起来了解一下吧!

    2022年3月20日
    45
  • qt服务器主动断开tcp连接_qtcpsocket 多线程

    qt服务器主动断开tcp连接_qtcpsocket 多线程简述对于一个C/S结构的程序,客户端有些时候需要实时得知与服务器的连接状态。而对于客户端与服务器断开连接的因素很多,现在就目前遇到的情况进行一下总结。分为下面六种不同情况客户端网线断开客户端网络断开客户端通过HTTP代理连接服务器,代理机器断开代理客户端通过HTTP代理连接服务器,代理机器的网络断开客户端通过HTTP代理连接服务器,代理机器的网线断开服务器断开同时对于以上六种情况又分为连接服务器之…

    2025年10月12日
    3
  • vue子组件调用父组件父页面的方法「建议收藏」

    vue子组件调用父组件父页面的方法「建议收藏」如图:选择城市后,页面重新请求数据!(城市选择是子组件)选择完后不刷新页面,重新调用接口渲染主页面!第一种方法是直接在子组件中通过this.$parent.event来调用父组件的方法父组件<template><div><child></child></div></template&gt…

    2022年10月2日
    6
  • 初中生学java行吗_初中生学Java开发有前景吗?「建议收藏」

    初中生学java行吗_初中生学Java开发有前景吗?「建议收藏」初中生学Java开发有前景吗?学Java手机开发有什么优势?Java手机开发工程师前景好吗?据北大青鸟华工学校的专业老师分析。作为唯一在互联网上开发的语言,Java平台以其移动性、安全性和开放性受到追捧。据IDC预计,自2001年起的其后5年内,采用Java的IT产品的价值将翻番,在2006年将达到4.53亿美元,年增长率为14.9%。截止到2003年5月,Java注册开发商超过300万人,对JR…

    2022年7月26日
    7
  • OpenResty 最佳实践学习–实战演习笔记(3)

    本篇总结来自 OpenResty(Nginx+Lua)开发入门 ,基本的代码没有改动,主要是自己实际动手操作,测试 Nginx Lua API !我们需要接收请求、处理并输出响应。而对于请求我们需要获取如请求参数、请求头、Body体等信息;而对于处理就是调用相应的Lua代码即可;输出响应需要进行响应状态码、响应头和响应内容体的输出。因此我们从如上几个点出发即可。接收请求我自己进行测试,在我安装的

    2022年2月26日
    51

发表回复

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

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