爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 61818|回复: 53

[作图] 添加底图的另一种方法——shp文件

[复制链接]
发表于 2016-8-28 10:05:47 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
其中红色区域就是添加shapefile的内容参见官网
;   Example script to produce plots for a WRFreal-data run,
;   with the ARW coordinate dynamics option.
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/cnmap/cnmap.ncl"
load"$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
; TheWRF ARW input file.  
; Thisneeds to have a ".nc" appended, so just do it.
  a =addfile("/home/fc/datafile/wrfoutdata/wrfout_d01_2016-08-17_18:00:00","r")
; Wegenerate plots, but what kind do we prefer?
; type ="x11"
; type ="pdf"
  type = "ps"
; type ="ncgm"
  wks = gsn_open_wks(type,"surface")
; Set somebasic resources
  res = True
  ;res@MainTitle                   = "REAL-TIME WRF"
  pltres = True
  pltres@PanelPlot = True    ;这句话不写,图片经纬度会出问题,这句话找了好久
  mpres = True
;;;;;;;;;;;;;;;;;setfor map;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
mpres@pmTickMarkDisplayMode= "Always"  ;显示经纬度
mpres@mpLimitMode           = "LatLon"
mpres@mpMinLatF= 15
mpres@mpMaxLatF= 55
mpres@mpMinLonF= 72
mpres@mpMaxLonF= 136
mpres@mpFillOn                = True
mpres@mpOutlineOn             = False  ; Use outlines from shapefile
mpres@mpDataBaseVersion       = "MediumRes"
mpres@mpDataSetName           = "Earth..4"
mpres@mpAreaMaskingOn         = True
mpres@mpMaskAreaSpecifiers    =(/"China","Taiwan","Disputed area between India andChina","India:Arunachal Pradesh"/)
mpres@mpLandFillColor         = "white"
mpres@mpInlandWaterFillColor  = "white"
mpres@mpOceanFillColor        = "white"
mpres@mpOutlineBoundarySets   = "NoBoundaries"  ;or National
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Whattimes and how many time steps are in the data set?
  times =wrf_user_getvar(a,"times",-1) ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  do it = 0,ntimes-1,2             ; TIME LOOP
    print("Working on time: " +times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Firstget the variables we will need        
    slp =wrf_user_getvar(a,"slp",it)   ; slp
      wrf_smooth_2d( slp, 3 )            ; smooth slp
    tc =wrf_user_getvar(a,"tc",it)     ; 3D tc
    td =wrf_user_getvar(a,"td",it)     ; 3D td
    u  =wrf_user_getvar(a,"ua",it)     ; 3D U at mass points
    v  =wrf_user_getvar(a,"va",it)     ; 3D V at mass points
    td2 = wrf_user_getvar(a,"td2",it)  ; Td2 in C
    tc2 =wrf_user_getvar(a,"T2",it)    ; T2 in Kelvin
       tc2 = tc2-273.16                  ; T2 in C
    u10 =wrf_user_getvar(a,"U10",it)   ; u at 10 m, mass point
    v10 =wrf_user_getvar(a,"V10",it)   ; v at 10 m, mass point
;    tf2 = 1.8*tc2+32.                    ; Turn temperature intoFahrenheit
;      tc2@description = "SurfaceTemperature"
;      tc2@units = "~S~o~N~C"
;    td_f = 1.8*td2+32.                   ; Turn temperature intoFahrenheit
;    td2@description = "Surface Dew Point Temp"
;    td2@units = "~S~o~N~C"
;    u10 = u10*1.94386                    ; Turn wind into knots
;    v10 =v10*1.94386
;      u10@units = "kts"
;      v10@units = "kts"
;       dt=tc2-td2
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; Plotting options for T               
      opts = res                        
      opts@cnFillOn = True  
      opts@ContourParameters = (/ 5., 40., 5./)
      opts@cnFillColors =(/"cadetblue1","steelblue2","green","green4",\
                           "gold2","orange","red","red3","red4"/)
      opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
      contour_tc = wrf_contour(a,wks,tc2,opts)
      delete(opts)
    ; MAKE PLOTS                                       
      plot =wrf_map_overlays(a,wks,(/contour_tc/),pltres,mpres)
;================add china map ==================================;  
;****************************************************************************  
;               add  shapefiles  
;****************************************************************************
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;国界
shp1="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/bou1_4l.shp"  
  lnres1        = True     
lnres1@gsLineColor      ="black"
lnres1@gsLineThicknessF = 2.0           ; 2x thickness  
  shp_plot1     = gsn_add_shapefile_polylines(wks,plot,shp1,lnres1)  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;省级界线
shp2="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/bou2_4p.shp"  
  lnres2        = True   
  lnres2@gsLineColor      ="black"
lnres2@gsLineThicknessF = 1.5              
  shp_plot2     = gsn_add_shapefile_polylines(wks,plot,shp2,lnres2)  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;市级界线
shp3="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/diquJie_polyline.shp"  
  lnres3        = True     
  lnres3@gsLineColor      ="black"
lnres3@gsLineThicknessF = 1.0        
shp_plot3     =gsn_add_shapefile_polylines(wks,plot,shp3,lnres3)  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;河流
;shp4="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/hyd1_4l.shp"  
; lnres4        = True
;  lnres4@gsLineColor      = "blue"
  ;lnres4@gsLineThicknessF = 1.0   
; shp_plot4     =gsn_add_shapefile_polylines(wks,plot,shp4,lnres4)  
maximize_output(wks,False)      这句话不写,add shapefiles 就没意义!
;================add china map ==================================;  
  end do       ; END OF TIME LOOP
end

tem.png

评分

参与人数 2金钱 +25 贡献 +11 收起 理由
mofangbao + 15 + 5
kongfeng0824 + 10 + 6 很给力!

查看全部评分

密码修改失败请联系微信:mofangbao
 楼主| 发表于 2018-9-3 20:42:25 | 显示全部楼层
yybubble 发表于 2018-9-3 19:35
您好,我有一个流域shp文件,我想wrfout数据只在这个shp范围内显示,其他区域不填充颜色。
仿照ncl网站关 ...

mpres@mpFillBoundarysets = "Noboundaries" 看看这个可以不
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

发表于 2016-8-28 11:02:18 | 显示全部楼层
你可以使用这个边界http://www.gadm.org/country
效果图如下: QQ截图20160828110119.jpg
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

 楼主| 发表于 2016-8-28 11:07:17 | 显示全部楼层
andrewsoong 发表于 2016-8-28 11:02
你可以使用这个边界http://www.gadm.org/country
效果图如下:

哇,nice,nice
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2016-8-28 11:09:32 | 显示全部楼层
andrewsoong 发表于 2016-8-28 11:02
你可以使用这个边界http://www.gadm.org/country
效果图如下:

咦。网站打不开啊
密码修改失败请联系微信:mofangbao
发表于 2016-8-28 11:28:09 | 显示全部楼层
1649518749 发表于 2016-8-28 11:09
咦。网站打不开啊

最近G20开会,打不开正常。传给你吧

CHN_adm.zip

21.01 MB, 阅读权限: 10, 下载次数: 414, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao
发表于 2016-8-28 11:28:38 | 显示全部楼层

另外又加上了湖泊
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2016-8-28 11:30:31 | 显示全部楼层
andrewsoong 发表于 2016-8-28 11:28
另外又加上了湖泊

哦,谢谢了,这个好像其他论坛看到过,忘了是grads还是surfer
密码修改失败请联系微信:mofangbao
发表于 2016-8-28 11:56:24 | 显示全部楼层
谢谢分享   
密码修改失败请联系微信:mofangbao
发表于 2016-8-28 13:07:14 | 显示全部楼层
{:5_213:}{:5_213:}
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

发表于 2016-8-28 21:09:48 | 显示全部楼层
好美的图  赞
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表