使用ncl在nc文件中提取不规则区域的数据

151170080

在处理nc数据时，有时候分析区域是规则的时候是比较容易处理的。但是有些时候我们分析的区域是不规则的，比如拿全球再分析资料分析中国区域（或者中国某个区域）时。这个时候我们就需要从一个nc数据中扣出我们分析区域的的数据。接下来，我总结了提取步骤供大家讨论使用。

我的上一篇帖子讲述了如何自己制作中国区域（或者某个省份）的shape文件，接下来的提出过程需要我们制作的shape文件，此外，为了能更好地理解脚本思路，我以提取中国区域的内蒙古地区数据为例，进行讲解。
   ——————>   
一、制作待提取区域的数据文件（区域内值为1，区域外值为0）
1. 我自己写的脚本
begin

; ... read shapefile data
      shp_dir = "***/SData/China_kang/"
      shp_file = systemfunc("ls "+shp_dir+"*.shp")
      ;print(""+shp_file)
      f = addfile(shp_file,"r")


; ... read 01D netCDF files 读取的是中国区域的向下短波辐射
      nc_dir = "***/FData/Y_clim/"
      nc_files = systemfunc("ls "+nc_dir+"*.nc")

      a = addfile(nc_files(0),"r")  ; read dlwrf file
      data = a->dlwrf(0,0,:,:)
      LAT = a->lat
      LON = a->lon


; ... Deternube the grid type. rectilinear is default. 将多维数组的经纬度转化为一维数组
      dims = dimsizes(data)
      lat1d = ndtooned(conform_dims(dims,data&$data!0$,0))
      lon1d = ndtooned(conform_dims(dims,data&$data!1$,1))


; ... Read data off the shapefile
      segments = f->segments  
      geometry = f->geometry  
      segsDims = dimsizes(segments)
      geomDims = dimsizes(geometry)

      geom_segIndex = f@geom_segIndex
      geom_numSegs = f@geom_numSegs
      segs_xyzIndex = f@segs_xyzIndex
      segs_numPnts = f@segs_numPnts
      numFeatures = geomDims(0)

      lat = f->y
      lon = f->x

; ... Use the whole shapefile
      min_lat_shp = min(lat)
      max_lat_shp = max(lat)
      min_lon_shp = min(lon)
      max_lon_shp = max(lon)

; ... lat/lon coordinates of data array
      min_lat_data = min(lat1d)
      max_lat_data = max(lat1d)
      min_lon_data = min(lon1d)
      max_lon_data = max(lon1d)

; ... Create array to hold new data mask
      data_mask_1d = new(dimsizes(lat1d),integer)
      data_mask_1d = 0 ; 0 is out, and 1 is in。

