Hi M.,
If you're interested specifically in plotting the METAR-reporting
stations, you're in luck - I have a script to do this (attached).
The script takes it's input from the list of all the METAR-reporting
surface stations in the world which is maintained at Greg Thompson's
excellent site:
http://www.rap.ucar.edu/weather/surface/stations.txt
Once you've downloaded that, you can run my script to plot all the
METAR-reporting stations for a particular region. You can adjust the
style and labeling to meet your needs.
At least this gives you a starting point.
Good luck!
Jonathan
marziyeh dadizadeh wrote:
> Hi,
> I want plot surface meteorological stations with their name on the
> regional map.
> Please guide me how can do it???
> Thanks Alot
> M.Dadizadseh
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
;********************************************************
; METAR_general.ncl
;
; The purpose of this script is to plot all METAR reporting
; stations for a specified region.
;
; This script takes general input from the stations.txt file
; available at from the UCAR-RAP Weather Page at:
; http://www.rap.ucar.edu/weather/surface/stations.txt
;
; Author: Jonathan Vigh, Colorado State University
;
; Program input:
; The program expects the following input parameters to be passed as command line arguments when the script is called:
; map_name desired filename of output map
; min_lat minimum latitude of plotting domain (positive north)
; max_lat maximum latitude of plotting domain
; min_lon minimum longitude of plotting domain (positive east)
; max_lon maximum longitude of plotting domain
;
; These can be easily set using the run_mapper script. If the command line arguments are not set, some defaults will be used (which makes a map of Colorado).
;
; Program History:
; circa 2004 Program created (originally read in a hand-edited file: stations.inp)
; 09/14/2007: Revised program input section to fix issues with reading directly from the stations.txt file (had been buggy before)
; 09/17/2007: Revised program to accept command line inputs - the NCL code is called from a script called run_mapper.
; 09/18/2007: Program now utilizes a much more efficient way of selecting the stations in the plotting domain.
; 09/19/2009: Fixed logic for indexing regional stations to account for a change in how NCL handles the _FillValue of logical variables.
;
; Suggestions for improvements:
; - *** Figure out how to take care of the overlapping labels problem (see the text example 10 code for ideas)
; - Put some place names on the map, like for island groups.
; - Have an option that can choose to plot all the stations in a specific country or state (use the two-letter code)
; - Add a few default mapping options? Or plot several types of maps (i.e. high resolution close-up vs. regular resolution?
; - Control font height of labels based on the size of the map (tie to span of lat/lon?)
; - Read in topography from a database, so we could get some nice terrain maps.
; - Have the program query the stations - only plot those which have given a recent report
;
;********************************************************
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/csm/shea_util.ncl"
;********************************************************
begin
msg_val = -999
nullChar = inttochar(0) ; could use this to test a single char to see if it's null
; Check to see if the command line arguments are present
if isvar("min_lat") then
use_RANGS = False ; set to True to enable high res coastlines and inland water boundaries, set to False to use NCL's Earth..4 dataset
minlat = stringtofloat(min_lat)
maxlat = stringtofloat(max_lat)
minlon = stringtofloat(min_lon)
maxlon = stringtofloat(max_lon)
print("map_name = "+map_name)
print("minlat = "+minlat)
print("maxlat = "+maxlat)
print("minlon = "+minlon)
print("maxlon = "+maxlon)
else ; otherwise set some default plotting parameters (if they are not passed in as command line arguements)
map_name = "Colorado"
use_RANGS = False ; set to True to enable high res coastlines and inland water boundaries, set to False to use NCL's Earth..4 dataset
minlat = 36.8
maxlat = 41.2
minlon = -109.2
maxlon = -101.8
end if
;********************************
; get data
;********************************
tStr = asciiread("stations.txt", -1,"string") ; Info includes type so cases can be weeded out
tChr = stringtochar(tStr)
nrows = dimsizes(tStr)
temptStr = new(nrows,"string")
; Read in the entire file, put all the lines which actually contain station data into a new array of strings
nglobalstations = 0
do irow = 0, nrows-1
FIRST = tChr(irow,0) ; next read the first letter in as a string
if ( (FIRST .ne. nullChar) .and. (FIRST .ne. "!") ) then ; only read the line if it's not a comment (!) or blank (test to see if the first character is the nullChar that was created using inttochar(0))
; if ( (strlen(tStr(irow)) .gt. 0) .and. (FIRST .ne. "!") ) then ; only read the line if it's not blank (string size greater than 0) - this is an alternative to checking the first character against the null chararacter
; print(tStr(irow))
; print(" The first character is: "+FIRST)
; print(" We know this line is not a comment or blank, continuing to read . . .")
DASH = tChr(irow,25:25) ; read the line to see if there is a dash in the 25th column (indicating a country or state identifier)
; print(" The 25th character is: "+DASH)
if (DASH .ne. "-") then ; don't read the line if it's a country/state identifier
if ( (tStr(irow) .ne. "CD STATION ICAO IATA SYNOP LAT LONG ELEV M N V U A C") .and. \
(tStr(irow) .ne. "CD STATION ICAO IATA SYNOP LAT LONG ELEV M N V U A X") ) then
; print(" It looks like we have a valid station line, entering data into the arrays . . .")
temptStr(nglobalstations) = tStr(irow)
nglobalstations = nglobalstations + 1
end if
end if
end if
end do
print(nrows+" rows were read from the data file.")
print(nglobalstations+" stations were read in for the entire world.")
; Make new array which only contains valid the valid station lines
; size = num(.not.ismissing(temptStr))
; print(" The number of nonmissing station lines is: "+size)
newtStr = temptStr(0:nglobalstations-1)
newtChr = stringtochar(newtStr)
; Read in the individual station values (note that if we are doing any conversions from strings to numerical type, we should check each individual value)
CD = chartostring(newtChr(:,0:1)) ; read the 2-letter state or province abbreviation (many international stations do not have this)
NAME = chartostring(newtChr(:,3:18))
ICAO = chartostring(newtChr(:,20:23))
IATA = chartostring(newtChr(:,26:28))
WMO_string = chartostring(newtChr(:,32:36))
LAT_HR_string = chartostring(newtChr(:,39:40))
LAT_MIN_string = chartostring(newtChr(:,42:43))
LON_HR_string = chartostring(newtChr(:,47:49))
LON_MIN_string = chartostring(newtChr(:,51:52))
ELEV_string = chartostring(newtChr(:,55:58))
METAR = chartostring(newtChr(:,62:62)) ; this flag specifies whether this is a METAR-reporting station
NEXRAD = chartostring(newtChr(:,65:65)) ; this flag specifies whether this is a NEXRAD (WSR-88D) Radar site
AVIATION = chartostring(newtChr(:,68:68)) ; this flag specifies some Aviation-specific stuff (V=AIRMET/SIGMET end point, A=ARTCC T=TAF U=T+V)
UPPERAIR = chartostring(newtChr(:,71:71)) ; this flag specifies an upper air (rawinsonde=X) or wind profiler (W) site
TYPE = chartostring(newtChr(:,74:74)) ; this flag specifies the station type: A = Auto (A=ASOS, W=AWOS, M=Meso, H=Human, G=Augmented) (H/G not yet impl.)
OFFICE = chartostring(newtChr(:,77:77)) ; this flag specifies the office type (F=WFO/R=RFC/C=NCEP Center)
NorS = chartostring(newtChr(:,44:44))
EorW = chartostring(newtChr(:,53:53))
WMO = new(nglobalstations,"integer",msg_val)
WMO = msg_val
LAT_HR = new(nglobalstations,"float",msg_val)
LAT_HR = msg_val
LAT_MIN = new(nglobalstations,"float",msg_val)
LAT_MIN = msg_val
LON_HR = new(nglobalstations,"float",msg_val)
LON_HR = msg_val
LON_MIN = new(nglobalstations,"float",msg_val)
LON_MIN = msg_val
ELEV = new(nglobalstations,"integer",msg_val)
ELEV = msg_val
lat = new(nglobalstations,"float",msg_val)
lat = msg_val
lon = new(nglobalstations,"float",msg_val)
lon = msg_val
inside = new(nglobalstations,"logical")
inside = False
IDEN = ICAO ; by default, use the ICAO as the station identifier, if this happens to be blank (rare), use the IATA (see code in loop below)
do i=0,nglobalstations-1
if ( (ICAO(i) .eq. " ") .and. (IATA(i) .ne. " ") ) then
IDEN(i) = IATA(i)
end if
if (WMO_string(i) .ne. " ") then
WMO(i) = stringtointeger(WMO_string(i))
end if
if (LAT_HR_string(i) .ne. " ") then
LAT_HR(i) = stringtofloat(LAT_HR_string(i))
end if
if (LAT_MIN_string(i) .ne. " ") then
LAT_MIN(i) = stringtofloat(LAT_MIN_string(i))
end if
if (LON_HR_string(i) .ne. " ") then
LON_HR(i) = stringtofloat(LON_HR_string(i))
end if
if (LON_MIN_string(i) .ne. " ") then
LON_MIN(i) = stringtofloat(LON_MIN_string(i))
end if
if (ELEV_string(i) .ne. " ") then
ELEV(i) = stringtointeger(ELEV_string(i))
end if
latMP = 1
if (NorS(i) .eq. "S") then
latMP = -1
end if
lonMP = 1
if (EorW(i) .eq. "W") then
lonMP = -1
end if
lat(i) = latMP*(LAT_HR(i) + LAT_MIN(i)/60.0)
lon(i) = lonMP*(LON_HR(i) + LON_MIN(i)/60.0)
; set a flag if the lat and lon are outside the desired map area (this logic is probably flawed for near the Dateline?)
if ( (lat(i) .ge. minlat) .and. (lat(i) .le. maxlat) .and. (lon(i) .ge. minlon) .and. (lon(i) .le. maxlon) ) then
inside(i) = True
end if
end do
;**************************************************************************************
; Mask section (get only METAR-reporting sites within the specified lat and lon ranges)
;**************************************************************************************
; newlat = lat
; newlon = lon
; newIDEN = IDEN
; newELEV = ELEV
; newlat@_FillValue = -999
; newlon@_FillValue = -999
; newIDEN@_FillValue = -999
; newELEV@_FILLValue = -999
; newlat = mask(lat,(inside .and. (METAR.eq."X")),True)
; newlon = mask(lon,(inside .and. (METAR.eq."X")),True)
; newIDEN = mask(IDEN,(inside .and. (METAR.eq."X")),True)
; newELEV = mask(ELEV,(inside .and. (METAR.eq."X")),True)
; Now pick off the indices of stations which are inside the plotting domain
localindex = ind( .not.ismissing(lat) .and. inside )
nlocalstations = dimsizes(localindex)
;***** Define some NICE color maps to use.
colors_6 = (/"White","Black","Black","MediumPurple1","MediumPurple3","Blue1",\
"CadetBlue3","Aquamarine2","SeaGreen2","LawnGreen","Green4", \
"GreenYellow","Gold","Tan1","Sienna1","Tomato","VioletRed1", \
"Yellow","DarkGreen","Grey37","Red","Orange","GoldenRod1","DarkOrange","SteelBlue1","SlateBlue1", \
"LightSlateBlue","DarkSeaGreen","Magenta","DodgerBlue","Moccasin","Blue"/)
;********************************
; create plot
;********************************
wks = gsn_open_wks("eps","map_"+map_name)
gsn_define_colormap(wks,colors_6)
res = True
res_at_tiMainString = "METAR-reporting stations for " + map_name
; res_at_gsnCenterString = "METAR-reporting stations"
res_at_gsnMaximize = True
; res_at_gsnPaperOrientation = "Landscape"
res_at_gsnDraw = False ; so we can add poly stuff
res_at_gsnFrame = False ; do not advance frame
if (use_RANGS) then
res_at_mpDataBaseVersion = "HighRes" ; set to HighRes to get nice coastlines
res_at_mpDataResolution = "Finest"
else
res_at_mpDataBaseVersion = "Ncarg4_1" ; set to HighRes to get nice coastlines
res_at_mpDataSetName = "Earth..4" ; we want county lines drawn (set to 3 for climate districts)
end if
res_at_mpProjection = "Mercator"
; res_at_mpProjection = "LambertConformal"
; res_at_mpLambertParallel1F = 20.0
; res_at_mpLambertParallel1F = 40.0
; res_at_mpLambertMeridianF = -60.0
res_at_mpLimitMode = "LatLon"
res_at_mpMinLatF = minlat
res_at_mpMaxLatF = maxlat
res_at_mpMinLonF = minlon
res_at_mpMaxLonF = maxlon
res_at_mpFillOn = True ; False to turn off gray continents
res_at_mpOutlineOn = True ; turn on continental outline
res_at_mpOutlineBoundarySets = "AllBoundaries"
res_at_mpLandFillColor = "GoldenRod1"
res_at_mpInlandWaterFillColor = "PaleTurquoise3"
res_at_mpOceanFillColor = "PaleTurquoise3"
res_at_mpGeophysicalLineColor = "Black"
res_at_mpGeophysicalLineThicknessF = 0.5
res_at_mpUSStateLineColor = "Red"
res_at_mpUSStateLineThicknessF = 2.5
res_at_mpNationalLineColor = "Black"
res_at_mpNationalLineThicknessF = 2.5
res_at_mpGridAndLimbOn = "True"
res_at_mpGridAndLimbDrawOrder = "Draw"
if (res_at_mpDataBaseVersion .ne. "HighRes") then
res_at_mpGridMaskMode = "MaskLand"
end if
res_at_mpGridSpacingF = 1.0
res_at_mpGridLineColor = "Blue"
res_at_tmXBLabelFontHeightF = 0.005
res_at_tmXBMajorLengthF = -0.001
res_at_pmTickMarkDisplayMode = "Always"
print("Now plotting the map")
;***************************************
; plot base map *
;***************************************
plot = gsn_csm_map(wks,res) ; draw one of eight map projections
; getvalues plot
; "mpAreaNames" : anames
; end getvalues
;
; print(anames)
;***************************************
; Draw METAR stations as polymarkers *
;***************************************
res_mark = True
res_mark_at_gsMarkerIndex = 1 ; polymarker style
res_mark_at_gsMarkerSizeF = 0.018 ; polymarker size
res_mark_at_gsMarkerColor = "Red" ; change marker color
dum = gsn_add_polymarker(wks,plot,lon(localindex),lat(localindex),res_mark)
print("Now plotting text labels")
;***************************
; Plot some text labels *
;***************************
txres = True
; Now draw the 4-letter ICAO station identifiers
txres_at_txFontHeightF = 0.008
txres_at_txJust = "TopCenter"
txres_at_txPerimOn = True
txres_at_txPerimColor = "Black"
txres_at_txPerimThicknessF = 0.3
txres_at_txPerimSpaceF = 0.4
txres_at_txBackgroundFillColor = "White"
; Calculate label offset as a function of map domain size
dy = (maxlat - minlat)/125.
; print(dy)
; Put the station labels onto the map as text strings
dumb = new(nlocalstations,graphic)
do j=1,nlocalstations-1
i = localindex(j) ; return the value of the global station number from our big array
; gsn_text(wks,plot,IDEN(i)+" "+ELEV(i),lon(i),lat(i)-dy,txres)
dumb(j) = gsn_add_text(wks,plot,IDEN(i)+" "+ELEV(i)+" m~C~"+NAME(i),lon(i),lat(i)-dy,txres) ; this is really slow
end do
print(newtStr(localindex))
; txres_at_txFontHeightF = 0.01
; txres_at_txPerimThicknessF = 1.5
; txres_at_txPerimSpaceF = 0.5
; txres_at_txConstantSpacingF = 1.0
; txres_at_txFontAspectF = 1.5
; txres_at_txFontThicknessF = 3.0
; gsn_text_ndc(wks,"Azore Islands",0.5,0.72,txres)
draw(plot)
frame(wks)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Jul 19 2009 - 11:52:18 MDT
This archive was generated by hypermail 2.2.0 : Thu Jul 23 2009 - 08:02:42 MDT