Re: US map with cities marked

From: Dennis Shea (shea AT XXXXXX)
Date: Mon Aug 16 2004 - 11:15:34 MDT

  • Next message: Dennis Shea: "Re:"

    >> >
    >> >I need to produce a conus map with some cities and their
    >> >surrounding points marked on the map. I have the city locations
    >> >in lat/long but the surrounding points are in RUC20 grid points.
    >> >[ie subscript indices]
    >> >
    >> >I was just wondering if someone out there has already done
    >> >something similar and can provide an NCL sample script.
    >> >
    >> >Thanks in advance.
    ===============================================================
    >
    >Here is a sample data file:
    >city % lat long x0 y0 x1 y1 x2 y2 x3 y3
    >--------------------------------------------------------
    > ABQ 5.82 35.05 106.62 111 86 111 87 112 86 112 87
    > ALB 6.16 42.75 73.80 253 134 253 135 254 134 254 135
    > ATL 8.98 33.75 84.38 212 78 212 79 213 78 213 79
    >
    >I want to mark each city location and its surrounding grid
    >points on the map with the city's name and the % value.
    >The xs and ys are RUC20 grid point locations. Will NCL
    >convert the grid point pair into lat/long?
    >

    Basically, the sequence is the following:

    [1] Read the ascii file. The 1st two lines [ie, rows]
        are header lines. The rest of the ascii file contains
        the data. The data are one column of ascii and the
        rest numeric.
        
        [a] use the "readAsciiTable" function located in
            contributed.ncl and the standard "asciiread" function
        [b] extract the pertinent info
        
    [2] Read appropriate RUC GRIB file that contians the
        appropriate lat/lon grid coordinates
        
    [3] Generate plot

        [a] generate map background
        [b] Use "gsn_polymarker" to generate a polymarker to indicate
            the city's location (filled circle used)
        [c] use "gsn_text" to generate the city's 3 alphabetic identifier
        [d] use "gsn_polymarker" to generate grid point locations
        
    Please look at the following:
        http://www.cgd.ucar.edu/csm/support/CSM_Graphics/polyg.shtml
        http://www.cgd.ucar.edu/csm/support/CSM_Graphics/text.shtml
        
    A sample script is attached.

    --
    One issue: are the given subscript indices 0-based like C/NCL
               or are they 1-based like fortran. If 1-based, you
               must subtract 1 from each subscript.
    

    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" begin ;------------------------------------------------------------- ; Read ascii file ;------------------------------------------------------------- dir = "./" fil = "city.ascii" ncol = 11 ; # col of numeric data nhead = 2 ; # header lines

    ; read all rows as string s = asciiread ( dir+fil, -1, "string") nrow = dimsizes(s) ;print (s) ;print (s(nhead:)) ; convert to character and create city string c = stringtochar(s(nhead:)) ; skip 1st two header lines city = chartostring(c(:,1:3)) print (city) ; read all numbers as float x = readAsciiTable( dir+fil, ncol, "float", nhead) printVarSummary (x) dimx = dimsizes(x) npts = dimx(0)

    ; For clarity assign members of x to separate variables pc = x(:,0) ; % latc = x(:,1) ; city lat lonc = -x(:,2) ; city lon (make W)

    ix0 = floattointeger(x(:,3)) iy0 = floattointeger(x(:,4)) ix1 = floattointeger(x(:,5)) iy1 = floattointeger(x(:,6)) ix2 = floattointeger(x(:,7)) iy2 = floattointeger(x(:,8)) ix3 = floattointeger(x(:,9)) iy3 = floattointeger(x(:,10))

    ;------------------------------------------------------------- ; Read RUC file ;-------------------------------------------------------------

    diri = "/cgd/cas/shea/celiaChen/" fili = "ruc2.bgrb.20020418.i12.f00.grb"

    f = addfile (diri+fili, "r") lat2d = f->gridlat_252 ; [gridx_252 | 225] x [gridy_252 | 301] lon2d = f->gridlon_252 ; [gridx_252 | 225] x [gridy_252 | 301] do n=0,npts-1 ; print surrounding pts print("n="+n+" "+city(n)+" "+lonc(n)+" "+latc(n)) print(" "+ix0(n)+" "+iy0(n)+" "+lon2d(iy0(n),ix0(n))+" "+lat2d(iy0(n),ix0(n)) ) print(" "+ix1(n)+" "+iy1(n)+" "+lon2d(iy1(n),ix1(n))+" "+lat2d(iy1(n),ix1(n)) ) print(" "+ix2(n)+" "+iy2(n)+" "+lon2d(iy2(n),ix2(n))+" "+lat2d(iy2(n),ix2(n)) ) print(" "+ix3(n)+" "+iy3(n)+" "+lon2d(iy3(n),ix3(n))+" "+lat2d(iy3(n),ix3(n)) ) end do

    ;------------------------------------------------------------- ; plot ;------------------------------------------------------------- wks = gsn_open_wks("ps","city") ; open a ps file rmap = True rmap@gsnMaximize = True ; maximize size rmap@gsnFrame = False ; manually advance

    rmap@mpProjection = "LambertConformal" rmap@mpLambertParallel1F = lat2d@Latin1 rmap@mpLambertParallel2F = lat2d@Latin2 rmap@mpLambertMeridianF = lat2d@Lov

    rmap@mpLimitMode = "Corners" ; choose range of map rmap@mpLeftCornerLatF = lat2d@corners(0) rmap@mpLeftCornerLonF = lon2d@corners(0) rmap@mpRightCornerLatF = lat2d@corners(2) rmap@mpRightCornerLonF = lon2d@corners(2)

    rmap@pmTickMarkDisplayMode = "Always" ; lat/lon tick marks ;res@tmXTOn = False ; turn off top tick marks ;res@tmYROn = False ; turn off right tick marks

    rmap@mpFillOn = False ; turn off map fill rmap@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last rmap@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries

    map = gsn_csm_map(wks,rmap) ; create a default plot

    res_poly = True res_poly@gsMarkerIndex = 16 ; polymarker style res_poly@gsMarkerSizeF = 6. ; polymarker size res_poly@gsMarkerColor ="red" gsn_polymarker(wks,map,lonc,latc,res_poly) res_text = True res_text@txJust = "BottomCenter" ;"CenterCenter" ; "CenterLeft","CenterRight" res_text@txFontHeightF = 0.0075 ; size gsn_text (wks,map,city,lonc,latc+0.5,res_text) ; offset lat

    ; plot local RUC grid pts

    res_grid = True res_grid@gsMarkerIndex = 2 ; polymarker style res_grid@gsMarkerSizeF = 6. ; polymarker size res_grid@gsMarkerColor ="green"

    do n=0,npts-1 gsn_polymarker(wks,map,lon2d(iy0(n),ix0(n)),lat2d(iy0(n),ix0(n)),res_grid) gsn_polymarker(wks,map,lon2d(iy1(n),ix1(n)),lat2d(iy1(n),ix1(n)),res_grid) gsn_polymarker(wks,map,lon2d(iy2(n),ix2(n)),lat2d(iy2(n),ix2(n)),res_grid) gsn_polymarker(wks,map,lon2d(iy3(n),ix3(n)),lat2d(iy3(n),ix3(n)),res_grid) end do

    frame(wks) end

    _______________________________________________ ncl-talk mailing list ncl-talk@ucar.edu http://mailman.ucar.edu/mailman/listinfo/ncl-talk



    This archive was generated by hypermail 2b29 : Mon Aug 16 2004 - 13:14:59 MDT