- 积分
- 3544
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-9-29
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
我想把除了shp文件的区域都填成白色，但为什么效果是反过来的？有高手帮忙看看吗？
代码：
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 "/home/zx/test/sansha/shapefile_utils.ncl"
;==========================================================================
;坐标属性
undef ("setProperty")
procedure setProperty(olon, olat)
begin
olon!0 = "lon"
olon@long_name = "lon"
olon@units = "degrees_east"
olon&lon = olon
olat!0 = "lat"
olat@long_name = "lat"
olat@units = "degrees_north"
olat&lat = olat
end
;数据赋坐标
undef ("setCoordinate")
procedure setCoordinate(data,lon,lat)
begin
data!0 = "lat"
data&lat = lat
data!1 = "lon"
data&lon = lon
end
function get_areas_of_interest(shp_file_name,shp_var_name,opt[1]:logical)
begin
;---Open the shapefile
f = addfile(shp_file_name,"r")
features = f->$shp_var_name$
if(opt.and.isatt(opt,"areas_to_exclude")) then
features@_FillValue = "missing"
do na=0,dimsizes(opt@areas_to_exclude)-1
ii := ind(features.eq.opt@areas_to_exclude(na))
if(.not.any(ismissing(ii))) then
features(ii) = features@_FillValue
end if
end do
return(features(ind(.not.ismissing(features))))
else
return(features)
end if
end
procedure add_shapefile_primitives_by_name(wks,plot,shp_file_name, \
shp_var_name,requested_features,\
opt[1]:logical)
local poly_type, ptres, f, geomDims, numFeatures, features, segments, \
geometry, segsDims, geom_segIndex, geom_numSegs, segs_xyzIndex,\
segs_numPnts, lat, lon, startSegment, numSegments, startPT, endPT
begin
polytype = get_res_value(opt,"polytype","polyline") ; "marker", "polygon"
valid_prim_types = (/"polymarker","polyline","polygon"/)
if(.not.any(polytype.eq.valid_prim_types)) then
print("add_shapefile_primitives_by_name: invalid polytype.")
print(" Must be "+str_join(valid_prim_types,","))
return
end if
;---Read data off the shapefile
f = addfile(shp_file_name,"r")
geomDims = getfilevardimsizes(f,"geometry")
numFeatures = geomDims(0)
features = f->$shp_var_name$
segments = f->segments
geometry = f->geometry
segsDims = dimsizes(segments)
;---Read global attributes
geom_segIndex = f@geom_segIndex
geom_numSegs = f@geom_numSegs
segs_xyzIndex = f@segs_xyzIndex
segs_numPnts = f@segs_numPnts
;---Section to attach polygons to plot.
lon = f->x
lat = f->y
;
; Set custom primitive resources. It doesn't hurt to set, say line
; color, even if you are just drawing markers. They will be ignored.
ptres = True
ptres@gsLineColor = get_res_value(opt,"gsLineColor","darkorchid4")
ptres@gsLineThicknessF = get_res_value(opt,"gsLineThicknessF",10.0)
ptres@gsMarkerIndex = get_res_value(opt,"gsMarkerIndex",16)
ptres@gsFillColor = get_res_value(opt,"gsFillColor","white")
do i=0,numFeatures-1
if(.not.any(features(i).eq.requested_features)) then
continue
end if
startSegment = geometry(i, geom_segIndex)
numSegments = geometry(i, geom_numSegs)
do seg=startSegment, startSegment+numSegments-1
startPT = segments(seg, segs_xyzIndex)
endPT = startPT + segments(seg, segs_numPnts) - 1
plot@$unique_string("primitive")$ = gsn_add_primitive(wks,plot,lon(startPT:endPT),lat(startPT:endPT),False,polytype,ptres)
end do
end do
end
begin
filepath = "/home/zx/test/sansha/config.txt"
outputfile="TP"
argu = asciiread(filepath,-1,"string")
lines = asciiread(argu(2)+argu(3),-1,"float")
delim=" "
data0 = stringtofloat(str_get_field(lines(22::),1,delim))
data=onedtond(data0,(/241,321/))
longrid = 0.125 ;经度格距
latgrid = -0.125 ;纬度格距
lon0 = 100 ;起始经度
lon1 = 140 ;结束经度
lat0 = 40 ; 起始纬度
lat1 = 10 ;结束纬度
nx=floattointeger(abs((lon1-lon0)/longrid)+1)
ny=floattointeger(abs((lat1-lat0)/latgrid)+1)
lon=fspan(lon0,lon1,nx)
lat=fspan(lat0,lat1,ny)
setProperty(lon,lat)
setCoordinate(data,lon,lat)
x0=lon0
x1=lon1
y0=lat1
y1=lat0
imgWidth=nx
imgHeight=ny
;产品生成目录
outputPath="/home/zx/test/sansha"
; checkOrCreateDir(outputPath)
colors = (/(/255,255,255/),(/ 0, 0, 0/),(/166,243,141/),(/58,186,65/),(/97,184,253/),(/1,2,253/),(/255,0,254/)/)*1.0
cmap = colors/255.
;工作台大小与图片类型设置
;rscan = (/1,.9,.8,.7,.6,.4,.1/)
;contour = gsn_csm_contour_map(wks,data({y0:y1},{x0:x1}),res)
; contour=gsn_csm_contour(wks, data({y0:y1},{x0:x1}),mpres)
shapefile = "/home/zx/test/sansha/zhangzhou/zhangzhou.shp"
areas_of_interest = (/"ZZ"/)
opt = True
opt@minlat = y0 ; Setting these four resources makes
opt@maxlat = y1 ; the function run slightly faster,
opt@minlon = x0 ; b/c it won't check the whole
opt@maxlon = x1 ; shapefile.
opt@debug = True
opt@shape_var = "NAME" ; var name that contains region names
opt@keep = False ; throw away points inside desert regions
opt@shape_names = areas_of_interest
data_mask = shapefile_mask_data(data,shapefile,False)
;data_mask = shapefile_mask_data(data,shapefile,True)
wks_type = "png"
wks_type@wkWidth = 1000 ;工作台宽度
wks_type@wkHeight = 1000 ;工作台高度
wks = gsn_open_wks(wks_type,outputPath+"/"+outputfile)
gsn_define_colormap(wks,cmap)
;色标设置
res = True
res@gsnFrame = False
res@gsnDraw = False
; 左上角坐标
res@vpXF = 0.1
res@vpYF = 0.9
; 调节矩形的长宽
res@vpWidthF = 0.8
res@vpHeightF = 0.8
;色斑图
res@cnFillOn = True
res@cnLinesOn = False
res@cnLineLabelsOn = False
res@cnLevelSpacingF = 0.2
res@cnFillMode ="RasterFill"
;去掉等值线信息标识
res@cnInfoLabelOn = False
res@cnConstFLabelOn = False
;自定义等值线
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/10,25,50,100/)
res@cnFillColors =(/2,3,4,5,6/)
;色标显示
res@lbLabelBarOn = False
res@lbLabelAutoStride = True
; 垂直色标条
res@lbOrientation = "Vertical"
; pmLabelBarWidthF, pmLabelBarHeightF, lbLabelFontHeightF, and lbPerimOn.
; 调整色标图的大小
;res@lbPerimOn = False
res@pmLabelBarWidthF = 0.05
;res@pmLabelBarHeightF = 0.4
res@lbLabelFontHeightF = 0.012
res@lbPerimFillColor = 0
res@lbPerimThicknessF = 0.1
res@lbLabelStride = 4 ;色标标签间隔
;地图设置
res@gsnAddCyclic = False
res@mpFillOn = True
res@mpOutlineOn = False; Use outlines from shapefile
res@mpDataBaseVersion = "MediumRes"
res@mpDataSetName = "Earth..4"
res@mpLandFillColor = 0
res@mpInlandWaterFillColor = 0
res@mpOceanFillColor = 0
res@mpOutlineBoundarySets = "NoBoundaries"
;网格线设置
res@mpProjection = "CylindricalEquidistant" ;"Mercator" ;"CylindricalEquidistant" ; ;
res@mpLimitMode = "LatLon"; "Angles"
res@mpMinLatF = 22 ;30;17 ; y为纬度
res@mpMaxLatF = 26 ;40;55
res@mpMinLonF = 116 ;100;72
res@mpMaxLonF = 120 ;120;136
;四周边框
res@tmXBBorderOn = False
res@tmXTBorderOn = False
res@tmYLBorderOn = False
res@tmYRBorderOn = False
res@tmXBOn = False
res@tmXTOn = False
res@tmYLOn = False
res@tmYROn = False
res@mpGeophysicalLineThicknessF= 2.
res@mpNationalLineThicknessF= 2.
;res@mpLimitMode = "LatLon"
;res@mpLambertParallel1F = .001
;res@mpLambertParallel2F = 89.999
res@cnFillOn = True
res@cnLinesOn = False
res@cnLineLabelsOn = False
res@pmTickMarkDisplayMode = "never"
map_mask = gsn_csm_contour_map(wks,data,res)
opt = True
opt@areas_to_exclude = areas_of_interest
areas_to_fill = get_areas_of_interest(shapefile,"NAME",opt)
popt = True
popt@polytype = "polygon"
add_shapefile_primitives_by_name(wks,map_mask,shapefile,"NAME",areas_to_fill,popt)
lnres = True
lnres@gsLineColor = "black"
lnres@gsLineThicknessF = 2.0
line_mask = gsn_add_shapefile_polylines(wks, map_mask, shapefile, lnres)
draw(map_mask)
frame(wks)
end
|
-
想要的效果是相反的
|
窥视卡
雷达卡
发表于 2016-9-23 09:23:43
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主