- 积分
- 541
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-10-19
- 最后登录
- 1970-1-1
![[小馒头面包_] 粉丝数:1 微博数:4 新浪微博达人勋](source/plugin/sina_login/img/light.png)
|

楼主 |
发表于 2017-3-29 18:22:08
|
显示全部楼层
u850_JJA@_FillValue = -9.96921e+36
v850_JJA@_FillValue = -9.96921e+36
MASK_INSIDE = False ; Whether to mask data inside or outside the given geographical area.
minlat = -25
maxlat = 45
minlon = 60
maxlon = 160
nlat = 71
nlon = 101
;---Create dummy 1D lat/lon arrays
lat1d = fspan(minlat,maxlat,nlat)
lon1d = fspan(minlon,maxlon,nlon)
lat1d@units = "degrees_north"
lon1d@units = "degrees_east"
printVarSummary(lat1d)
printVarSummary(lon1d)
;Open shapefile and read Qinghai Tibet Plateau lat/lon values.
f = addfile("/lustre/home/niehw/moshi/wind1/Qinghai-Tibet_Plateau/Qinghai-Tibet_Plateau.shp","r")
;dir = "/lustre/home/niehw/moshi/wind1/DBATP/"
;filename = dir + "DBATP_Polygon.shp"
mrb_lon = f->x
mrb_lat = f->y
printVarSummary(mrb_lon)
printVarSummary(mrb_lat)
nmrb = dimsizes(mrb_lon)
min_mrb_lat = min(mrb_lat)
max_mrb_lat = max(mrb_lat)
min_mrb_lon = min(mrb_lon)
max_mrb_lon = max(mrb_lon)
if(MASK_INSIDE) then
;---Start with data filled in.
data_mask1 = u850_JJA
data_mask2 = v850_JJA
else
;---Start with data all missing
data_mask1 = new(dimsizes(u850_JJA),typeof(u850_JJA),u850_JJA@_FillValue)
copy_VarCoords(u850_JJA,data_mask1)
data_mask2 = new(dimsizes(v850_JJA),typeof(v850_JJA),v850_JJA@_FillValue)
copy_VarCoords(v850_JJA,data_mask2)
end if
printVarSummary(data_mask1)
printVarSummary(data_mask2)
ilt_mn = ind(min_mrb_lat.gt.lat1d)
ilt_mx = ind(max_mrb_lat.lt.lat1d)
iln_mn = ind(min_mrb_lon.gt.lon1d)
iln_mx = ind(max_mrb_lon.lt.lon1d)
ilt1 = ilt_mn(dimsizes(ilt_mn)-1) ; Start of lat box
iln1 = iln_mn(dimsizes(iln_mn)-1) ; Start of lon box
ilt2 = ilt_mx(0) ; End of lat box
iln2 = iln_mx(0) ; End of lon box
if(MASK_INSIDE) then
;---Put missing values in the areas that we want masked.
do ilt=ilt1,ilt2
do iln=iln1,iln2
if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
data_mask1(ilt,iln) = u850_JJA@_FillValue
data_mask2(ilt,iln) = v850_JJA@_FillValue
end if
end do
end do
else
;---Put data back in the areas that we don't want masked.
do ilt=ilt1,ilt2
do iln=iln1,iln2
if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
data_mask1(ilt,iln) = u850_JJA(ilt,iln)
data_mask2(ilt,iln) = v850_JJA(ilt,iln)
end if
end do
end do
end if
|
|