province-basemap.gs文件_provincebasemap的用法grads

province-basemap.gs文件（见下面的帖子的附件，下载后拷贝到pcgrads/lib下）和32个省份的边界文件（见下面的帖子的附件）以及1个重新修改了的全国边界文件。

实现这个需要的文件除了本省的数据文件外，还要有边界文件和地图文件（各省的地图可从aq9807制作的全国各省地图，见

http://www.lasg.ac.cn/cgi-bin/fo ... um=3&topic=4556）下载，然后只需要按照的上面原则重新把各省的改名字修改一下就可以了。

注意：各省的文件名字均采用汉语拼音的全部（如xizang对应西藏；guangdong对应广东；不过只有shanxi对应山西，shannxi对应陕西是例外，广东边界文件是guangdong_out.txt。地图文件是guangdong，对应陕西的边界文件是shannxi_out.txt，陕西的地图边界文件是shannxi）。

需要在你的pcgrads目录下新建一个out文件夹，然后把边界文件拷贝到该文件下。

* Written 12/2005 by tibet tibet-weather@163.com
****************************************************************************************
function main(args)
dxy=0.01
* There are defaults for the colors and resolution
* User must specify which mask
***********************   USAGE   *********************************************
usage1="Usage: p-map  provincname varname < color > "
usage2="color: colors of line (Default=0)"
*******************************************************************************
say args
  if (args='')
  say usage1 ;say usage2
  return
  endif
  map = subwrd(args,1)
* map is provincename
  varname = subwrd(args,2)
  color   = subwrd(args,3)
if (color = '')   ; color = 0  ; endif
file = 'c:/PCGrADS/out/'map'_out.txt';say "read file name="file
'set mpdset 'map  
  'd 'varname

* Make sure there's a plot already drawn
  'q gxinfo'
  line5 = sublin(result,5)
  line6 = sublin(result,6)
  xaxis = subwrd(line5,3)
  yaxis = subwrd(line5,6)
  proj  = subwrd(line6,3)

'q config'
  line = sublin(result,1)
  word = subwrd(line,2)
  version = substr(word,2,3)
  if (version >= 1.8)
    newgrads = 1
  else
    newgrads = 0
  endif
  
* See if map projection will be supported
  if (newgrads = 0)
    if (proj != 1 & proj != 2)
      say 'Error: Only scaled or latlon projections are supported with GrADS v'version
      return
    endif
  endif
* Clip image accordingly
  'q gxinfo'
  line3 = sublin(result,3)
  line4 = sublin(result,4)
  x1 = subwrd(line3,4)
  x2 = subwrd(line3,6)
  y1 = subwrd(line4,4)
  y2 = subwrd(line4,6)
  'set clip 'x1+dxy' 'x2-dxy' 'y1+dxy' 'y2-dxy

* Read the first record from the polygon file
  result = read(file)
  rc = sublin(result,1)
  rc = subwrd(rc,1)
 
  if (rc!=0)
    say 'Error reading 'file
    return
  endif
  nwcmd = sublin(result,2)
flag = 1
  while (flag)
    ignore = 0
    wcmd = nwcmd
    while(1)
      result = read(file)
      rc = sublin(result,1)
      rc = subwrd(rc,1)
      if (rc!=0)
        flag = 0
        break
      else
        nwcmd = sublin(result,2)
        if (subwrd(nwcmd,5) != 'draw')
          wcmd = wcmd % nwcmd
        else
          break
        endif
      endif
    endwhile

*   Get the lat/lon range of the current dimension environment
    'q dims'
    line1 = sublin(result,2)
    line2 = sublin(result,3)
    minlon = subwrd(line1,6)
    maxlon = subwrd(line1,8)
    minlat = subwrd(line2,6)
    maxlat = subwrd(line2,8)
*   The range of the polygon is specified in the first four words of the record
    minwx = subwrd(wcmd,1)
    maxwx = subwrd(wcmd,2)
    minwy = subwrd(wcmd,3)
    maxwy = subwrd(wcmd,4)

*   If the polygon is outside the current dimension, ignore it
    if (minwx >= maxlon) ; ignore = 1 ; endif
    if (maxwx <= minlon) ; ignore = 1 ; endif
    if (minwy >= maxlat) ; ignore = 1 ; endif
    if (maxwy <= minlat) ; ignore = 1 ; endif

    if (!ignore)   
      count = 6
      nvert = 1
      if (newgrads) 
        cmd = 'draw mappoly '
       else
        cmd = 'draw polyf '
      endif
      while (1)
        countx = count
        county = count + 1
        wx = subwrd(wcmd,countx)
        wy = subwrd(wcmd,county)
        if ((wx = '') | (wy = ''))
          break
        endif

