请教！NCL中用landsea mask去掉海洋格点，遇到未知问题

liao0126

最后在师兄的指点下，用海洋网格缺测的气候态数据除自身作为mask文件，成功把海洋网格mask掉了！
但是原来遇到的问题仍然没有解决，如果有遇到同样问题的朋友，希望能一起探讨！

                               
登录/注册后可看大图

                               
登录/注册后可看大图

各位好，想请教一下关于NCL中用landsea mask去掉海洋格点的问题。

我的资料是1901-2018年中国区域（lat：10.25N-59.75N | lon: 60.25E-149.25E ）逐月的温度数据，分辨率为0.5°*0.5°。
landsea mask文件用的是老师给的0.25°*0.25°的文件。

我选用了温度资料中1961-2018年的部分，用landsea mask去掉海洋部分格点，再对每一年求平均得到逐年数据，最后对每个格点进行一元线性回归得到斜率即温度变化率（*10得十年变化率）。

我一开始用的landsea mask是ncl自带的mask文件，分辨率是1*1，画出来的图沿岸边缘呈锯齿状，影响美观和分析。

后来改用老师给的0.25°*0.25°的文件，mask之后发现画出来的温度变化率图基本是全部空白，只有一两个格点有数值和颜色。计算出来的平均温度变化率也不符合认知，但是全程没有报错。找了很久没有找到问题在哪，希望各位能帮我看看！

我还尝试把mask完还没处理的温度数据ano(year,month,lat,lon)随意取了一个月份，画出来的结果也有很多本来应该也有数据的网格点没有数据。



两个资料ncl_filedump之后的信息附在这里：

①温度
Variable: f
Type: file
filename:        idw_tavg_ano_m4n10
path:        /cygdrive/d/lls/NCL/yanxi/nc/idw_tavg_ano_m4n10.nc
   file global attributes:
   dimensions:
      lat = 100
      lon = 180
      month = 12
      year = 118
   variables:
      double ano ( year, month, lon, lat )
         units :        degC
         _FillValue :        -999
         long_name :        ano

      double lat ( lat )
         units :        degress_north
         long_name :        latitude

      double lon ( lon )
         units :        degress_east
         long_name :        longitude

      double month ( month )
         units :        month
         long_name :        month

      double year ( year )
         units :        year
         long_name :        year


②landsea mask
Variable: f
Type: file
filename:        land-sea_mask
path:        /cygdrive/d/lls/NCL/yanxi/nc/land-sea_mask.nc
   file global attributes:
      Conventions : CF-1.6
      history : 2020-06-05 05:26:21 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data2/adaptor.mars.internal-1591334780.811916-17544-21-ffa08836-6e9f-4d3e-89f4-5e71f1a263ab.nc /cache/tmp/ffa08836-6e9f-4d3e-89f4-5e71f1a263ab-adaptor.mars.internal-1591334780.812469-17544-5-tmp.grib
   dimensions:
      longitude = 1440
      latitude = 721
      time = 1
   variables:
      float longitude ( longitude )
         units :        degrees_east
         long_name :        longitude

      float latitude ( latitude )
         units :        degrees_north
         long_name :        latitude

      integer time ( time )
         units :        hours since 1900-01-01 00:00:00.0
         long_name :        time
         calendar :        gregorian

      short lsm ( time, latitude, longitude )
         scale_factor :        1.525948758640685e-05
         add_offset :        0.4999923702562068
         _FillValue :        -32767
         missing_value :        -32767
         units :        (0 - 1)
         long_name :        Land-sea mask
         standard_name :        land_binary_mask



③下附计算和画图过程：

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" ; 如需将台站资料插值至格点资料，则必须加载该ncl文件
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


begin

  latS     =  10.25               
  latN     =  59.75
  lonL     =  60.25
  lonR     =  149.75
  ntStrt   =  1961
  ntLast   =  2018

;--------------读取数据 （数据文件1961.1-2018.12，逐月数据）

  f = addfile("idw_tavg_ano_m4n10.nc","r")
  ano   = f->ano({year|ntStrt:ntLast},month|:,{lon|lonL:lonR},{lat|latS:latN})   ;ano(year,month,lon,lat)
  year  = f->year({year|ntStrt:ntLast})
  month = f->month(month|:)
  lon   = f->lon({lon|lonL:lonR})
  lat   = f->lat({lat|latS:latN})

  printVarSummary(ano)                      ;ano(year,month,lon,lat)
  printVarSummary(year)


;--------------mask掉海洋格点的数值

