对EASE-Grid Snow Water Equivalent的nc数据读取、网格转换和作图

逝风丝痕

对EASE-Grid Snow Water Equivalent的nc数据

----读取、网格转换和作图

一．数据

NSIDC的swe数据网址Index of /swe(globsnow.info)

该数据产品格式为Equal-AreaScalable Earth Grid (EASE-Grid)，其中archive_v3.0的数据缺少lat,lon属性，且每个格点的经度和纬度均有一个二维数组指定，这种类型数据属于结构网格中的曲线网格，所以在使用前需要进行预处理，将其转换成直线网格，如常见的等经纬度网格，以方便后续计算。

二．步骤

1 批量读取源数据

    path                 = "./data/"

    p1                    = systemfunc("ls " + path + "*_northern_hemisphere_monthly_swe_0.25grid.nc")

    print(p1)

    f1                    = addfiles(p1, "r")

    ListSetType(f1, "join")

    swe                  = f1[:]->swe(0,:,:)

    printVarSummary(swe)

2 添加样本数据的lat，lon

初步判断archive_v2.0的对应属性可以简单套用，通过NCL相关函数可以解决；

    difi                       = ("./sample/data/")  

    sample                 = addfile(difi+"GlobSnow_SWE_L3B_monthly_201002_v2.0.nc", "r")

    swe2                    = sample->SWE_avg

    easelat                 = sample->lat

    easelat@units       = "degrees_north"

    easelon                = sample->lon

    easelon@units      = "degrees_east"

    printVarSummary(easelat)

    printVarSummary(easelon)

    swe_new              = tofloat(swe)

    copy_VarCoords(swe, swe_new)

    copy_VarMeta(swe2, swe_new)

    swe_new@lat2d    = easelat

    swe_new@lon2d   = easelon

3 借助Earth SystemModeling Framework (ESMF)进行网格转换

需要load"$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

网格转换所需属性和参数：

InterpMethod：插值方法。

SrcFileName：该属性为源网格的描述文件，该文件会自动创建，默认值为“source_grid_file.nc”。

DstFileName：该属性为目标网格的描述文件，该文件会自动创建，默认值为“destination_grid_file.nc”。

WgtFileName：插值权重文件，该文件会自动创建，默认值为“weights_file.nc”

ForceOverwrite：强制覆盖重写。

SrcMask2D：将需要遮盖的区域设为0，否则为1。用于对缺省值的处理。

SrcGridLat：源网格的纬度数组。

SrcGridLon：源网格的经度数组。

DstGridType：目标网格的类型。

DstLLCorner：目标网格的左下角点经纬度。

DstURCorner：目标网格的右上角点经纬度。

本次处理主要用到以上属性，其他属性详见官网http://www.ncl.ucar.edu/Document/Functions/ESMF/ESMF_regrid.shtml或施宁等《NCL数据处理与绘图基础教程》P78-P81相关讲解。

    Opt                   = True             ; Options for regridding

    Opt@InterpMethod      = "patch"          ; regrid method

    Opt@SrcFileName       = "EASE_SCRIP0.nc"  ; Output files

    Opt@DstFileName       = "NH_SCRIP0.nc"

    Opt@WgtFileName       = "EASE_0.25_SCRIP_patch0.nc"

    Opt@ForceOverwrite    = True

    Opt@SrcMask2D         = where(.not.ismissing(easelat),1,0)

    Opt@SrcGridLat        = easelat

    Opt@SrcGridLon        = easelon

    Opt@DstGridType       = "0.25deg"        ; Destination grid

    Opt@DstTitle          = "Northern Hemisphere 0.25 resolution"

    Opt@DstLLCorner       = (/ 0.0d,   0.0d/)

    Opt@DstURCorner       = (/90.0d, 359.75d/)

    Opt@Debug             = True

    swe_regrid            = ESMF_regrid(swe_new, Opt)    ; Regrid swe

    printVarSummary(swe_regrid)

4 借助插值权重文件

在第3步的基础上，可以简化后续批量处理的工作

      swe_regrid = ESMF_regrid_with_weights(swe_new,"EASE_0.25_SCRIP_patch0.nc",False)

      printVarSummary(swe_regrid)

5 作图

可参考官网例子www.ncl.ucar.edu/Applications/Scripts/ESMF_regrid_13.ncl

6 批量处理

见附件

Ebony

想问问您V3的坐标变成了一维的X、Y，要如何处理呢？

逝风丝痕

Ebony 发表于 2021-9-8 22:14
想问问您V3的坐标变成了一维的X、Y，要如何处理呢？

如果比较一下V2，您会发现V3少了lat，lon。而V2里的lat，lon是由x，y决定的。V3也是如此，即x，y并不是格点的坐标而是lat，lon的坐标，所以处理V3前需要先把lat，lon赋上，再进行后续的regridding。脚本介绍里已经给了详细的过程。

Ebony

逝风丝痕 发表于 2021-9-9 12:57
如果比较一下V2，您会发现V3少了lat，lon。而V2里的lat，lon是由x，y决定的。V3也是如此，即x，y并不是格 ...

非常感谢！

stepdance

赞一个！！！

meng0109

请问有matlab代码吗？

myj0219

感谢分享！！！