*       Convert world coordinates to screen coordinates if necessary
        if (newgrads) 
          sx = wx
          sy = wy
       else
          'q w2xy 'wx' 'wy
          sx = subwrd(result,3)
          sy = subwrd(result,6)
        endif

*       Append the coordinates to the draw command
        cmd = cmd%sx' 'sy' '
       count = count + 2
      endwhile  

*     Draw the polygon
      'set line 'color
       cmd
    endif
  endwhile
'draw map'
  rc = close(file)
  return

cnbasemap.gs
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* chinamap.gs
*
* Download all polygon files from ftp://grads.iges.org/grads/scripts/
*     
* Written 12/2000 by Jennifer M. Adams, jma@cola.iges.org
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
function main(args)
'set mpdset cnworld'

dxy=0.01
* There are defaults for the colors and resolution
* User must specify which mask
***********************   USAGE   *********************************************
usage1="Usage: chinamap varname < color > "
usage2="color: colors of line (Default=0)"
*******************************************************************************
say args
  if (args='')
   say usage1 ;say usage2
  return
  endif
  varname = subwrd(args,1)
color   = subwrd(args,2)
if (color = '')   ; color = 0  ; endif  
  'd 'varname

  file = 'C:/PCGrADS/out/cn_out.txt'

* Make sure there's a plot already drawn
  'q gxinfo'
  line5 = sublin(result,5)
  line6 = sublin(result,6)
  xaxis = subwrd(line5,3)
  yaxis = subwrd(line5,6)
  proj  = subwrd(line6,3)

'q config'
  line = sublin(result,1)
  word = subwrd(line,2)
  version = substr(word,2,3)
  if (version >= 1.8)
    newgrads = 1
  else
    newgrads = 0
  endif
  
* See if map projection will be supported
  if (newgrads = 0)
    if (proj != 1 & proj != 2)
      say 'Error: Only scaled or latlon projections are supported with GrADS v'version
      return
    endif
  endif
* Clip image accordingly
  'q gxinfo'
  line3 = sublin(result,3)
  line4 = sublin(result,4)
  x1 = subwrd(line3,4)
  x2 = subwrd(line3,6)
  y1 = subwrd(line4,4)
  y2 = subwrd(line4,6)
  'set clip 'x1+dxy' 'x2-dxy' 'y1+dxy' 'y2-dxy

* Read the first record from the polygon file
  result = read(file)
  rc = sublin(result,1)
  rc = subwrd(rc,1)
  if (rc!=0)
    say 'Error reading 'file
    return
  endif
  nwcmd = sublin(result,2)
*say nwcmd
flag = 1
  while (flag)
    ignore = 0
    wcmd = nwcmd
    while(1)
      result = read(file)
      rc = sublin(result,1)
      rc = subwrd(rc,1)
      if (rc!=0)
        flag = 0
        break
      else
        nwcmd = sublin(result,2)
        if (subwrd(nwcmd,5) != 'draw')
          wcmd = wcmd % nwcmd
        else
          break
        endif
      endif
    endwhile

*   Get the lat/lon range of the current dimension environment
    'q dims'
    line1 = sublin(result,2)
    line2 = sublin(result,3)
    minlon = subwrd(line1,6)
    maxlon = subwrd(line1,8)
    minlat = subwrd(line2,6)
    maxlat = subwrd(line2,8)

*   The range of the polygon is specified in the first four words of the record
    minwx = subwrd(wcmd,1)
    maxwx = subwrd(wcmd,2)
    minwy = subwrd(wcmd,3)
    maxwy = subwrd(wcmd,4)

*   If the polygon is outside the current dimension, ignore it
    if (minwx >= maxlon) ; ignore = 1 ; endif
    if (maxwx <= minlon) ; ignore = 1 ; endif
    if (minwy >= maxlat) ; ignore = 1 ; endif
    if (maxwy <= minlat) ; ignore = 1 ; endif

    if (!ignore)   
      count = 7
      nvert = 1
      if (newgrads) 
        cmd = 'draw mappoly '
      else
        cmd = 'draw polyf '  
      endif
      while (1)
        countx = count
        county = count + 1
        wx = subwrd(wcmd,countx)
        wy = subwrd(wcmd,county)
        if ((wx = '') | (wy = ''))
          break
        endif

*       Convert world coordinates to screen coordinates if necessary
        if (newgrads) 
          sx = wx
          sy = wy
        else
          'q w2xy 'wx' 'wy
          sx = subwrd(result,3)
          sy = subwrd(result,6)
        endif

*       Append the coordinates to the draw command
        cmd = cmd%sx' 'sy' '
        count = count + 2
      endwhile  

*     Draw the polygon
      'set line 'color
      cmd
    endif
  endwhile
'draw map'
  rc = close(file)
  return