;------读取landsea mask文件，short
  a      = addfile("land-sea_mask.nc","r")    ;landsea mask   0.25*0.25
  lsdata = shorttobyte (a->lsm(0,:,:))        ;lsm(time,latitude{90:-90},longitude{0:359.75})    lsdata(latitude,longitude)
  printVarSummary(lsdata)
        
  longitude  = a->longitude
  latitude   = a->latitude
  printVarSummary(latitude)
  printVarSummary(longitude)

;------shorttobyte时不提取metadata，需加上
  lsdata!0   = "lat"
  lsdata&lat = latitude

  lsdata!1   = "lon"
  lsdata&lon = longitude

  lsdata&lat@units = "degrees_north"
  lsdata&lon@units = "degrees_east"

  printVarSummary(lsdata)                  

;------mask掉海洋格点后得到ano1(year,month,lat,lon)
  lsmask = landsea_mask(lsdata,ano&lat,ano&lon)  
  printVarSummary(lsmask)
  ano1 = mask(ano(year|:,month|:,lat|:,lon|:),conform(ano(year|:,month|:,lat|:,lon|:),lsmask,(/2,3/)).ge.1,False)   ;ano1(year,month,lat,lon)

  printVarSummary(ano1)


;--------------对月求平均，将逐月数据转为逐年数据

  anoY = dim_avg_n_Wrap(ano1,1)                             ;anoY(year,lat,lon)
  anoY@_FillValue=-999                              
  printVarSummary(anoY)


;--------------计算每一个网格点的一元线性拟合后的斜率，即十年增温率

  trends = new( (/dimsizes(lat),dimsizes(lon)/) , double)   ;trends(lat,lon)
  trends@_FillValue=-999

  trends  = regCoef_n(year,anoY,0,0)                        ;trends(lat,lon)
  trends  = trends*10                                       ;得十年增温率
  printVarSummary(trends)


;-------------设置trends属性

  trends!0          = "lat"
  trends&lat        = lat
  trends!1          = "lon"
  trends&lon        = lon

  trends&lat@units  = "degrees_north"
  trends&lon@units  = "degrees_east"
  printVarSummary(trends)


;--------------计算最大、最小增温率

  print("max trends = " + max(trends))                     
  print("min trends = " + min(trends))

;--------------计算平均增温率

avg_trends=dim_avg_n_Wrap(dim_avg_n_Wrap(trends,1),0)
print("avg trends = " + avg_trends)


;--------------画图

  wks  = gsn_open_wks("png","D:/lls/NCL/yanxi/avg/trends1961to2018" )

  res=True
  res@gsnDraw      = False
  res@gsnFrame     = False
  res@gsnAddCyclic = False         
  res@gsnMaximize  = True
  res@gsnCenterString = "average trends: "+ tofloat(avg_trends)        ; add title


  res@tiMainString =  "China land surface air temperature trends, " + (ntStrt) +"-"+ (ntLast)


  res@cnFillOn             =  True
  res@cnFillDrawOrder      = "PreDraw"
  res@cnFillPalette        = "WhiteYellowOrangeRed"

  res@cnLevelSelectionMode = "ManualLevels"

  res@cnMinLevelValF       =  0.05
  res@cnMaxLevelValF       =  0.45

  res@cnLevelSpacingF      =  0.05

  res@cnLinesOn             = False


  res@mpMaxLatF = latN  
  res@mpMinLatF = latS
  res@mpMaxLonF = lonR
  res@mpMinLonF = lonL

  res@mpLandFillColor = "Transparent"

  res@mpOceanFillColor = "Gray40"

  res@mpOutlineBoundarySets = "National"


  plot=gsn_csm_contour_map(wks,trends,res)    ; create the plot


  draw(plot)
  frame(wks)
  delete(plot)

end
下面第一个图是一开始用ncl自带的1*1的mask画的，第二个图是用现在0.25*0.25的mask画的。

1518982

你好，请问可以分享一下0.25*0.25的land-sea_mask.nc文件吗，我最近在处理海洋数据，需要去掉陆地格点，但是官网上landsea_mask 分辨率太低，万分感谢

wangzhen

您好，请问NCL自带的mask文件在哪里可以下载呢？

点hewen

1518982 发表于 2020-11-24 21:55
你好，请问可以分享一下0.25*0.25的land-sea_mask.nc文件吗，我最近在处理海洋数据，需要去掉陆地格点，但 ...

你好，我最近也在处理海洋数据，需要把陆地格点去掉，想请问一下你用了什么方法，谢谢！

liao0126

wangzhen 发表于 2021-6-30 21:29
您好，请问NCL自带的mask文件在哪里可以下载呢？

这个应该安装NCL的时候就有了？系统自带的 没有另外下载