地形图上加站点

当时年少春衫薄

undef("create_topo_map")
function create_topo_map(wks,minlat,maxlat,minlon,maxlon)
local topo_file, a, elev, res, cmap
begin

  topo_file =  "/home/fanweiwei/NCEP/dbjr/pictures/ht/shuju/ETOPO2v2c_f4.nc"
  a         = addfile(topo_file,"r")
  elev      = a->z({minlat:maxlat},{minlon:maxlon})
printVarSummary(a)
;  elev = where(elev.lt.-100.,elev@_FillValue,elev)
printVarSummary(elev)
elev!0      = "lat"
elev!1      = "lon"

elev&lat =  a->y({minlat:maxlat})
elev&lon =a->x({minlon:maxlon})
elev&lat@units = "degrees_north"  ; 赋予lat方向信息，这个是必须有的
elev&lon@units = "degrees_east"  ; 赋予lon方向信息，这个是必须有的
printVarSummary(elev)
  cmap = read_colormap_file("MPL_BrBG")


  res                    = True

  res@gsnMaximize        = True         ; maximize plot in frame

  res@gsnDraw            = False
  res@gsnFrame           = False

  res@cnFillOn           = True         ; turn on contour fill
; res@cnFillMode         = "MeshFill"   ; for faster draw
  res@cnFillPalette      = cmap
   res@cnFillMode         = "RasterFill"   ; for faster draw
  res@cnLinesOn          = False        ; turn off contour lines
  res@cnLineLabelsOn     = False        ; turn off line labels
  res@cnInfoLabelOn      = False        ; turn off info label
  res@lbBoxLinesOn       = False        ; turn off labelbar box lines
  res@lbTitleString      = "elevation (meters)"    ; add a labelbar title
  res@lbTitleFontHeightF = 0.015
  res@lbTitlePosition    = "Bottom"
  res@pmLabelBarOrthogonalPosF = 0.15

;---Pick "nice" contour levels
  mnmxint = nice_mnmxintvl( min(elev), max(elev), 18, False)
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF       = mnmxint(0)
  res@cnMaxLevelValF       = mnmxint(1)
  res@cnLevelSpacingF      = 50   ; Increase the number of levels
                                             ; by choosing a smaller spacing.
;---Zoom in on map
  res@mpMinLatF            = minlat
  res@mpMaxLatF            = maxlat
  res@mpMinLonF            = minlon
  res@mpMaxLonF            = maxlon
  res@mpCenterLonF         = (res@mpMinLonF + res@mpMaxLonF) / 2.

  res@mpDataBaseVersion    = "MediumRes"

  res@mpFillOn             = False
  res@mpOutlineOn          = True
  res@mpOutlineBoundarySets = "AllBoundaries"
  res@gsnAddCyclic         = False           ; don't add longitude cyclic point

  res@tiMainString          = "Rivers of Colorado"
  res@gsnLeftString         = ""
  res@gsnRightString        = ""
  res@pmTickMarkDisplayMode = "Always"
  res@pmTitleZone           = 4              ; move main title down a little
;  res@mpMaskAreaSpecifiers = (/"China:states"/)
;res@mpMaskAreaSpecifiers =(/"China:states","Taiwan"/)
res@mpDataSetName         = "Earth..4"   ; This new database contains
;---Create map and return it.
  plot = gsn_csm_contour_map(wks,elev,res)
  return(plot)
end


begin
  wks = gsn_open_wks("png","topo")          ; send graphics to PNG file

;---Lat/lon limits for Colorado
  minlat   = 25
  maxlat   = 40
  minlon   = 70
  maxlon   = 105

;---Create topo map of Colorado
  topo_map = create_topo_map(wks,minlat,maxlat,minlon,maxlon)
  station             = asciiread("/home/fanweiwei/NCEP/dbjr/pictures/ht/shuju/zhandian.txt",(/61,2/),"float")  

lat1  = station(:,0)
  lon1   = station(:,1)
  pmres                 = True                       

  pmres@minlat           = minlat
  pmres@maxlat           = maxlat
  pmres@minlon           = minlon
  pmres@maxlon           = maxlon


pmres@gsMarkerIndex    = 16
  pmres@gsMarkerSizeF    = 0.01
  pmres@gsMarkerColor    = "black"
   ;
d=gsn_add_polymarker(wks,topo_map,lon1,lat1,pmres)

;---Drawing the topo map will also draw shapefile outlines
  draw(topo_map)
  frame(wks)
end

当时年少春衫薄

以上是正确代码

当时年少春衫薄

谢谢帮助～

当时年少春衫薄

zly4814624 发表于 2016-10-31 08:33
Data either does not contain a valid latitude coordinate array or doesn't contain one at all
数据的 ...

您好我画出来了，非常感谢您的帮助，谢谢

wqywgq

学习了，这个好

zly4814624

当时年少春衫薄 发表于 2016-10-31 13:28
您好我画出来了，非常感谢您的帮助，谢谢

那太好了。请问具体是什么问题呢。

当时年少春衫薄

zly4814624 发表于 2016-10-31 14:28
那太好了。请问具体是什么问题呢。

您好，按照你的要求设置了一下地形数据的经纬度信息
a         = addfile(topo_file,"r")
  elev      = a->z({minlat:maxlat},{minlon:maxlon})
printVarSummary(a)
;  elev = where(elev.lt.-100.,elev@_FillValue,elev)
printVarSummary(elev)
elev!0      = "lat"
elev!1      = "lon"

elev&lat =  a->y({minlat:maxlat})
elev&lon =a->x({minlon:maxlon})
elev&lat@units = "degrees_north"  
elev&lon@units = "degrees_east"

平流层的萝卜

科罗拉多河咋还在西藏呢

zly4814624

当时年少春衫薄 发表于 2016-10-31 15:37
您好，按照你的要求设置了一下地形数据的经纬度信息
a         = addfile(topo_file,"r")
  elev      ...

太棒了。看来我蒙对了哈哈。

当时年少春衫薄

平流层的萝卜 发表于 2016-10-31 15:54
科罗拉多河咋还在西藏呢

哈哈，原来不在西藏的{:5_250:}