c******************************************************** c Name: genhrap.f c c Purpose: Write the HRAP coordinates for a selected region of c cells to a file and create a subsequent file in geographic c coordinates in a suitable format to serve as input to the c GENERATE (polygon) command in ARC/INFO. This program is designed c to be followed by genhrap.aml. c c Two options are available for defining the region of cells to be c created --(1) Define the latitude and longitude extent of the region c to be mapped, or (2) Specify the SW corner of the grid to be created c and the number of colums and rows of cells to be created. With c either option, the program computes the HRAP coordinates of the c SW corner (if necessary) and generates grid cells starting with c the bottom row, moving from left to right, and then c moving to the next row up and repeating. c c Comments: Only the output files hrap.COD.dat and inputgc.COD are c required as input to genhrap.aml. Intermediate files and optional c files that were created in an earlier version are also listed below. c Code can be compiled on a UNIX F77 compiler. I do not know whether c the syntax is compatible with other compilers. c c Calls subroutines: wll, topoly, crdat(numx,numy,xstart,ystart) c c Inputs: none c Output: "COD" is a user defined suffix c hrap.COD = file of hrap coordinates /*temporary c geoc.COD = file of geocentric coordinates /*temporary c hrap.COD.dat = file containing HRAP coordinates in a format that c can be attached to the polygon attribute table c *pster.COD = file of polar stereographic coordinates c inputgc.COD = input file of geoc. coordinates to make a polgon c coverage c *inpster.COD = input file of polar stereographic coordinates to c make a polygon coverage c *inhrap.COD = input file of HRAP coordinates to make a polygon c coverage c c A * denotes optional files -- lines corresponding to the creation c of optional files have been commented out c in this version of the code. c********************************************************************** program genhrap c <<< Variable Declaration >>> parameter (maxcol = 336, maxrow = 160) c *** maxcol and maxrow are limited to the extent of HRAP cells for c which data is available in the Arkans.-Red River Basin integer xstart,ystart,numx,numy,numpts,numx1,numy1 double precision xhrap(maxcol), yhrap(maxrow) integer count,bool,rfunit,wfunit c *** rfunit and wfunit store the readfile unit number and the c *** writefile unit number to be passed to the subroutine topoly. character suff*3,file1*8,file2*8,file3*12,file5*11 c *** c <<< End of variable declaration >>> c *** Allow two options for defining the study region. write(*,*) 'Enter 1 to specify your region by latitudes' write(*,*) 'and longitudes of its corners.' write(*,*) '' write(*,*) '-OR-' write(*,*) '' write(*,*) 'Enter 2 to specify your region by hrap coordinates' write(*,*) 'and number of columns and rows.' read(*,*) bool if (bool.eq.1) then call llinput(xstart,ystart,numx1,numy1) else write(*,*) 'Enter the hrap(x,y) for the lower left hand' write(*,*) 'corner of the region of interest:' read(*,*) xstart,ystart write(*,*) 'Enter the number of grid columns and rows to be' write(*,*) 'created:' read(*,*) numx1,numy1 endif c *** Number of points to write is one greater than the number of c *** columns or rows. The name numx1 can be thought of as number of c *** x coordinates - 1. numx = numx1 + 1 numy = numy1 + 1 write(*,*) 'Enter a 3 character suffix to uniquely identify' write(*,*) 'your grid:' read(*,*) suff c ***Create names for all of the output files. c *** file1 = file of hrap coordinates c *** file2 = file of geocentric coordinates c *** file3 = file containing HRAP coordinates in a format that can be c attached to the polygon attribute table c *** file4 = file of polar stereographic coordinates c *** file5 = input file of geoc. coordinates to make a polgon c coverage c *** file6 = input file of p. stereographic coordinates to make a c polgon coverage c *** file7 = input file of hrap coordinates to make a polgon coverage file1 = 'hrap.'//suff file2 = 'geoc.'//suff file3 = 'hrap.'//suff//'.dat' c file4 = 'pster.'//suff file5 = 'inputgc.'//suff c file6 = 'inpster.'//suff c file7 = 'inhrap.'//suff open(unit = 10, file = file1, status = 'unknown') open (unit = 20, file = file2, status = 'unknown') open (unit = 30, file = file3, status = 'unknown') c open (unit = 40, file = file4, status = 'unknown') open (unit = 50, file = file5, status = 'unknown') c open (unit = 60, file = file6, status = 'unknown') c open (unit = 70, file = file7, status = 'unknown') c *** Compute the total number of cell corners numpts = numx*numy xnew = xstart do 100 i=1,numx xhrap(i) = xnew xtemp = xnew + 1.0 xnew = xtemp 100 continue ynew = ystart do 200 j=1,numy yhrap(j) = ynew ytemp = ynew + 1.0 ynew = ytemp 200 continue count = 1 do 300 j=1,numy do 400 i=1,numx c ** with "count" in the file, this file can be used as input c ** for creating an ARC/INFO point coverage write(10,*) count,xhrap(i),yhrap(j) count = count + 1 400 continue 300 continue c write(*,*) ystart call wll(numpts) rfunit = 20 wfunit = 50 call topoly(numx,numy,rfunit,wfunit) c rfunit = 40 c wfunit = 60 c call topoly(numx,numy,rfunit,wfunit) c rfunit = 10 c wfunit = 70 c call topoly(numx,numy,rfunit,wfunit) call crdat(numx,numy,xstart,ystart) close(10) close(20) close(30) close(50) stop end c<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>> c ********************************************************************* c c Purpose: Convert HRAP coordinates contained in file "hrap.COD" to c lat-long coordinates based on a spherical earth and write them to an c output file that can be used to generate a point coverage. The c parameter "numpts" stores the number of entries that will be expected c from "hrap.COD" c c Inputs: file hrap.COD c Outputs: file geoc.COD c c********************************************************************* subroutine wll(numpts) double precision xhrap, yhrap, x, y double precision bigr, arg, latd, lond, ang double precision stlatd,earthr,mesh,stlond integer rec,numpts c c*** Define constants stlond = -105.0 stlatd = 60.0 c*** earthr,mesh, x, and y are in meters. earthr = 6371200.0 mesh = 4762.5 rewind(unit=10) do 100 i=1,numpts read(10,*) rec,xhrap, yhrap x = (xhrap - 401.0)*mesh y = (yhrap - 1601.0)*mesh bigr = (x*x + y*y)**0.5 arg = bigr/(earthr*(1 + dsind(stlatd))) latd = 90.0 - 2*datand(arg) ang = datan2d(y,x) if (y.gt.0) then ang = 270.0-stlond-ang else ang = -90.0-stlond-ang endif if (ang.lt.180) then lond = -1 * ang else lond = 360.0 - ang endif c*** Write polar stereographic coordinates and geocentric c*** coordinates to a file. c write(40,*) i,x,y write(20,*) i,lond, latd 100 continue return end c123456789012345678901234567890123456789012345678901234567890123456789012 c********************************************************************** c Purpose: Given a list of corner points for a grid (can be (ID,x,y) or c (ID, lon,lat) in which the coordinates for the bottom row are c listed one per line followed by the coordinates for the next c row up, create a file that can be used to generate a polygon c coverage of the grid cells. c c Input: File of corner points (ID,x,y), c Ouput: File with lines: "poly-id, ll,lr,ur,ul,ll,end" -- repeated for c each polygon. ll = lower left, lr = lower right, ur = upper c right,ul = upper left c c********************************************************************** subroutine topoly(numx,numy,rfunit,wfunit) c <<< Variable Declaration >>> c parameter (numx = 20, numy = 20) c*** The old number of x-coordinates was 336. c*** The old number of y-coordinates was 160. double precision xrowa(336),yrowa(336),xrowb(336),yrowb(336) c ** xrowa, yrowa are x and y coordinates of points in row a character*3 end integer i,l,rcount,r,polynum,numx,numy integer rfunit,wfunit c <<< End of Variable Declaration >>> end = 'end' rewind(unit=rfunit) rcount = 1 polynum = 1 do 200 i=1,numx read(rfunit,*) rec,xrowa(i),yrowa(i) 200 continue 100 if (rcount.lt.numy) then do 250 i=1,numx read(rfunit,*) rec,xrowb(i),yrowb(i) 250 continue l = 1 300 if (l.lt.numx) then r = l + 1 write(wfunit,*) polynum, xrowa(l), yrowa(l) write(wfunit,*) xrowa(l),yrowa(l) write(wfunit,*) xrowa(r),yrowa(r) write(wfunit,*) xrowb(r),yrowb(r) write(wfunit,*) xrowb(l),yrowb(l) write(wfunit,*) xrowa(l),yrowa(l) write(wfunit,*) end l = l + 1 polynum = polynum + 1 goto 300 endif rcount = rcount + 1 do 350 i=1,numx xrowa(i) = xrowb(i) yrowa(i) = yrowb(i) 350 continue goto 100 endif write(wfunit,*) end return end c ********************************************************************** c Purpose: This subprogram will create a data file that can be joined to c the projected "hrap" polygon coverage so that "hrap" coordinates of the c lower left hand corner of each polygon will be added to the appropriate c line in the PAT. c c Note: The only difference between "hrap.COD.dat" produced by this c subroutine and "hrap.COD" produced by the main program is that c hrap.COD.dat does not contain entries for the last column and last c row of points. c ********************************************************************** subroutine crdat(numx,numy,xstart,ystart) c*** Old value of numx was 336 c*** Old value of numy was 160 double precision xhrap(336), yhrap(160) integer count,numx,numy,xstart,ystart,numx1,numy1 numx1 = numx - 1 numy1 = numy - 1 xnew = xstart do 100 i=1,numx1 xhrap(i) = xnew xtemp = xnew + 1.0 xnew = xtemp 100 continue ynew = ystart do 200 j=1,numy1 yhrap(j) = ynew ytemp = ynew + 1.0 ynew = ytemp 200 continue count = 1 do 300 j=1,numy1 do 400 i=1,numx1 write(30,*) count,xhrap(i),yhrap(j) count = count + 1 400 continue 300 continue return end c********************************************************************** c At user's request, allow the user to input the latitude and c longitude of the four corners that are of interest in the c study. c c Note: The user should input geodetic coordinates. These c geodetic coordinates will be interpreted as geocentric coordinates c to be consistent with methodology used by the c National Weather Service. c********************************************************************** subroutine llinput(xstart,ystart,numx1,numy1) c <<< Variable Declaration >>> parameter (stlat = 60.0) c*** clon is a constant used to account for the standard longitude c*** see eqn. in "Geographic Positioning of the HRAP" parameter (clon = 15.0) parameter (rad = 6371.2) integer xstart,ystart,numx1,numy1 real lon(4), lat(4) real sfactor,R,x,y,hrapx(4),hrapy(4) c*** Declare variables llhrapx and llhrapy to pick the hrap coordinates of c*** the lower left hand coordinates desired. real minhx,minhy,maxhx,maxhy c <<< End Variable Declaration >>> write(*,*) 'Enter the latitudes and longitudes of four' write(*,*) 'corners of a rectangle that encloses the' write(*,*) 'study region (in decimal degrees).' write(*,*) '' write(*,*) 'Enter a longitude value and then a space' write(*,*) 'and then a latitude value. Hit return after' write(*,*) 'each coordinate. Remember to input West' write(*,*) 'longitude values as negative numbers.' do 100 i = 1,4 read(*,*) lon(i),lat(i) sfactor = (1+sind(stlat))/(1+sind(lat(i))) c** x and y are in km R = rad*cosd(lat(i))*sfactor x = R*cosd(lon(i)+clon) y = R*sind(lon(i)+clon) hrapx(i) = x/4.7625 + 401 hrapy(i) = y/4.7625 + 1601 write(*,*) 'hrapx, hrapy:', hrapx(i), hrapy(i) 100 continue minhx = hrapx(1) minhy = hrapy(1) maxhx = hrapx(1) maxhy = hrapy(1) do 200 j = 2,4 if (hrapx(j).lt.minhx) then minhx = hrapx(j) endif if (hrapy(j).lt.minhy) then minhy = hrapy(j) endif if (hrapx(j).gt.maxhx) then maxhx = hrapx(j) endif if (hrapy(j).gt.maxhy) then maxhy = hrapy(j) endif 200 continue xstart = minhx ystart = minhy numx1 = maxhx - minhx numy1 = maxhy - minhy write(*,*) 'Lower left, num rows, num columns' write(*,*) xstart,ystart,numx1,numy1 return end c*********************************************************************