; ... Check in each feature
      do nf=0,numFeatures-1
            startSegment = geometry(nf,geom_segIndex)
            numSegments = geometry(nf,geom_numSegs)
            do seg=startSegment, startSegment+numSegments-1
                  startPT = segments(seg,segs_xyzIndex)
                  endPT   = startPT + segments(seg,segs_numPnts)-1
                  lat_sub := lat(startPT:endPT)
                  lon_sub := lon(startPT:endPT)


                  ;print(""+lat_sub(0)+" "+lat_sub(dimsizes(lat_sub)-1))
                  ;print(""+seg+" "+dimsizes(lat_sub))


                  ii_latlon := ind(min(lat_sub) .le. lat1d .and. lat1d .le. max(lat_sub) .and.\
                              min(lon_sub) .le. lon1d .and. lon1d .le. max(lon_sub)）

                  if(any(ismissing(ii_latlon))) then
                        continue
                  end if
                  nii = dimsizes(ii_latlon)


                  do n=0,nii-1
                        iltln = ii_latlon(n)
                        if (gc_inout(lat1d(iltln),lon1d(iltln),lat_sub,lon_sub)) then
                              print("lat "+lat1d(iltln)+"lon "+lon1d(iltln))
                        end if
                  end do
            end do
      end do




; ... Convert a multidimensional array to a one-dimensional array
      data_mask = onedtond(data_mask_1d,dims)
      data_mask!0 = "lat"
      data_mask!1 = "lon"
      data_mask&lat = LAT
      data_mask&lon = LON

; ... Write to netCDF
      system("rm CN_mask.nc")


      fout = addfile("CN_mask.nc","c")


      dim_names = getvardims(data_mask)
      dims_1 = dimsizes(data_mask)
      dim_unlimited = (/False,False/)
      filedimdef(fout,dim_names,dims_1,dim_unlimited)


      filevardef(fout,"data_mask",typeof(data_mask),dim_names)
      filevardef(fout,"lat",typeof(LAT),"lat")
      filevardef(fout,"lon",typeof(LON),"lon")


      filevarattdef(fout,"data_mask",data_mask)
      filevarattdef(fout,"lat",LAT)
      filevarattdef(fout,"lon",LON)


      fout->data_mask = (/data_mask/)
      fout->lat = (/LAT/)
      fout->lon = (/LON/)


      delete(fout)
end


几点解释：
（1）脚本里重要部分的就是根据shape文件中的segments进行循环判断，为什么要这样做呢？
主要是考虑到像内蒙古这么大的polygon里面会有很多个polygons拼接而成的，如果这么大的一个区域仅仅作为一个polygon来判断的话会出现区域提取不正确的问题。
（2）为什么脚本创建了一个nc文件？
这个nc文件中的data_mask变量的分辨率和我要提取的数据的分辨率是一致的。输出文件的原因是为了后面的提取数据节省时间。上面的脚本在在segments循环判断区域内外的时候需要的时间比较长，而我们需要提取的变量通常还会有时间维度，有的甚至有高度维度。如果每一个时间都要提取一边的话，这个时间成本显然是我们不能接受的，所以为了后面提取方便，需要单独输出出来。

2. 利用ncl官网上提供的函数
ncl官方提供了一个补充函数shapefile_mask_data(data,shapefile_name,res),用于提取不规则区域的数据。这个函数是由脚本shapefile_utils.ncl提供的，需要单独下载。详情请参考：http://www.ncl.ucar.edu/Applications/Scripts/shapefiles_20.ncl。我同样建议先提取某一时间的数据，输出为nc文件。以方面后面使用。
函数中使用的res设置推荐如下：
      shres                                 = True
      shres@debug                     = True
      shres@return_mask            = True  ; 0 is out, 1 is in.



二、使用提取的nc文件扣取长时间的数据


begin

; ... read 01D netCDF files 读取的是中国区域的向下短波辐射
      nc_dir = "***/FData/Y_clim/"
      nc_files = systemfunc("ls "+nc_dir+"*.nc")

      a = addfile(nc_files(0),"r") ; read dlwrf file dlwrf(time,alt,lat,lon)
      data = a->dlwrf

; ... read mask data
      mask_f = addfile("***/CN_mask.nc","r")
      mask_data = mask_f->data_mask


      dims = dimsizes(data)
      do nt=0,dims(0)-1
           do na=0,dims(1)-1
                data(nt,na,:,:) = where(mask_data .eq. 1,data(nt,na,:,:),data@_FillValue)
           end do
       end do

; ... 使用mask函数可以替代上面的两层循环
      ;data_nmg = mask(data,mask_data,1)

; ... your analysis
      ...
      ...
      ...

end

由于之前输出了mask data的nc文件，上面这个脚本会很快的执行，大大节省了时间成本。另外，mask函数可以替代上面的两层循环，脚本为：data_nmg = mask(data,mask_data,1)


以上就是我自己扣取不规则区域数据文件的经验，写的过程中难免在表述中会存在问题，如果你使用了上述方法，但是没有准确地扣取到数据，可以在下方留言交流。





来自群组: 南京大学风云英华

大到暴雨

你好，程序运行过程中出现了这个错误，请问你遇到过吗，如何解决呢，gc_inout: the rightmost dimension of lat/lon must be at least three.

TaoYao

你好，程序运行过程中遇到了gc_inout:the rightmost dimension of lat/lon must be at least three  该怎么解决呢

TaoYao

TaoYao 发表于 2021-11-27 10:55
你好，程序运行过程中遇到了gc_inout:the rightmost dimension of lat/lon must be at least three  该怎么 ...

您好，我也遇到了相同的问题，知道怎么解决了吗

dei阿珺

TaoYao 发表于 2021-11-27 10:55
你好，程序运行过程中遇到了gc_inout:the rightmost dimension of lat/lon must be at least three  该怎么 ...

我也是 一脸懵

blue.

我也遇到了这个问题。gc_inout:the rightmost dimension of lat/lon must be at least three。换个shp文件就可以了。

大气环流

blue. 发表于 2023-5-24 11:08
我也遇到了这个问题。gc_inout:the rightmost dimension of lat/lon must be at least three。换个shp文件 ...

你好，遇到了这个问题说明你截取的范围内不闭合的区域，似乎是因为数据里面包含了九段线

novelbean

貌似你从例子中截取了内蒙古区域，然后写入到nc文件，它经纬度不应该是中国区域的LAT和LON，而应该是内蒙古区域的lat和lon吧

; ... Convert a multidimensional array to a one-dimensional array
      data_mask = onedtond(data_mask_1d,dims)
      data_mask!0 = "lat"
      data_mask!1 = "lon"
      data_mask&lat = LAT
      data_mask&lon = LON
你mask后区域应该比原始的区域小吧，它的经纬度值应该如下
data_mask&lat =lat
data_mask&lon = lon
猜对吧？