BASEMAP User defined mask .asc file
Frank Bandle
frank at BANDLE.COM
Fri Oct 14 03:46:22 EDT 2005
Hi Luigi
here is program in written in Fortran by Dr. Rupa Kumar Kolli .
I guess that is what yor are looking for.
see description and program above
regards
Frank
PROGRAM GRADS_MAP
c
c Program to prepare customized background map for GrADS from a text
c file containing line segments in (longitude,latitude) coordinates.
c
c By Rupa Kumar Kolli, Indian Institute of Tropical Meteorology
c kolli at tropmet.res.in
c created: December 10, 1997
c modified: March 12, 1998
c PC Fortran support added: January 26, 1999
c minor adjustments: March 20, 2002
c ifort compile option indicated: June 3, 2005
c Please do contact me for comments/feedback etc. on this program.
c
c This program converts map files in ASCII into GrADS readable
format.
c The resolution of the map background, of course, depends on the
c basic data that you can get. The format of the GrADS map data set
c prepared by this program conforms to that used in 'lowres'.
While it
c should work for any kind of resolution, it may be a bit slower
than
c the standard 'hires' for huge files.
c
c Input File Format:
c The input file can contain any number of line segments.
c Each segment starts with a header record containing the
c the number of pairs of coordinates and associated map type
c (see format statement 1)
c (the map type is used to selectively display map lines using
c the 'set mpt ...' command)
c The coordinates should be listed from the subsequent record,
c as longitude,latitude pair, one pair per line
c (see format statement 2).
c
c Example:
c
c 5 2
c 209.830 -76.580
c 211.500 -76.750
c 212.830 -76.420
c 211.830 -76.170
c 209.830 -76.580
c
c VERY IMPORTANT: NOTE THE FOLLOWING !!!
c
c 1. Longitudes should be in the range 0 to 360 (No negative
values !).
c
c 2. Latitudes should be in the range -90 (South Pole) to 90
(North Pole).
c
c 3. DON'T CROSS THE GREENWICH MERIDIAN (0.0 or 360.0) IN ANY
SEGMENT.
c For lines across this meridian, split them into two parts, with
c a common point between them. Use the value 360.0 for Greenwich
c meridian if you are moving towards it from west to east, and
the
c value 0.0 if you are moving from east to west. That is, you
c cannot go from 359.0 to 0.0; this will wrap the line across the
c whole globe !
c
c 4. The maximum number of points (pairs of coordinates) in a single
c segment is 255.
c
c 5. Don't use GrADS standard file names (lowres, mres, hires,
nam) for
c the output file (e.g., rupres). Preserve them for standard use.
c
c 6. After generating the GrADS map background file, copy it to a
path
c where GrADS looks for its support files (generally the path
defined
c in $GADDIR)
c
c 7. Finally, start GrADS and give the command 'set mpdset
filename':
c e.g., set mpdset rupres. Note that you can keep multiple map
c data sets open in GrADS 1.7; thus, you can create separate map
c data files and selectively use subsets of the files. For this
c purpose, use the command like 'set mpdset map1 map2', whereupon
c GrADS uses both the files for drawing the map.
c
c 8. This program was initially developed on a HP 9000/720
machine with
c HP-UX 9.05 operating system. I have successfully tested it
in dos/
c Windows as well as in RedHat Linux (g77). It should work on
other
c platforms as well. For using this program on a PC (dos/
windows/
c linux), a slight change is necessary, as indicated in the
comments
c preceding the output statements (marked by --->) in the
program.
c
c 9. The GrADS map file generated by this program is platform-
independent,
c and it can be used across platforms (I didn't check Mac,
though).
c
c 10. If you have successfully generated a specialized GrADS map
that may
c be of interest to other GrADS fans, please consider offering
it for
c sharing, possibly by making a suitable announcement in GrADS
mailing
c list.
c
c 11. If you are using recent versions of Intel FORTRAN, you will
have to
c compile this program with 'assume byterecl' option.
c i.e., ifort -assume byterecl gradsmap.f
c
c ENJOY HAVING YOUR OWN CUSTOMIZED GrADS MAP !!!
byte ihead(3),ilon(4),ilat(4)
integer lon,lat
equivalence (lon,ilon), (lat,ilat)
c
c Open the input and output files
c
open(10,file='spain.txt',
> status='old')
open(20,file='spain.map',
> form='unformatted',access='direct',recl=3)
c
c Initialize the first header-record byte, ihead(1): fixed to be 1
c A value of 2 can be used to skip unwanted sections of the file,
c (mres and hires were created like this), but this facility
c is not implemented in this program. So don't change this value
c unless you want to implement skipping (could be fairly
complicated !)
c
ihead(1)=1
irec=0
nsegs=0
c
c Read the number of pairs and map type for this segment.
c Leave the loop if end of file is reached.
c
99 read(10,1,end=100)npairs,ihead(2)
1 format(2i6)
print *,npairs,ihead(2)
nsegs=nsegs+1
if (npairs.gt.255) then
print *,'Number of pairs in Segment No.',nsegs,' exceeds 255 !'
stop 'Job Aborted'
endif
c
c Write the header record (3 bytes)
c The third byte indicates the number of pairs
c of values following (maximum=255)
c
ihead(3)=npairs
irec=irec+1
write(20,rec=irec)ihead
c
c Loop over the current segment
c
do 10 i=1,npairs
c
c Read one pair of coordinates (Longitude,Latitude)
c
read(10,2)alon,alat
2 format(2f7.3)
c
c Convert the values to integers, in GrADS readable format
c
lon=int(alon*10000.)+0.5
lat=int((alat+90.)*10000.)+0.5
c
c For use with PCs running DOS/Windows/Linux and Alpha workstations:
c Write the first three-bytes of lon and lat on the output file,
c in reverse order.
c Comment out the statements meant for Other Workstations.
c
c For use with Other Workstations (e.g., HP-UX, IRIX, SOLARIS):
c Write the last three-bytes of lon and lat on the output file.
c Comment out the statements meant for PCs and Alphas.
c
c Note that ilon and ilat are byte equivalents of lon and lat.
c
irec=irec+1
c ---> For PC (dos/windows/linux) and Alpha workstations:
write(20,rec=irec)(ilon(k),k=3,1,-1)
c ---> For Other Workstations:
c write(20,rec=irec)(ilon(k),k=2,4)
irec=irec+1
c ---> For PC (dos/windows/linux) and Alpha workstations:
write(20,rec=irec)(ilat(k),k=3,1,-1)
c ---> For Other Workstations:
c write(20,rec=irec)(ilat(k),k=2,4)
10 continue
c
c Go back to read next line segment
c
goto 99
100 print *,'Total Number of Segments Processed=',nsegs
close (10)
close (20)
stop 'Normal Termination'
end
Am 29.09.2005 um 16:30 schrieb Luigi Romolo:
> Hi everybody;
>
> I am about to create my own user defined GrADS basemaps.
> Essentially, what
> I am trying to do is create state maps that I can bring into the
> GrADS map
> library. I would also like to create my own poly.asc files for use
> with
> BASEMAP.gs so that I can mask out everything except my state
> basemap. Does
> anyone have a program that would do this ???? Any progress toward a
> solution to this dillema is greatly appreciated
>
> Cheers !
>
> Luigi Romolo
>
Dipl. Met. Frank Bandle
frank at bandle.com * FeMax Wetteragentur IT & Medienberatung
More information about the gradsusr
mailing list