ΛgPΛRATPg ftWP MfiTHOP FOR C9MgTRPCTIW9

A MOSAIC OP DATA

A portion of the disclosure of this patent document contains material which is subject to copyright protection. The copyright owner has no objection to the facsimile reproduction by anyone of the patent document or the patent disclosure, as it appears in the U. S. Patent & Trademark Office patent files or records, but otherwise reserves all copyright rights whatsoever.

FIELD OF THE INVENTION The present invention relates to display of radar images, and more particularly to producing such a display by combining a plurality of radar images into a mosaic of those images. The present invention also relates to electrical computers and data processing systems, and more particularly to applications of same to earth science such as weather. The present invention further relates to information processing system organization, and more particularly to applications using artificial intelligence with details of the artificial intelligence system such as for earth sciences such as weather. The present invention still further relates to communications, directive radio wave systems and devices (e.g. radar or radio navigation) , and more particularly to weather radar, plural radar, display circuits for image production or scan conversion for same, and directive radio wave systems and devices for position indicating (e.g. triangulation) for storm location.

BACKGROUND OF THE INVENTION There exist applications in which there is a need to combine several radar images into a single mosaic image covering a geographical area. The present invention fulfills that need.

There is a large number of weather radar sites distributed across the U.S.A. There is a national weather radar network consisting of the new NEXRAD SR-88D weather surveillance Doppler radars and the preexisting WSR-57 and WSR-74 non-Doppler weather radars. Each such site can essentially operate independently, but the ranges of these sites can overlap. It is necessary to assemble the data from these sites into a single mosaic showing weather conditions in a particular region or across the country. The present invention fulfills this need.

A mosaic of radar data can provide location and intensity of precipitation and severe weather throughout a large geographic area and is, therefore, useful to meteorologists and other persons responsible for monitoring and forecasting weather conditions. Creation of a mosaic of radar data usually requires conversion of data points from a local geographic coordinate system (wherein the positions there are defined relative to an individual radar site) to the geographic coordinate system of the mosaic. Conversion of data points from one coordinate system to another coordinate system (reprojection) is a computationally intensive task, but is necessary to ensure accuracy. Radar data from even one site, however, can consist of over 100,000 individual data values, each of which must be reprojected. A mosaic could be made up of data from dozens of such radar sites. There is therefore a need to accurately combine all of that data from all desired sites fast enough to provide a current depiction of the weather. The present invention fulfills this need.

SUMMARY OF THE INVENTION

Accordingly, an object of the present invention is to provide apparatus and method for combining data from a plurality of sites into a single image covering a desired geographical area.

Another object of the present invention is to provide apparatus and method capable of combining, in real time,

radar data (such as weather radar data) in digital format from multiple radar sites into a mosaic covering a geographical area.

A further object of the present invention is to provide apparatus and method for combining data from a plurality of sources into a single image that can be run on a general purpose computer such as a personal computer.

Still another object of the present invention is to provide apparatus and method capable of producing a mosaic of radar data that can provide location and intensity of precipitation and other meteorological phenomena throughout a large geographic area.

A still further object of the present invention is to provide apparatus and method capable of accurately combining a substantial amount of data from a plurality of desired sites fast enough to provide a current depiction of the weather or other phenomena.

Yet another object of the present invention is to provide apparatus and method capable of accurately yet efficiently building mosaics of radar data.

Briefly, these and other objects of the present invention are accomplished by apparatus and method for converting each of a plurality of radar data bins to a corresponding grid location or box in a mosaic image. However, while each radar data bin is mapped to a grid box for the mosaic, not every such grid box within the coverage area of the radar may receive a bin. This may occur because there may be more grid boxes available than there are radar bins, or because the radar data bins do not overlap perfectly. After the mapping from bins to boxes is completed, the inventive apparatus and method then goes back and fills those holes in the following manner. Bins of radar data are mapped to a predetermined grid for a portion of the resulting display. Each time that a grid element of that grid receives a bin, then a corresponding element in the temporary array is "set". The temporary array is then checked for any "unset" or blank elements.

If an array element is found without a bin, then the locationally closest data bin is found. That closest data bin and its corresponding grid location are used to fill the hole in the mosaic grid.

Other objects, advantages and novel features of the invention will become apparent from the following detailed description of the invention when considered in conjunction with the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

In the drawings,

Fig. 1 is a diagrammatic representation of a two dimensional array for the display of a single site radar product;

Fig. 2 is a diagrammatic representation of one example of mapping of individual weather radar product data into a larger mosaic grid;

Fig. 3 is a diagrammatic representation of an overview of a process according to the present invention to indirectly map individual radar data values to a mosaic grid according to the present invention;

Fig. 4 is a diagrammatic representation of an overview of a process to directly map individual radar data values to a mosaic grid according to a preferred embodiment of the present invention;

Fig. 5 illustrates the format of a lookup table that can be utilized in the method of Fig. 4;

Fig. 6 shows one illustrative example of a display resulting from radar data mapping that has holes present where the source data grids and the destination data grid do not match up;

Figs. 7A and 7B together provide a diagrammatic representation of one embodiment of a method of filling any one of the holes of Fig. 6 according to the present invention;

Figs. 8A, 8B, 8C, 8D and 8E together show a flowchart for one method of building the lookup table of Figs. 4 and 5; and

Figs. 9A, 9B, 9C, 9D and 9E together show a flowchart for one method of building a mosaic of radar data according to the present invention using the lookup tables of Figs. 4, 5 and 8A through 8E.

DETAILED DESCRIPTION

Referring now to the drawings, wherein like reference characters designate like or corresponding parts throughout the several views, there is shown in Fig. 1 an array 10 made up of a plurality of elements or pixels, of which one is identified by reference numeral 12. Each array element (such as element 12) contains a single value of radar return. Each such element represents a finite area on a horizontal plane essentially tangent to the earth's surface. Array 10 contains one example of a decoded single site radar data product 14. In array 10, product representation 14 is bounded by circle 16 representing the extent of weather radar data return within the product. Within product 14, an individual radar data bin location (such as location 18) is identified by its x, y coordinates in the array, with the bin latitude and longitude being computed based on its distance relative to the radar site. Central bin 20 represents the radar site location. For a fixed radar site, the latitude and longitude of the radar site is known, so that referencing an individual bin with respect to location of the central bin will provide precise position information for that individual bin also. The single site radar data is received in digital format which allows the data values to be assigned to a computer program's memory in a two-dimensional array 10. Each array element (such as element 12) is a data "bin" which represents an atmospheric parameter (e.g. precipitation intensity) within a relatively small finite area such as 1 kilometer by 1 kilometer. Thus, the coordinates of any one

bin (such as bin 18) can provide the geographic x, distance of that bin from the radar site.

However, often more than one radar product must be used in a single grid array or visual representation of weather conditions in a region. This may be because the range for two radar sites may overlap, or because it is desired to show, in one display, radar data for a geographic area that is sufficiently large to include a plurality of radar sites. In such situations, a mosaic is built by taking the individual data bins from the radar data for a single site and assigning the values for those bins to the corresponding locations in a large scale (e.g. regional or national) grid which preferably also resides in a computer memory as a two-dimensional array 22. A simple example of formation of such a mosaic is shown in Fig. 2. Fig. 2 shows one example of mapping of individual weather radar product data from arrays 10, 24 and 26, such as from different radar sites, into a larger scale mosaic grid 22. Each array 10, 24 and 26 here represents radar product from a single site A, B, or C respectively. Each such array 10, 24 and 26 is similar to array 10 of Fig. 1. As shown in Fig. 2, the data from array 10 is combined with data from other such arrays 24 and 26 to construct a mosaic grid 22 which covers the geographic area that includes all three such radar sites. Thus, in the example of Fig. 2, array 10 is mapped onto portion or subarray 28 of grid 22; array 24 is mapped onto portion or subarray 30 of grid 22; and array 26 is mapped onto portion or subarray 32 of grid 22. In the example of Fig. 2, the coordinate system of the radar data from a single radar radar does not necessarily match the coordinate system of the mosaic grid 22. That is, the bins of arrays 10, 24 and 26 and the grid boxes of grid 22 are not necessarily the same size, and the respective coordinate systems are not necessarily oriented the same way. For this reason, the coordinates of each bin from an array 10, 24 or 26 must be converted (i.e. reprojected) to the coordinates of grid 22.

One method of reprojecting the coordinates of an individual data point such as for transfer from array 10 to grid 22 according to the present invention is shown in Fig. 3. Fig. 3 shows an overview of a process according to the present invention to indirectly map individual radar data values such as from array 10 to a mosaic grid such as grid 22. In single site radar product array 10, data values are assigned to bins. In the method of Fig. 3, each bin is converted to latitude and longitude and then that latitude and longitude is converted to the corresponding grid location. At step 34, each set of bin coordinates x,y is converted to distance from radar site 20 on the x and y axes. At step 36, this distance in the form of x and y values is then converted to geographic latitude and longitude. At step 38, these values of latitude and longitude are then converted to x, y distance from the grid 22 center, now using the x and y axes of grid 22. At step 40, this x, y distance is converted to the corresponding x and y coordinate values for grid 22. In this manner, data values from an array 10 are mapped to mosaic grid 22. A similar procedure can be utilized for mapping from other arrays 24, 26 to grid 22. Although this approach is accurate, it is not preferred because it is computationally intensive (see steps 34, 36, 38 and 40) and therefore would not be desirable for real-time operation on a general purpose computer such as a personal computer.

An alternative, more preferred method for mapping from a single site grid or product array 10 to a mosaic grid 22 is illustrated in Fig. 4. The method of Fig. 4 is more efficient in real time than that of Fig. 3 because the method of Fig. 4 maps the radar data bins directly to the grid 22 using prebuilt look-up tables. These look-up tables are built using computationally intensive algorithms, illustrated in Fig. 8 and shown in the corresponding code for that figure, to compute latitude and longitude of each radar data bin and to convert those latitude and longitude values to box locations of grid 22.

Weather data radar sites usually have fixed locations, that the location and orientation of the product array 1 24, or 26 for each such site with respect to that of gr 22 is known. These look-up tables are preferably bui offline and are preferably stored as files in a comput database 46. The program of Fig. 9 which builds mosai then reads in the appropriate look-up table file as i processes data from each site. Each look-up table contain a list of bin x,y values and the corresponding grid x, values preceded by the number of entries in that table. preferred format for such a look-up table is shown in Fig 5. Each entry in the look-up table is then processed b assigning the radar data value from bin x,y to th corresponding grid x,y location.

Fig. 4 shows an overview of a preferred process t directly map individual radar data values to a mosaic gri using a prebuilt look-up table. In Fig. 4, in single sit radar product array 10, data values are assigned to bins i that array. In the method of Fig. 4, a direct mapping of for example, bin 42 to box 44 is accomplished in th following manner. A look-up table for each site is store in database 46. When a mapping of, say, array 10 to gri 22 is to take place, then the look-up table for the sit for array 10 is read in from database 46 at step 48. Then each bin (e.g. bin 42) of array 10 is mapped to a box (e.g. box 44) of grid 22 by converting that bin's x and coordinate values directly to x, y coordinate values fo grid 22 by finding the former in the look-up table an reading out the corresponding grid box location from tha table, at step 50. Grid 22 is the mosaic grid to whic data values such as from array 10 are mapped. In effect, with the method of Fig. 4, steps 34, 36, 38 and 40 of Fig. 3 are performed in advance offline to reduce the amount of processing required during real-time operation.

Fig. 5 shows one example of a file format for a look¬ up table of Fig. 4. In Fig. 5, all entry values are integers. Each file is site specific. Also, each file is

specific to the coordinate system and spatial resolution of the individual product (e.g. product 14) and of the mosaic grid (e.g. grid 22). A given bin x, y coordinate can be paired with more than one grid x, y coordinate. Likewise, a given grid x, y coordinate can be paired with more than one bin x, y coordinate. However, there is at least one table entry for each bin x, y coordinate, and there is at least one entry for each grid x, y coordinate. However, the total number of entries likely would not exceed the larger of the number of bins or the number of boxes to be mapped.

As discussed above, the array of product data bins (e.g. array 10, 24 or 26) will not necessarily be perfectly aligned with mosaic grid 22. As a result, some grid boxes will receive more than one bin. For similar reasons, some grid boxes, within the radar data range, will not receive any bins. In the latter case, "holes" would remain in the completed mosaic, such as can be seen in Fig. 6. Fig. 6 shows one example of radar data mapping having one or more holes 52 where the source and destination grids, i.e. the radar bins and mosaic grid, do not match up. The following additional processing is utilized in the present invention to fill such holes 52 in grid 22. This process is accomplished in the following manner. In the process to build look-up tables as shown in Fig. 8, the location of each such "hole" in the grid is identified and the closest bin is assigned to the grid box having that hole. This is done by using the temporary array of Fig. 7 to find unfilled grid boxes within the radar coverage area. Each time a bin x, y coordinate is paired with a grid x, y position, an element in temporary array 54 is set. Then, when processing of all bin x, y locations for a site is complete, a "backward calculation" is performed for all unset or unfilled grid boxes in temporary array 54 within the range of the radar for which a look-up table is then being constructed. In such a backward calculation, the grid box x, y location is converted to latitude and

longitude, and then that latitude and longitude i converted to the corresponding bin x, y coordinate to find the bin closest in location to that bin x, y coordinate for the hole. Each so "back calculated" bin x, y and grid x, y pair is added as an entry to the look-up table for that radar site. This process ensures that every data bin is assigned to the mosaic grid and that every grid box within the range of the radar then being considered is assigned a data value.

Fig. 7 shows use of a temporary array 54 to fill a hole in grid 22. Each element of temporary array 54 can contain any one of three values: a first, initial "blank" value (such as -1) ; a second "unset" value (such as 0) for each corresponding grid element or box falling within the range of a radar; and a third "set" value (such as +i) indicating that a radar bin has been mapped to the corresponding grid element. In a subarray (such as subarray 32) of grid 22, bins of weather radar data from the corresponding product array (e.g. array 26) are mapped to grid 22 with the orientation of grid 22 and not that of the source array then being mapped. In this mapping process, each time that a grid 22 element receives a bin, then a corresponding element in temporary array 54 is "set" to so indicate. Temporary array 54 is then checked for any "unset" or unfilled elements. In temporary array 54, circle 56 represents the range of weather data in the weather radar product then being mapped. All array 54 elements within circle 56 are to receive a data bin or value from the corresponding product array being mapped to the grid 22. If an array element of array 54 is found that has not been so "set", then a back calculation is performed to find the filled data bin that is closest in location to that array element. That data bin x, y and its corresponding grid location x, y are added to the look-up table to fill the hole that would otherwise result in mosaic grid 22.

To perform the process of Fig. 4, a two-part process to build mosaics is needed. The first of these two parts, shown in Figs. 8A-8E, can be executed off-line and is responsible for the generation of look-up tables which are built and stored in a local database 46. The second part, which is designed to run in real time and uses the pre- built look-up tables to actually build the mosaic(s) , is shown in Figs. 9A-9E.

In Figs. 8A-8E, the process to build site-specific look-up tables can be summarized as follows. The following process is conducted for each radar site, with a separate table being constructed for each radar site. Each element in temporary array ^" 54 is first initialized to reflect a "blank" or cleared value. For each element in the temporary array, the element x, y coordinates are converted to distance or range from the radar site; if this range is within the radar data coverage area defined by circle 56, then a marker accompanying that temporary array element in the database is given an "unset" indication. For each radar data bin location of the radar product then being mapped, the bin x, y coordinate is converted to latitude and longitude; the resulting bin latitude and longitude is converted to a grid x, y location; the resulting bin x, y and grid x, y pair so found is assigned to an appropriate location in the look-up table for that site; the temporary array element x, y location is computed for that pair; and the temporary array element at that x, y location so found is given a "set" indication, showing that no hole is present at that location. Thereafter, the following steps are performed for each element in the temporary array 54. It is first determined whether that temporary array element is unset and thus is in the radar data coverage area but without a radar data bin assigned to it. If so, then the temporary array element x, y coordinates are converted to grid x, y coordinates; the resulting grid x, y coordinates are converted to latitude and longitude; the resulting latitude and longitude values are converted to radar bin x.

y coordinates; and the resulting bin x, y coordinate and grid x, y coordinate pair are assigned to the look-up table for that single site radar product, to fill or to partially fill a hole in the mosaic grid. Finally, the resulting look-up table is assigned to a file in the database 46.

The process of Fig. 9 to build a radar data mosaic using the look-up tables can be summarized as follows. First, for each location or element in the large grid 22, that grid element is initialized to a suitable predetermined value. Then, the following steps are performed for each radar site. The latest weather radar data product for that site is read in from the database. That radar product is then decoded, and individual data values from that product are assigned to a two-dimensional array. The corresponding prebuilt look-up table is then read in from the database 46. Then, for each entry in that look-up table, extract the corresponding bin x, y coordinate and its corresponding grid x, y coordinate from the look-up table. The product data value at that bin x, y location is then extracted from the decoded radar product. That product data value is then assigned to the corresponding grid x, y location obtained from the look-up table. These last three steps are repeated for each entry in that look-up table. Finally, the resulting combined radar data in the grid is encoded and saved to file in the database 46.

Figs. 8A, 8B, 8C, 8D and 8E together show a structure diagram or flowchart for a method of building look-up tables. These look-up tables are essential components for the method of Fig. 4 and of its corresponding system for building mosaics of radar data in real time such as on a general purpose digital computer. The method of Figs. 8A- 8E can be implemented as a computer software program that preferably can be run offline when time is not critical, so that time-consuming computations are performed in the look¬ up table building process and thus well before the results of those computations would be needed.

15504

- 13 -

As s.- ^• wn in Fig. 8A, the program uses a manually created configuration file 101 that contains information on one or more sites for which a look-up table is to be built. Configuration file 101 contains the number of sites, and the latitude/longitude and a site I.D. number of each site. At step 102, the program reads the number of sites from file 101, and at step 103 reads in the site information from file 101 and stores that in the computer's memory for later use. Then, at step 104, the program enters a loop to process each such site on that list. This loop is executed until all sites are processed (determined at step 105) , at which time the program terminates (step 106) . The processing by that loop for each site begins at step 107 with the initialization of the temporary array by assigning the value of -1 to all elements of that array.

Next, the temporary array is set up by the program, as shown by Fig. 8B, to indicate which elements of that temporary array are located within the radar coverage area. At step 110, this program loops for each x coordinate in the temporary array until all x coordinates are processed (determined at step ill) . in this second loop, the temporary array x coordinate then being processed is converted at step 112 to the corresponding grid array x coordinate. At step 113, it is then determined whether the resulting grid x coordinate is outside of the boundary of the grid; if so, then the program skips to the next temporary array x coordinate (steps 110, ill and 112). This could occur for sites located near the east or west edges of the grid (for example of the region of interest) where the temporary array 54 extends beyond the edge of the grid 22. Otherwise, for a valid grid x coordinate, the program proceeds to another loop starting at step 114 for each y coordinate in the temporary array until all y coordinates are processed (step 115) . At step 116, the temporary array y coordinate then being considered by this third loop is converted to the corresponding grid array y coordinate. If at step 117 the resulting grid y coordinate

is found to be outside of the boundary of the grid, then the program skips to the next temporary array y coordinate (steps 114, 115 and 116). That could occur for sites located near the north or south edges of the grid (or of the geographic region of interest) where the temporary array extends beyond the edge of the grid. Otherwise, for a valid grid y coordinate, at step 118 the grid x, y position found by steps 112 and 116 is converted to latitude and longitude by using the projection equations that define the particular grid being used. For example, Lambert conformal projection equations could be used for this purpose, since a grid is treated as being flat, but is used to represent a curved surface, some projection will be used to transform the grid into the curved surface or vice versa; various ways to accomplish such projection are well known to cartography. At step 119, the latitude and longitude is then converted to x and y coordinate values for the radar coordinate system, i.e. the coordinate system for the single site radar product. For step 119, an appropriate cartographic (or the like) conversion from latitude and longitude to the coordinate system of the radar product array is needed. The opposite such conversion was accomplished at step 118, except that as discussed above, the grid coordinates would not necessarily correspond to those of the single site radar product array. The x coordinate and y coordinate found at step 119 is then used at step 120 to compute the distance (in the radar product array) of the location defined by the x, y coordinate from the radar site. At step 121, it is determined whether the distance found at step 120 is within the radar coverage area such as is defined by circle 16 of Fig. l; if so, then the value of the temporary array element, defined by the x, y coordinate utilized at steps 112 and 116, is set to zero at step 122, to indicate that the corresponding radar bin value is defined by the single site radar product array. Either way, the program goes on to the next y coordinate at step 114.

Once the temporary array is set up by the method of Fig. 8B to identify the elements that fall within the radar coverage area, generation of the look-up table proceeds, as shown by Fig. 8C. As shown in Fig. 8C, generation of the look-up table begins at step 130, where the number of look¬ up table entries is set to zero. Starting at step 131, the program then performs a loop for each y coordinate in the radar product raster array until all y coordinates in that array are processed (step 132). For each iteration of the loop beginning at step 131, the program performs a loop starting at step 133 for each x coordinate in the radar product raster array; this loop is exited at step 134 to return to step 131 when all x coordinates in the radar product raster array are processed by the loop beginning at step 133. In the loop beginning at step 133, using the current product raster array x, y coordinate, the distance of that coordinate (in that array) from the radar site is determined at step 135. At step 136, it is determined whether the distance computed at step 135 is outside of the radar coverage area; if so, then the program skips to the next x, y position by returning to step 133. Otherwise, the position is within the radar coverage area and further processing is performed, as shown in Fig. 8C. In this further processing, at step 137 the x, y coordinate then being processed is first converted to latitude/longitude. As used herein, coordinate or position refers to an x-value and a y-value that together (such as Cartesian coordinates) define a location in a system, array, grid, surface, or the like. At step 138, that resulting latitude/longitude is then converted to the corresponding grid coordinate system x, y position. At step 139, the resulting grid x, y coordinate is checked for validity. If at step 139 it is found that the grid x, y position found at step 138 is invalid (not within the preset grid boundaries) , then those coordinates do not belong in the look-up table then being constructed, and the program skips, to the next element in the radar product array being processed, by returning to

step 133. Otherwise, at step 140 the radar product raster array x, y coordinate, and the corresponding grid x, y coordinate found at step 138 for that radar x, y coordinate are, assigned to the look-up table array then being formed in the computer's memory. At step 141, the total number of look-up table entries for that look-up table is then incremented by one; recall that before the program entered the step 131 loop, that number had been initialized to zero at step 130. After step 141, the grid x, y coordinate found at step 138 and entered in the look-up table at step 140 is used at step 142 to compute the temporary array x, y coordinate corresponding to that grid x, y coordinate. At step 143, the value of the temporary array element at that x, y position is set to the value of plus one, to indicate that the corresponding grid box received at least one data value from the radar product raster array. The program then returns to step 133. If at step 132 it is found that all of the y coordinates in the radar product array have been processed by the loop beginning at step 131, then the program proceeds to the process shown at Fig. 8D.

After all elements in the radar product raster array have been processed by the method of. Fig. 8C, then the program locates all grid boxes within the radar coverage area that were not matched up with a coordinate or position from the radar product raster array, as shown in Fig. 8D. This is accomplished by searching the temporary array for all unset elements given the value of zero at step 122 of Fig. 8B. In Fig. 8D, a loop beginning at step 150 is performed for each x coordinate in the temporary array until there are no more x coordinates to process (step 151) . Within that loop, another loop beginning at step 152 is performed for each y coordinate for the x coordinate then being considered by the step 150 loop until all such y coordinates are processed (step 153) . In this other loop, after it is determined at step 153 whether there still remains any y coordinates to so process, at step 154

the value at the current x, y location in a temporary array is evaluated. At step 154, it is determined whether the temporary array element at the x, y location then being addressed by the two Fig. 8D loops is equal to zero. If the value of that location is non-zero, then it is either outside of the radar coverage area or else it has been matched with a radar data bin location. Otherwise, when the value at that x, y location is still unset (in this example, assigned a value of zero) , then a hole exists in the mosaic. For each such hole, a match is needed to the radar data bin coordinate closest in position to that. If such a hole exists, then the grid x, y coordinate is determined at step 155 from the temporary array x, y coordinate considered at step 154. At step 156, the resulting grid x, y coordinate is converted to latitude/longitude using the equations (discussed above) that define the grid projection. That latitude/longitude found in step 156 is then converted at step 157 to the closest x, y coordinate in the radar product array coordinate system. At step 158, it is determined whether the radar product x, y coordinate found at step 157 is valid; if not, then the process returns to step 153. At step 159, it is then determined whether the grid x, y coordinate determined at step 155 is valid; if not, then the process returns to step 153. If the radar product x, y coordinate is valid and the grid x, y coordinate is valid (i.e., they fall within the region of interest), then at step 160 those coordinates are assigned to the lookup table array in the next available location in that array. Should an invalid such coordinate occur, the process then skips to the next temporary array element at step 153. After the radar product x, y coordinate and the corresponding grid x, y coordinate are added to the lookup table array, at step 161 the total number of entries in the lookup table is incremented by one. From step 161, the process returns to step 153. If no further y coordinate(s) is (are) present, then the process returns to step 150 so that the next x

coordinate can be addressed in the manner described above for Fig. 8D. If no more such x coordinates remain to be processed, then the program goes from step 151 to step 170 of Fig. 8E.

Once the lookup table is completed in the manner shown in Fig. 8D, the program saves the table so determined in a file such as on the computer's disk in the manner shown in Fig. 8E. First, a unique file name for the lookup table is created at step 170 associating the current site ID number with the file. At step 170, the file is then opened for output. At step 172, the program then builds a lookup table header in memory containing, among other information, the number of entries in the table, determined at steps 130, 141 and 161. The program then writes to the new lookup table file 173 at steps 174 and 175. First, the header built at step 172 is written at step 174 to the lookup table file 173, followed at step 175 by the complete lookup table. At step 176, the lookup table file 173 is then closed. The program then returns to step 104 of Fig. 8A to continue processing until lookup tables for all sites in the site list are complete. It is only necessary to rerun the program of Figs. 8A-8E if one or more sites are moved or added to the mosaic; real time running of the program of Figs. 8A-8E is not needed.

Construction of a look-up table according to the present invention can be implemented in software. An example of such software, written in the Microsoft version of the C language, is given in the following pages.

However, other programming languages can be used for this purpose. Figs. 8A-8E together generally show a flowchart for this software. ..--_.__. * __•-.-----•-• -.--*-

*• UNISYS WZS MOSAIC PROCESSING SYSTEM ** Copyright 1994 ϋniβyβ Corporation

**

* * PROGRAM NAME : CR_LUT. EXE

** ^—

** DESCRIPTION:Read in the configuration file and build a lut for each _itβ. ** Map* Conpoaite Reflectivity Product 36 and 38 to a 4x4km National grid.

**

•• DATA FILES NEEDED:

*• CONFIG.DAT ( Binary data file of aitea for which lut file* ar

** to be built. Created by the config.exe program.

** Located in the local directory.)

**

** OUTPUT FILES CREATED1

** CR4 ???.LUT (Maps NEXRAD Compoβite Ref productβ to 4x4km grid.

** ^* ~ There muat be one for each NEXRAD aite.) The ???

** repreββntβ the NEXRAD site id number.

_*

********* ***** * **** /

/include <math.h> /include <_tdio.h> /include <βtdlib.h> /include <βtring.h> /include <limits.h> /include <metnory.h> /include <doe.h> /include /include /include /include

/define MAXWRITE 65534

/define NODATA 0 /define MIN NUM 0 /define MIN~VALUE 0 /define NO_VALUE 0 /define MIN X 0 /define MIN f 0

/define GRID X SIZB_4KM 1280 /define GRID~y~SIZE_4KM 896

/define LUT_MAX 190000

/define GRID_X_HALP_4KM 640 /define GRID_Y_HALP_4KM 448

/define PR_BINS_CR4 232 /define PR ROWS CR4 232 /define TMP_X_SΪZB_CR4 250 /define TMP_Y_SIZB CR4 250 /define TMP_HALF_CR4 125

/define BIN_X_IDX 0 /define BIN_Y IDX 1 /define GRID X_IDX 2 /define GRID~Y_IDX 3 /define LUT_WΣDTH 4

/define M_PI 3.14159265358979323846 /define M_PI_2 1.57079632679489661923 /define M_PI_4 0.78539816339744830962 long lut_key » {123456}; typβdβf struct{ long keyword; long num_entries; short header size; short site i3; short versTon; short datβ_built; unsigned long time_built; short grid_reβ; short grid_width; short grid_height;

short bin_res; short bin^width; short bin~~height; short projection; double grid_center_phi; double grid~center lamda; double radar_phi; ~ double radar_la_αda; } lut_header_format; lut_header_format lut_header; typedef struct{ double dx; double dy; double phi org; double la_-3a_org; double phi value; double lam3a_value; short »; ^" ~ short y; short pr id; short grld_id; } coords; ~

/define MAX_SITES 200 long site_key - {24680}; short nexrad_βite ■ {1}; short rrwdβ_βite * {2}; short yes - {1}; short no ■ {0}; typedef struct{ long keyword; short num_siteβ; } βite_header_format; βita_header_foπnat βitβ_file_headβr;

βite_contentβ_format βite_file_contentβ[MAX_SITBS]; double max_rangβ » {464.0}; double phi center_deg * 38.0; double lam3a_center_dβg * -98.0;

/*

** Prototypes.

*/ short rnd of (double fl); int 1 c kό_21atlon( coords *l_coordβ ); int l~c~latlon2km( coords *l~coordβ ); int l ^" "c latlon2grid( coords ⁷ l_coordβ ); int l~c~grid21atlon( coords *l~coordβ ); int nβx~ka_21atlon( coords *l_cδordβ );

- 21 - int nex_latlon2km( coords *1 coords ); int nex Latlon2bin( coords *T_coordβ ); int nex~bin21atlon( coords *l~coordβ ); long hugβ_fwrite( void huge~*buffer, long length, βize_t count, int fEw) ;

/*

** Main program.

*/ int main()

PILE * fpr; int fhw;

/*

** Initialize the local_coordβ "pointer" so that

** it points somewhere.-

*/ local_coordβ ■ &dummy_coordβ ;

/*

** Open the configuration file for "read binary".

*/ fpr-f open ( "config.dat" , "rb" ) ; if ( fpr— NULL)

₍ voi _{d )} printf ₍ " _E rror opening the configuration fil ^β . \n" ) ; return(l) ; }

/*

** Read in the configuration file header.

*/

- 22 - nread«fread(fiaite file header,sizeof(site header format),1,fpr); if (nrβad«*NODATA)~ ~ ~

<

(void)printf ("Error reading the configuration file header. \n") ; return ( 1 ) ; }

/*

^* ** Check the configuration file header for the correct keyword

** and exit the program if the keyword does not agree.

*/ if ( site file header . keyword 1- site key ) {

(void)print ( "Invalid configuration file keyword\n" ) ; return(l) ; >

/*

** Extract fehe number of sites from the configuration file header */ num_βites * βitβ_file_hβadβr.num_βiteβ;

/*

** Print the maximum number of sites.

*/

(void)printf("\nbuilding mosaic look up tables for %d

/*

** Read in the contents of the configuration file. The amount to be read in

** is determined by the number of sites and the size of the structure for

** each site.

*/ nread*fread(βite file contents, num βiteβ*βizeof (βitβ_contentβ_f ormat) , l, fpr) ;

) /

/*

** Close the configuration file.

*/ status - fcloβe(fpr); if( status )

{

(void)printf("Error closing the configuration file.\n"); return(l);

}

/*

*• Loop through all of the sites and build a lut for each.

*/ for (j«MIN_NUM; j<num_βiteβ; j++) if ( (site file contents! j ) . site flag »- nexrad site)££

(βite~file~contβntβ(j J .βitβ ^" build lut — yes) ) {

/*

** Initialize the temporary array with the following.

** The array is initialized by setting every array element.

** This array is used to determine where holes exist in the

** initial mapping of NEXRAD data bins to the national mosaic. */ for (tx-MIN X; tx<TMP X SIZE CR4; tx++) { ^" for(ty«MIN Y; ty<TMP Y SIZE CR4; ty++)

{ tmp{txJ [t ]»tmp init;

) }

/*

** Extract the radar site latitude and longitude and assign the

** values to simple variables and assign to the local coords structure.

*/ phi radar_deg«βite file contents[jJ.site lat; lam3a_radar_deg-βite_fiTe_contentβ(jJ.βite_lon; locarl_coordβ->phi value - phi_radar_deg; local~coorda->lam3a_value ■ lamda_radar_deg;

/•

** Print the site id and radar latitude/longitude.

*/

(void)printf("file/ %3.3d NEXRAD site: %3.3d ",j,βite file contβntβfj) .βite_id); (void)printf(" radar lat/lon: %f/%f (deg)\n", phi_radar_deg,lamda_radar_deg);

/*

** Convert the radar latitude/longitude of the current site to an x,y

** coordinate in the national grid with the following routine

** (which uses the Lambert conformal projection). The x and y

** coordinates are stored in the simple variables: grid radar x, grid radar y.

*/ ^" statuaa-l c latlon2grid( local coords ); if(status) ^{" ""}

(void)printf("Error from l_c_latlon2grid: ,status); return(1) ; > grid_radar_x ■ local_coordβ->x grid_radar~y ■ loca_ coordβ->y

/*

** Next, unset the array elements which fall inside of the

** range of the radar. This is done because this is the region

** which must be checked for "holes" left by mapping NEXRAD bins to

** the grid. First, loop through the temporary array in the x direction.

*/ for (tx-MIN X; tx<TMP X SIZE CR4; tx++) {

/*

** Convert the temporary array x coordinate to the corresponding

** grid box x coordinate. */ grid_x-(grid_radar_x - TMP_HALF_CR4 + tx);

/*

•• Check that the x grid coordinate falls in the array.

** Otherwise, the position is just skipped.

*/

if((grid x>-MIN X)£_(grid x<GRID X SIZE 4KM) ) _{ ~ - - -

/*

** Loop through all of the y bins in the temporary array.

*/ for(ty«MIN Y; ty<TMP Y SIZE CR4; ty++) { ^"

/*

** Convert the temporary array y coordinate to the

** corresponding grid box y coordinate.

*/ grid_y«(grid_radar_y - TMP_HALF_CR4 + ty);

/*

** Check that the y grid coordinate falls in the array.

** Otherwise, the position is just skipped.

*/ if((grid y>«MIN Y)fi&(grid y<GRID Y SIZE 4KM)) { ~

/*

** Convert the x,y grid box coordinate to

** a latitude longitude value.

*/ local_coordβ->x ■ grid_x; loca_2coordβ->y ■ grid~y; status ■ l_c_grid21atlon( local_coordβ ); phi deg -> local_coordβ->phi_valuβ; lam3a_deg ■ local_coorda->lamda_value;

/*

** Then convert that latitude longitude position

** to an x,y position in tan's from the radar in

** the NEXRAD coordinate system.

*/ local_coordβ->phi value « phi_deg; local~coords->lam3a_value - lamda_deg; local~coordβ->phi_org ■ phi_radar_dβg; local2coordβ->lamda_org ■ lamda_radar_deg; status * nex_latlon2km( local_coordβ );

/*

** Then use the x,y coordinate in km'β to determine the

** actual distance from the radar in km'β. Unset the

** temporary array value if the distance is within the

** range of the product.

*/ nx ■ local_coordβ->dx ny ■ local~coordβ->dy range»βqrt((nx*nx)+(ny*ny)); if(range<max range)

{ tmp[tx] (ty]«tmp_unset;

}

/*

** The look up table is built here. The first step is to "map"

*• each NEXRAD bin to a grid x,y coordinate. The overall approach

** is to convert the bin location to latitude/longitude and then

** convert that latitude/longitude to a grid x,y coordinate

** (which is in the Lambert conformal projection) . The look up

** table consists of a list of entries, each containing a NEXRAD

** bin and row coordinate and the corresponding grid x and y coordinate.

** Later on, in the second step, holes in the grid are found and the

** closest NEXRAD bin is mapped to the grid x,y coordinate of the hole.

/•

** Initially, there are no entries in the look up table. */ entry»NO_VALUE;

/*

** Af-βign the radar latitude longitude value of the current site

** to the local coords structure.

*/ local_coordβ->phi org - phi_radar_dβg; local_coordβ->lam3a_org ■ lamda_radar_deg;

/*

** Loop for each NEXRAD bin in the product. The outer loop

** handles each row of data.

*/ for(bin_y-MIN_Y; bin_y<PR_ROWS CR4; bin y++) {

/*

** The inner loop handles each bin in a row.

*/ for(bin_x-MIN X; bin_x<PR BINS CR4; bin x++) {

/*

•• For the current bin, which is in the NEXRAD

** coordinate system, determine the latitude longitude.

*/ local_coordβ->x - bin_x; local_coordβ->y ■ bin~y; βtatuβ->nex_bin21atlon( local_coordβ ) ; phi deg » local_coordβ->phi_value la_o3a_deg _* local_coordβ->larada_value; if ( βtatuβ _* «good_βtatuβ ) {

/*

** For a valid x and y bin, take that latitude/longitude

** and compute the corresponding grid x,y coordinate.

*/ local_coordβ->phi value ■ phi_dβg; local~coordβ->lam3a_valuβ - lamda_dβg; βtatua-l_c_latlon2grid( local_coordβ ); grid_x ■ local_coordβ->x grid~y - local~coordβ->y

96/15504

- 26 - if(βtatus««good status)

{

/*

** For a valid latitude/longitude, the corresponding

** bin x,y and grid x,y coordinates are assigned to the

** look up table array and the number of entries

** is incremented.

V lut[entry][BIN X IDX] • bin x; lut[entry][BIN~Y~IDXJ » bin~y; lut[entry][GRID X IDX] » grϊd x; lut[entryJ[GRID tD J - grid ^' y; entry++;

I*

** Then the grid x,y coordinate is used to compute the

** temporary array x,y coordinate. If the temporary

** array x,y coordinate is valid for the array, then

** the array element is set. Any unset array elements

** represent "holes" which will be handled later.

*/ tx-grid x - grid radar x + TMP HALF CR4; ty«grid~y - grid~radar~y + TMP~HALF~CR4; if((tx»M_N X)&S(tx<TMP X SIZE CR4)S&

(ty>«MIN Y7sϊ(ty<ΪMP Y SIZE CR4)) . - tmp[tx][ty]-tmp set; >

}

/*

** In the next step, go through the temporary array and find any unset

** array elements which represent grid boxes within the product

** coverage area that did not already receive a bin. For each of those

** grid coordinates, find the closest NEXRAD data bin and "map" it to

** the grid. This is where the holes are filled. Additional entries

*• are made to the look up table for each set of NEXRAD bin,row and

** corresponding grid x,y coordinatea.

*/ for (tx-MIN X; tx<TMP_X_SIZE_CR4; tx++)

{ for(ty-MIN Y; ty<TMP_Y_SIZB CR4; ty++)

{

/*

** Check the current temporary array element to determine if

** there is a hole with the following. (Otherwise, the loop jus

** goes on to the next element.)

*/ if( tmp[tx][ty]--trap unset )

{

/ ^•

** Convert the temporary array x,y coordinate to a grid x,y

** coordinate.

*/ grid x-(grid radar x - TMP HALF CR4 + tx); grid ^r -(gri ^adar Z " TMP~HALF~CR4 ^♦ ty);

/* Convert the grid x,y coordinate to the corresponding *• latitude/longitude.

* / local_coordβ->x - grid_x; local~coordβ->y - grid~y; status > l_c_grid21atlon( local_coordβ ); phi deg « local_coords->phi_value lam3a_deg « local_coordβ->l__αda_value; if(status — good status)

{

/*

** For a valid grid x,y coordinate, convert the

** latitude/longitude to the corresponding (closest)

** NEXRAD data bin (store in: bin x, bin y) .

*/ ^"* local_coordβ->phi_value - phi_dβg; local2coords->lamda_value - lamda_deg; local~coordβ->phi_org - phi_radar~deg; local coordβ->lamda_org - lamda_radar_deg; status - nβx_latlon2bin( local_coordβ ); bin_x - local_coordβ->x bin_y - local~coordβ->y if(status « good status)

{ if((bin x>«MIN X)M(biπ x<PR BINS CR4)&&

Jbin y>-MIN Y)S_.7bin y<PR_RO S CR4)&6 (gri3 x>-MlR X)fi6(grld x<GRID X SIZE 4KM)&& (grid~y>-MIN~Y)-:fi(grid~y<GRID~Y~SIZB~41M))

{

/ ^•

** When the bin coordinates and the grid x,y ** coordinates are valid, then assign those •* coordinate pairs to the-look up table array ** and increment the number of entries.

lut(entry][BIN X IDX]-bin x; lut[entry] [BIN~Y~IDX]-bin~y; lut[entry][GRID X ^* IDX]-grid x; lut(entry] (GRIDJf~IDX)-grid~y; entry++; ~

/*

** Check the number of entries to determine if the

** unlikely condition has occurred that the look up table

** is full. In that case, print a message and break

** out of the loop.

*/ if(entry -- LUT MAX)

{

(void)printf("LUT is full: entry-%d\n",entry); goto lut full;

>

} lut_full:; /* The goto branches here to break from loop. Do nothing. */

/*

** Create the look up table (lut) file name in which the site id number

** is embedded with the following.

*/

(void)βprintf(lut_name,"cr4_%3.3d.lut\0",βite_file_cont βntβ[j].βite_id);

/*

** Print the lut file name.

*/

(void)printf(" look up table written to lut_namβ);

** Create and Open the lut file for "write".

*/

I

_fmode - 0_BINARY; fhw-creat(lut name,S IWRITE); if(fhw—-1)

{

(void)printf("Error creating lut file: %β\n",lut_name); return(1); } lut_header.keyword - lut_key; lut~header.num. entries - ^" entry; lut~header.hea3er size - βizeof(lut_header_format) ; lut~header.βite i3 - βite_file_contentβ[j]7βitβ_id; lut~header.version - 1; ~ ~ ~ lu- ^" header.date_built - 0; lut_header.time~built - 0; lut~header.grid~reβ - 4; lut~header.grid~width - GRID X SIZE 4KM; lut_headβr.grid~height - GRID___SIZB_4KM; lut~header.bin res - 4; ~ ~ lut_header.bin~width - PR BINS CR4; lut_hβader.bin~hβight - PR_ROWS_CR4; lut~header.projection - 1; lut~header.grid_centβr_phi - phi_center_dβg; lut_header.grid_cβnter_lamda - lamda_center_deg; lut~header.radar_phi -~phi_radar_deg7 ~ lut~header.radar_lamda - lamda_radar_deg;

/*

** Write the lut file header structure to the lut file

** with the following.

*/ doe status - _doβ_write(fhw,&lut_header, (size_t)βizeof(lut_headβr_format) , if(3os statue ^* ) ^" ~ ~

{

(void)printf("Error writing lut file header: %β\n",lut_namβ); return(l); > if(nwrite—NODATA)

{void)printf("Error writing lut file header, nwrite - %d\n",nwrite); return(1); >

** Compute the size of the lut file contents based on the number

** of entries and number of bytes per entry.

*/ lut_βize - entry*LUT_ IDTH*βizeof(ehort);

/*

** Write the lut file contents to the lut file

** with the following.

*/ βize_write -huge_f rite(lut, lut_βize,l,fhw);

If(size write—NODATA) {

(void)printf("Error writing lut file contentβ.\n") ; return(1); } doβ_βtatuβ-close(fhw);

;

/*

** Print a message indicating that procesβing is complete.

*/

(void) printf ( "normal termination -> look up tables are complete\n\n" ) ; return ( 0 ) ; >

/ ^• **** * * ******

* function: rnd_off ( )

* Round off a real number and assign to a short. * * * / short rnd off (double fl)

{ short i ; i« (fKO.O) 7 (short)ceil(f1-0.5) : (short)floor(f1+0.5) ; return (i) ;

}

/* ***** ***** * ------- ../ int nex km21atlon( coords *1 coordβ )

/*

* Name: nβx_tan21a lon

_ ~"

* Description:

* Thie routine converts the x and y distance of an "object" in

* the NEXRAD plane to the corresponding latitude and longitude

* values in degrees.

*

*/

/* — _—_■_—_,_,_,— — input parameters———————— */

/*

* phi radar deg "Radar" latitude (degreeβ) .

* lamda_radar_deg "Radar" longitude (degreeβ) .

* ^— ~"

* zx "Object" x diβtance from radar ( in km' β)

- 30 - where x increaβeβ towards the eaβt.

* zy "Object" y diβtance from radar (in km'β) where y increaβeβ towards the north. */

/* -Output para eterβ

/*

* phi deg "Object" latitude (degreeβ).

* lam3a_deg "Object" longitude (degreeβ).

*/ /*

{

/* — Internal Variableβ */ double phi; /* "Object" latitude (radianβ). */ double lamda; /• "Object" longitude (radianβ). •/ double phi_radar; /* "Radar" latitude (radianβ). */ double lamda_radar; /* "Radar" longitude (radianβ). */ double radiuβ - {6380.}; /* Average earth radiuβ (in kilometers) */ double constant - {135.}; /* Constant. */ double deg2radians - {M_PI/180. > /"Converts degreeβ to radianβ. */ double zp; /* Intermediate term. */ double zβin_the a; /* Intermediate term. */ double zcoe~theta; /• Intermediate term. */ double zβin~β; /* Intermediate term. */ double zcoβ~β; /* Intermediate term. */ double zβin~phi; /* Intermediate term. */ double zcoβ_phi; /* Intermediate term. */ double lamda delta; /* Intermediate term. */

/* */ double zx; double zy; double phi radar_dβg; double 1am3a_radar_deg; double phi deg; ~ double lam~a_deg;

/*

** Extract the radar latitude/longitude (in degreeβ) from the l_coords

** βtructure. ~

*/ phi radar_dβg - l_coordβ->phi_org la__3a_radar_deg -~l_coordβ->lamda_org

/*

** Extract the object poβition (in kilometreβ) from the l_coordβ

** structure.

*/ zx - l_coordβ->dx zy - l_coordβ->dy

/* Check the x and y diβtanceβ from the radar location. */ if( ((zx>-0.1)fifi(zx<0.1)) &β ((zy>-0.1)fi&(zy<0.1)) )

{

/*

* Thie handlββ the caβe when the x and y diβtancoβ are both

• near zero (the object iβ essentially at the radar).

* Set the object latitude and longitude (in degreeβ) to the radar

• latitude and longitude (in degreeβ).

- 31 - g; r deg; {

/*

* Thie handleβ the caβe when the x or y diβtanceβ are not

* near zero (the object iβ not at the radar). Compute

* the object latitude and longitude. First, convert the

* radar latitude and longitude from degreeβ to radianβ. */ phi radar > phi_radar_deg * deg2radianβ; lam3a_radar - lamda_radar_dβg * deg2radianβ;

/* compute the diβtance to the object (in kilometere). */ zp - βqrt( (zx*zx) + (zy*zy) );

/* Compute the βin and cos of the angle to the object. */ zβin_theta - (zx/zp); zcoβ_theta - (zy/zp);

/*

* Compute the intermediate βin and coβine termβ of the angular great

* circle diβtance from the radar poβition to the object poβition. */ zβin β - ( (zp/radiuβ) * ( 1.0-(( conβtant * zp )/ ~ (radiua • radiuβ))) ); zcoβ_β - βqrt( 1.0 - (zβin_β * zβin_β) );

/*

* Compute the intermediate βin and coβine termβ of the object

* latitude. */ zβin_phi - (sin(phi_radar)*zcoβ_β + coa(phi_radar)*zain_β*zeoβ_theta); zcoβ_phi « βqrt( 1.0 - (zβin_ρhi*zβin_phi) );

/* Compute the difference between object and radar longitude. */ lamda_delta - aβin( zβin_β • zβin_theta / zcoβ_phi );

/* Compute the object latitude and longitude (in radianβ). */ phi - atan2( zβin_phi,zcoβ_phi ); lamda - lamda_radar+lamda_delta;

/• Convert the latitude and longitude from radianβ to degreeβ. */ phi deg _* (phi/deg2radianβ); la__3a_dβg _* (lamda/dβg2radianβ);

.

/'

** Aeeign the latitude/longitude (in degreeβ) to the l_coordβ

•* etructure.

*/ l_coordβ->phi__value - phi__deg;

O 96/15504

- 32 - l_coordβ->lamda_value = lamda_deg;

/ ^•

** Return a zero value.

*/ return ( 0 ) ;

}

/_ _ _ _•_._ _ __, int nex_latlon2km( coords *l_coordβ )

/*

* Name: nex_latlon2km

*

* Deβcription:

* Thie routine converts the latitude and longitude of an "object" in

* the NEXRAD plane to the x and y componentβ of diβtance

* in kilometers.

/

{ - Input parameters double phi_deg; /* "Object" latitude (degreeβ). double lamda_deg; /* "Object" longitude (degreeβ) double phi_radar_dβg; /• "Radar" latitude (degreeβ). double lamda radar deg; /* "Radar" longitude (degreeβ).

/ ^*

/* -Output parameters */ double zx; /• "Object" x diβtance from radar (in km'β) */

/* where x increaβeβ towardβ the eaβt. •/ double zy; /* "Object" y diβtance from radar (in km'β) */

/* where y increaβes towardβ the north. */ * . /

I- Internal Variableβ */ double phi; /• "Object" latitude (radianβ). */ double lamda; /• "Object" longitude (radianβ). */ double phi radar; /* "Radar" latitude (radianβ). */ double lamda radar; /* "Radar" longitude (radianβ). */ double radiuβ - {6380.}; /* Average earth radiuβ (in kilometerβ). */ double constant - {135.}; /* Constant. */ double deg2radianβ - {M_PI/180.}; /"Convβrtβ degreeβ to radianβ. */ double zterma; /* Intermediate term A. */ double ztermb; /* Intermediate term B. */ double zβin β; /* Intermediate term: βin of angular great circle. •/ double zcoβ ^' β; /* Intermediate term: coβ of angular great circle. */ double zd; ~ /* Intermediate term D. */

/. / phi_dβg « l_coordβ->phi_valuβ; lam3a_dβg -~l_coordβ->l-__da_valuβ; phi radar__deg - l_coordβ->phi_org l__nda_radar_deg - l_coorda->lamda_org

/* Convert the radar latitude from degreee to radianβ. */ phi_radar - dβg2radianβ * phi_radar_dβg;

/* _C onvert the radar longitude from degreeβ to radianβ. */ lamda_radar - deg2radianβ * lamda_radar_dβg;

/* _C onvert the object latitude from degreeβ to radian ^β . */

phi - deg2radianβ * phi_deg;

/* Convert the object longitude from degreeβ to radianβ. */ lamda - deg2radianβ * lamda_deg;

/* Compute intermediate termβ a and b. */ zterma - cos(phi) * βin( lamda - lamda_radar ); ztβrmb - ( (cos(phi radar) * βin(phi))-

(βin(phi_radar) * coβ(phi) * coβ(lamda-lamda_radar) ) );

/* Compute intermediate βinβ and coβine of the angular great circle dietance from the radar poβition to the object poβition. */ zβin_β - βqrt( (zterma*zterma)+(ztermb*ztermb) ); zcoβ_β - aqrt( 1 - (zβin_β«zβin_β) );

/* Compute intermediate term D. */ zd - (constant * zsin_β) + radius;

/* Compute the x distance (in kilometers). */ zx - zd * cos(phi) * βin(lamda-lamda_radar) ;

/* Compute the y diβtance (in kilometerβ). */ zy - zd * (ein(phi) - (βin(phi_radar)*zcoβ_β) ) / coβ(phi_radar); l__coordβ->dx - zx; l ^" coordβ->dy « zy;

/*

** Return a zero value.

*/ return(0) ;

}

/***************** ** *** ^* ** int nex_latlon2bin( coordβ *1 coorde )

{ ehort bin_x, bin_y; int statue; ~ double zx,zy;

/*

•* Convert the lattitude/longitude to kilometerβ in the NEXRAD

** coordinate system with the following.

*/ status - nex_latlon2km( 1 coordβ ) ; if (status)

{ return(l);

> zx - l_coorde->dx zy - l_coorde->dy

/*

** Convert the x and y components of diβtance from kilometerβ to

** data bin coordinateβ. Round off the floating point valueβ aβ

** they are aββigned to integer variableβ.

*/

- 34 - bin x - rnd off( (zx+462. )/4.0 ); bin~y - rnd~off( (462.-zy)/4.0 );

/•

** Assign the x and y bin coordinates to the 1 coordβ βtructure.

*/ l_coordβ->x - bin_x; l~coordβ->y - bin~y;

/*

** Check that the x,y bin coordinate is within the valid range

** which determines the return value. A nonzero return value

** indicateβ that the latitude longitude corresponds to a

** nonvalid x,y bin coordinate.

^• / if( (bin x»MIN X)&&(bin x<PR BINS CR4)ft£ _{ Tbin-y>-MΪN-Y)SS(biE-y<PR-ROWS-CR4) )

/*

** This handleβ valid x,y bin valueβ. A zero value is returned.

*/ return(0); } return(2);

)

/******* ** * / int 1 c latlon2grid( coordβ *1 coordβ ) ehort grid_x, grid_y; int etatue7 ~ double zx,zy; double boxβize - {4.0};

/*

** Convert the incoming latitude longitude to a distance in

** (nominal) kilometerβ uβing the Lambert conformal projection.

*/ statue • 1 c latlon2km( 1 coordβ ) ; if(βtatuβ) ~

{

/*

*• If a nonvalid βtatuβ occurβ then return from thie routine

** with a nonzero βtatuβ.

*/ return(1);

/*

** Extract the x,y componentβ of diβtance from the l_coorda βtructure.

/ zx • l_coordβ->dx zy - l_coordβ->dy

/*

*• Convert the x,y componentβ of diβtance to grid x,y coordinateβ.

** The valueβ are rounded off aβ they are converted from floating

•* point to integer.

*/ grid x - rnd Off(zx/boxβizβ)+GRID X HALF 4KM; grid~y - GRlff_Y_HALF_4KM-rnd_off(zy7boxβIzβ);

/*

** Assign the grid x,y coordinates to the l_coords structure

** with the following.

*/ l_coordβ->x » grid_x; l_coordβ->y - grid_y;

/*

** Check that the grid x,y coordinate valueβ repreeent a valid

** array index. Nonvalid values result in a nonzero return.

*/ if ( (grid x>-MIN X) fifi(grid x<GRID X SIZE_4KM)&&

( grTd_y>-MΪN_Y ) && ( gr!d_y<GRΪD~Y_S I ZE_4RM ) )

/*

** Thie handleβ the valid grid x,y valueβ. A zero value is returned.

*/ return(0) ; > return(2) ; } * *** * * *** / int nex bin21atlon( coordβ *1 coordβ )

{ ehort bin_x, bin y; int etatuϊ; ~ double zx,zy; double range;

/*

** Extract the NEXRAD data bin coordinateβ from the 1 coordβ βtructure.

*/ bin_x - l_coordβ->x bin_y - l~coorda->y

/*

** Check the x,y coordinate pair for a valid array index value.

** A non valid coordinate reulte in a non zero return from the routine.

*/ if( (bin x>-MIN X)66(bin x<PR BINS CR4)&£(bin y»MIN Y)SS(bin y<PR_ROWS CR4) )

/* Convert the bin x,y coordinate to the x,y componentβ of diβtance

** where the reeult is a pair of floating point valueβ.

*/ zx - ((bin x*4.0)-462.0); zy - (-((bin_y*4.0)-462.0));

/*

** Use the componentβ of dietance to compute the range in kilometerβ .

*/ range - βqrt ( ( zx*zx) + ( zy*zy) ) ;

/*

** Check the range to inβure that it iβ in the NEXRAD coverage area.

** Othβrwiββ a nonzero value iβ returned.

*/ i (range<max range) {

/*

- 36 -

** Thie handleβ the coordinateβ located within the NEXRAD coverage ** area. Aββign the x,y components of distance to the 1 coordβ arra */ l_coordβ->dx « zx; l~coordβ->dy • zy;

/*

** Compute the latitude/longitude of the point with the following

** routine which convertβ x,y diβtanceβ in the NEXRAD plane to

** latitude longitude (in degreeβ).

*/ βtatuβ » nex_km21atlon( l_coordβ ); if(βtatuβ)

{

/* This handleβ a nonzero βtatuβ returned indicating non valid ** coordinateβ were paββed. Return a nonzero βtatuβ value. */ return(3) ₍

/*

** Thie indicateβ that valid latitude/longitude valueβ

** were computed. A zero value iβ returned.

*/ return(0);

> return(l); > return(2); )

/********** * * * * int 1 c grid21atlon( coordβ *1 coordβ ) { short grid_x, grid_y; int βtatuβ; ~ double zx,zy; double boxeize > {4.0};

/*

** Extract the grid x,y coordinate from the 1 coordβ βtructure.

*/ grid_x - l_coordβ->x grid~y - l~coordβ->y

/*

** Check that the x,y grid coordinate iβ a valid array index.

** Otherwise, a nonzero βtatuβ iβ returned.

*/ i ( (grid_x>-MIN X) fifi(grid x<GRID X SIZB_4KM) £fi

(grϊd_y>-MIN_Ϋ)££(grϊd~y<GRID_Y_SIZE_4KM) )

/*

*• Thie handleβ a valid x,y grid coordinate. Convert the x,y grid

** coordinate to x,y diβtance componentβ in (nominal) kilometers.

*• The y ordinate iβ reversed. Then aββign the diβtance components

*• to the 1 coordβ βtructure.

*/ zx - (grid x - GRID X HALF 4KM)*boxβize; zy - (GRID~Y_HALF_4~M~- grld_y)*boxβizβ;

- 37 - l__coordβ->dx - zx; l~coordβ->dy - zy;

/*

** Convert the x,y componentβ of diβtance (in kilometerβ) to

** latitude/longitude uβing the Lambert conformal projection

** with the following routine.

*/ atatuβ«l c km21atlon( 1 coordβ ) ; if(βtatuβ) ^"" ~

{

/* Thie handleβ the case when a nonvalid condition haβ occurred.

** A non zero value iβ returned.

*/ return(2); }

/*

** Thie handleβ the valid condition and a zero value iβ returned.

*/ return(0) ;

} return(1);

}

/************* * .. ..... ... .... .. ..... .. . ....

* function: l_c latlon2km( )

* Convert latitu3e/longitude to x/y uβing Lambert Conformal.

/* Convert latitude/longitude to x/y */ int l_c_latlon2km( coordβ *1 coordβ ) { etatic double deg2radianβ - {M_PI/180.}; /* Convertβ degreeβ to radianβ. •/ static double phi_l_deg - {33}; /* Standard latitude 1 (in degrees). •/ static double phi_2_deg - {45}; /* Standard latitude 2 (in degreeβ). •/ etatic double radiuβ - {6380.}; /* Average earth radiuβ (in km'β). */ etatic double phi 0 deg - {38.0}; /* Grid origin latitude (in degrees). */ static double lam3a_0_deg - {-98.0}; /• Grid origin longitude (in degrees). */

/*

** The following conetantβ were calculated for the grid center, earth radiuβ,

** and βtandard latitudeβ defined above. DO NOT uβe thie routine for other

** centerβ without recalculating the conβtantβ.

*/ static double phi 1- {0.575958653158129}; /* Standard latitudeβ (in radianβ). */ static double phi~2« {0.785398163397448}; /* Standard latitudes (in radianβ). */ etatic double phi 0- {0.663225115757845}; /* Grid origin latitude (in radianβ). */ etatic double lam3a_0- {-1.710422666954443}; /* Grid origin longitude (in radian ^β ). •/ etatic double rho 0- {7931.823624772932817}; static double n _* 70.630477697315427}; etatic double f- {1.955000201593793}; etatic double rf-{12472.901286168398656}; /• Intermediate term -(radiuβ*f) */ double phi deg; double lam3a_deg;

double zx; double zy; double rho; double theta; double phi; /* Current latitude value (radianβ). */ double lamda; /* Current longitude value (radianβ). */

/*

** Extract the latitude and longitude of the object (in degreeβ) from

** the 1 coσrdβ βtructure.

*/ phi_deg - l_coordβ->phi_value lamda_deg *~l_coordβ->lamda_value;

/*

*• Convert the latitude and longitude of the object from

*• degreeβ to radianβ.

*/ phi - phi_deg*deg2radiane; lamda - lamda_deg*deg2radianβ; rho - ( rf/(pow(tan(M_PI_4 + (phi/2)),n)) ); theta - n*(lamda - lamda_0);

/*

*• Compute the x and y componentβ of diβtance (nominally in kilometerβ).

*/ zx - (rho*βin(theta)); zy - (rho_0 - rho*coβ(theta) );

/*

** Aββign the x,y distance componentβ (nominally in kilometerβ)

** to the 1 coordβ βtructure.

*/ l_coordβ->dx • zx; l_coordβ->dy « zy;

/*

** Return a zero value.

*/ return(O);

) »... • • / int l_c_km21atlon( coordβ *l_coordβ )

Name: l_c_km21atlon

Deecription:

This routine convertβ the x and y diβtance componentβ of an "object' in the Lambert Conformal projection to the latitude and longitude in degreee.

/* —_._ Input parameters-- ——- */

/ _* -zx "Object" x distance from radar (in km'β) */

/• where x increaβeβ towardβ the eaβt. */

/* *

- 39 -

/* *«y "Object" y diβtance from radar (in km'β) */

/* where y increaβeβ towardβ the north. */

/* */

/* Output parameterβ— */

/* phi deg "Object" latitude (degreeβ). */ /* laπ~a deg "Object" longitude (degreeβ). */ / ^" */

etatic double phi 1- {0.575958653158129}; /* Standard latitudeβ (in radianβ) */ static double phi_2- {0.785398163397448}; /* Standard latitudeβ (in radianβ) */ etatic double phi 0- {0.663225115757845}; /* Grid origin latitude (in radianβ). •/ static double lam3a_0- {-1.710422666954443}; /* Grid origin longitude (in radianβ). * etatic double rho 0- {7931.823624772932817}; etatic double n- 70.630477697315427}; static double f _* {1.955000201593793}; etatic double rf-{12472.901286168398656}; /* Intermediate term -(radiuβ*f) static double radiuβ - {6380.}; /* Average earth radiuβ (in kilometerβ). */ static double deg2radianβ - {M PI/180.}; /*Conver β degreeβ to radianβ. V etatic double phi_l_deg - {33}; /* standard latitude 1 (in degreeβ) */ etatic double phi~2~deg - {45}; /* Standard latitude 2 (in degreeβ). */ etatic double phi 0_dβg - {38.0}; etatic double lam3a_0_dβg > {-98.0};

/* Internal Variablββ——— */ double phi; /* "Object" latitude (radianβ). */ double lamda; /* "Object" longitude (radians). double rho; double theta; double zx, zy; double phi deg; double laπκ_a_dβg;

/'

/*

*• Extract the object diβtance (in km'β) from the 1 coordβ βtructure.

*/ zx - l_coorde->dx/ zy - l~coordβ->dy

rho - eqrt( ( zx-zx) + ( (rho_0-zy)*(rho__0-zy)) );

/*

** Compute the latitude in radianβ.

*/ phi « {(2*atan( pow( (rf/rho) , (1/n) ))) - M_PI_2); theta - atan2( zx, (rho_0 - zy) );

/*

** Compute the longitude in radianβ.

*/ lamda - ( (thβta/n) + lamda_0);

/*

- 40 -

** Convert the latitude/longitude from radianβ to degreeβ. */ phi deg-phi/deg2radianβ; l-_n3a_deg-lamda/deg2radiana ;

/*

** λββign the latitude and longitude (in degreeβ) to the 1 coordβ βtructur

*/ l_coordβ->phi_value - phi_dβg; l_coordβ->l__nda_valuβ - lamda_deg;

/*

** Return a zero value.

*/ return(0);

>

/*

** Function to write contents of huge arrays to a file.

*/

/ * ** * * *****/ long huge_fwrite( void huge *buffer, long length, βize t count, int fHw) /*

** Deβcription: Function to write contentβ of huge arrayβ to a file. */ { char huge *pnt_buffβr; long remain; ~ βize_t amount; βizβ~t nwrite; long total; uneigned int doβ_βtatuβ; pnt_buffer - buffer; remain- length ; total-0; while ( remain>0 ) { if ( remain>MAXWRITE )

{ amount-MAXWRITB; remain — (long)amount; doe βtatue - _doβ_write(fhw,pnt_buffer,amount,finwrite); if(3oβ βtatuβj ~

{

(void)printf("doβ_write bad βtatuβ\n"); return(O); ~ > pnt_buffer +- amount; total +- nwrite; if(nwrite — 0)

{ return(total);

} } elee amount-(βizβ_t)remain;

- 41 - remain-0; doe βtatuβ - _doβ_write(fhw,pnt buffer,amount,&nwrite) ; if(3oe βtatuβ) ~

{

(void)printf("doe write bad βtatuβ\n"); return(0); ~ ) pnt_buffer +« amount; total +« (long)nwrite; if(nwrite — 0)

{ return(total) ;

> > > return(total) ; }

Figs. 9A, 9B, ?C, 9D and 9E together show a structure diagram or flowchart depicting a program capable of building mosaics of radar data. The program of Figs. 9A-9E is designed to run in real time by using the lookup tables built by the offline program of Figs. 8A-8E.

The program of Figs. 9A-9E, when started, runs as a continuous process. As shown in Fig. 9A, a manually created configuration file 200 containing a list of sites to be used in the mosaic is used to initialize the process at steps 201 and 202. First, the process reads in, from the configuration file 200, the number of sites (step 201) and the site list (step 202) that is then stored in memory for later use. Next, a loop is entered (step 203) that runs continually (shown in Figure 9B as a so-called "DO FOREVER" loop) . Each execution of the loop results in the creation of a new mosaic that is a compilation of the most current data from the individual sites specified in the site list. In preparation for the creation of the mosaic, all elements of the two-dimensional grid within the computer's memory are set (step 204) to the same predetermined value such as zero (i.e., the grid is cleared of any existing data) . Next, an index "to the current site in the site list" is set (step 205) to point to the first site on the list. Then, another loop is entered that processes (step 206) the radar sites on the site list one

at a time until all of the sites have been handled (st 207) .

The processing of each radar site begins, as shown Figure 9B, by extracting the site ID number from the si list (step 220) and then incrementing the site index by o (step 221) . The site ID number is used by the program locate (step 222) data for the current site within t computer's local database 46. If no data exists for t current site (step 223) then the program skips to the ne site. Otherwise, the program identifies the most curre data for the site. This data (called a product) resides a file (radar product data 225) in the local database 4 First, the file is opened (step 224) and the produ "header" is read (step 226) into the program's memory. T header contains information such as data date and time a product data length. If, by examining the product date a time, the data is determined to be too old (step 227) , th the program skips that site and goes on to the next sit Otherwise, using the data length, the remainder of the fi 225 is read (step 228) into the computer's memory, aft which the product file 225 is closed (step 229) .

Next, as shown by Figure 9C, the product data decoded and stored (step 241) in a two-dimensional "raste array. If any errors are detected during the decodi process, then the program skips to the next site (st 242) . When the product is successfully decoded, t program uses the current radar site ID number to determi the location of the appropriate look-up table file in t computer's local database 46 (step 243) . First, the looku table file 245 is opened (step 244) and the heade (containing the number of table entries) is read (step 246 into the program's memory. The program then uses th number of table entries to read in (step 247) the complet look-up table from file 245, which table is then stored i the program's memory. The look-up table file 245 i subsequently closed (step 248) .

- 43 -

Figure 9D shows how, once the look-up table is read in, the single site radar data (which resides in the product raster array) is assigned to the two-dimensional grid. A loop is entered to process all of the entries in the look-up table (step 260) . When all processing is complete, the program goes on to the next site (step 261) . First, the x and y indices into the two-dimensional product raster array are extracted from the look-up table (step 262) and checked for validity (step 263). (An invalid entry should never occur in practice, but if it should occur the program skips that entry and goes to the next entry.) The product data value is then extracted (step 264) from the product data array. Then, the x and y indices into the two-dimensional mosaic grid are extracted from the look-up table (step 265) . Those coordinates are also checked for validity (step 266) . (Again, an invalid entry should never occur in practice, but if one should be detected then the program skips to the next entry.) Next, the product data value is checked for validity (step 267) and any invalid data values cause the program to skip to the next look-up table entry. Then, the product value is compared (step 268) with the existing grid value and, if greater, is assigned (step 269) to the grid at the x,y location, and, if not greater, is skipped. Processing continues in that way for all remaining entries in that look-up table.

After all available sites have been so processed, the data in the grid is converted to a new product as shown by Figure 9E. First, the values in the grid are converted (step 280) to the output product format and assigned to a buffer in the computer's memory. Then, a header is built (step 281) , also in the computer's memory, containing the current date and time as well as the product data length. A unique output filename is created (step 282) and a file of that name is opened for output (step 283) . Next, the data is actually written to the output product file 284. First, the product header is written (step 285) followed by

- 44 - the product data itself (step 286) . Finally, the out product file 284 is closed (step 287) and the program g on to start the next mosaic.

Direct mapping of individual radar data values t mosaic grid, using one or more prebuilt lookup tab constructed with the software given above, can implemented in software. An example of such softwa written in the ANSI version of the C language usable o Sun Microsystems workstation is given in the follow pages. However, other programming languages can be u for this purpose. Figs. 9A-9E together generally sho flowchart for this software. The software appearing in following pages implements steps 48 and 50 of Fig. 4. software shown above, and the software shown bel together implement steps 34, 36, 38 and 40 of Fig. 3

----•---____-_-_**-_•_ ___.

** UNISYS WIS MOSAIC PROCESSING SYSTEM *• Copyrigth 1994 Uniβyβ Corporation

**

** PROGRAM NAME : BUILD . EXE

- •

** DESCRIPTION: Read in individual radar productβ and corresponding ** Look Up Tableβ and map the product data to a grid.

** Then take the reβulting grid valueβ and produce a

** new product.

»»

** ASSUMPTIONS: Productβ for the desired βiteβ reβide in the βpecifiβd ** directoriββ. Look Up tableβ for the deβired βitββ exiβt

** in the "lut" directory.

*•

*»_-_-**_-_»*___#*_»___**_*_*_***-***_***_-**_-**

/include <raath.h> /include <errno.h> /include <time.h> /include <etdio.h> /include <βtdlib.h> /include <βtring.h> /include <ctyρe.h>

/*

** Prototype of eubroutine.

*/ int gβt_lateβt(ehort , ehort , char *);

/define BIN X IDX 0 /define BIN~Y~IDX 1 /de ine GRID X IDX 2 /define GRID~Y~IDX 3 /define LUT_WIDTH 4

/define MAX DATA SIZE 300000

/define NUM_DATλ_LEVELS 16

/define NODATA 0 /define MIN NUM 0 /define MIN ^' VALUE 0 /define NO VALUE 0 /define Ml X 0 /define MINJf 0

/define GRID X SIZE 1280 /define GRID~Y~SIZE 896

/define LUT MAX 50000 /define PR BINS 232 /define PR~ROWS 232 long lut_key ■ {123456}; typedef struct{ long keyword; long num_entrieβ; short header_βize; short βite id; ehort version; short date_built; unβigned long time_built; short grid_reβ; ~ short grid~width; short grid_height; short bin_reβ; short bin_width; short bin_height; short projection; double grid_centβr_phi; double grid_center_lamda; double radar_phi; double radar_lamda; } lut_hβader_format; lut_header_format lut_headβr; typedef struct{ double dx; double dy; double phi org; double lax_3a_org; double phi value; double lam3a_value; ehort x; ~ ehort y; ehort pr id; ehort grld_id; } coordβ; long aite cey - {24680}; short nβxrad_βitβ « {1}; ehort rrwdβ_βite - {2}; short yee -~{1}; short no - {0}; typedef struct{ long keyword; ehort num βiteβ; } βitβ_hβadβr_format; βitβ__hβedβr_format βitβ_filβ_headβr; typedef etruct{

- 46 - ehort aite_id; ehort βite~flag; ehort βite_build_lut; ehort βite~uβe; ~ double aite_lat; double βite~lon; }aite__contentβ_7ormat;

/define MAX_SITES 200 βitβ_contentβ_format βitβ_file_contentβ[MAX_SITBS];

FILE * fpw; FILE * fpr; FILE * fpr_prod;

static int errno - 0;

}βitβ_path_format;

/define NUM_SITB_PATH 30

/*

** Example βite path information.

*/ etatic βitβ_path_format βitβ_path(NUM_SITE_PATH] -

354,-AL_maxwβll AFB/", 395, -AR_littlerόck/", 524 , "ΛZ_phoenix/ " , 347, "CO denvβr/", 351,"DE~dovβr/", 302 , "PL~melbourne/ " , 728,-_T_~mi__ni/", 337,-IL~chicago/", 381, -IN~indianapoliβ/", 350, -KS_dodgβcity/ " , 366, "KS goodland/", 554,"KS~topβka/", 562 , -__S~wichita/" , 349, -MI~dβtroit/", 308, -MO_βtlθuiβ/", 342, "MS columbuβ/", 382 , -MS ^" j ackβon/ " , 367, "N_f ^" haβtingβ/", 372," Y~gri fiβ AFB/", 340 , -OH~clβvβland/ " , 305 , -OK~f rβdβrick/ " , 301 , "OK~norman/ " , 557,"OK~tulea/*\ 525 , "PAjitteburgh/ " , 313 , "TX_amarillo/ " , 332, "TX central/", 353,"TX~dyβee AFB/", 378, "TX ^"" houeton/", 303 , "VA~βterling/ " , 000, "UNDEFINED/", )/

- 47 - typedef struct{ short prod_id; char *prod~text; }prod_jpath_format;

/define NUM_PROD_PATH 1 static prod__path_format prod path[NUM PROD PATH]

{

38,"COMP REF/",

>

/ ^•

** Main routine to build mosaic. */ main()

long entry; long nuπ_entry; ehort prβv_col; ehort nbytee; long layer length; long layer length; long block_Tength; long meββage_length; long total_lenςrth; ehort βcale_x, βcale y; uneigned ehort *nbytee_indβx; uneigned ehort *mβw_layer len_index; unβigned ehort *lβw~layer~len_indβx; unβigned ehort *mβw_layβr2_len_index; uneigned ehort *lβw~layer2~len_index; unβigned βhαrt *mβw block len_index; uneigned ehort *lβw_block_len_indβx; uneigned ehort *mβw_mβg_len_indβx; unβigned ehort *lβw~msg~len~index; uneigned ehort mβw_layβr_lβn; uneigned ehort lβw_layer_lβn; uneigned ehort maw_layer2_len; uneigned ehort lβw_layβr2~len; uneigned ehort mβw_block_len; uneigned short lβw~block_len; uneigned ehort msw~mβg_len; uneigned short lβw_mβg_lβn; ehort ref; ~ ~ ehort dl;

/*

** Data level codee to aββign to output product header.

*/ etatic unβigned ehort data_level 16[NUM_DATA LEVELS] - {0x8002,0x0005,0x000a, OxOOOf,0x0014,0x0019,OxOOle,0x0023,0x0028,0x002d, 0x0032,0x0037,0x003c, 0x0041,0x0046,0x004b }; ehort pr_typβ; uneigned ehort product color; unβigned ehort grid_coTor; ehort flag_off - {0}; ehort flag_on - {1}; ehort num; ehort prid; ehort dββtination_id; ehort num blockβ;~ short divider; long latitude origin; long longitudϊ origin; double phi center deg - {38.0}; double lam3a_cβntβr_dβg - {-98.0}; double lat_origin_dβg; double lon~origin_dβg; uneigned βEort mβw_latitudβ; unβigned ehort lew latitude;

- 49 - uneigned ehort mew_longitude; uneigned ehort lβw_longitudβ; ehort hβight_valuβ7 ehort eequence_num; ehort βcan_num7 ehort βlβvation_num; ehort βym_block~id; ehort num_layerβ; uneigned ehort raβter_opcodβ; uneigned ehort raβtβr~conβt1; uneigned ehort raβtβr_conet2; uneigned ehort raβter_conet ; char *out βtr home - "/home/"; char *out~etr~βitβ - "NATIONAL/"; ehort num_includβd; ehort num_excluded; ehort n; ~ ehort tab included[MAX_SITES]; ehort tab~excluded[MAX_SITES] ; ehort year, onth,day; long time_eec; long eecoήde; long cur_ββcondβ; ehort hour,minute; ehort hh,mm; long ee; βtruct tm Jan_01_1970_βtruct - {0}; /* Initialize all fieldβ to 0. ; timβ_t total_timβ; βtruct tm gmt_βtruct; ehort j_datβ; short cur_j_datβ; short mβg~date; long meg_time; uneigned~βhort mβw_mβg_timβ; uneigned ehort lβw_mβg_timβ; ehort ecan_date; long ecan__time; unsigned ehort mew_ecan_time; uneigned ehort lβw_βcan~timβ; ehort gβneration_datβ; long generation_time; uneigned ehort mβw_genβration_timβ; uneigned ehort law~gβnβration_timβ; union{ long l_valuβ; ehort β_valuβ[2]; } long2ehort;

/*

** Obtain a buffer into which βxiβting product ^β are read

** and new producte are built and then etored from.

*/ buf-(cbJ_r *)malloc(MAX_DATA SIZE); if( buf—NULL)

- 50 -

{

(void)printf("error: malloc failed for buf, errno«%d\n",errno); return(l); >

/*

** Open the configuration file and read in the header.

*/ fpr ^» fopen( "conf ig.dat" , "rb" ) ; if(fpr—NULL)

{

( oid)printf("fopen error for configuration file header, errno-%d\n",er return(1); } nread-fread(fieite_file header,βizeof(βite file header),1,fpr); if(nread—NODATA) ~ ~

{

(void)printf("fread error with βitβ file header, errno-%d\n",errno); return(1); >

/*

** Check the configuration file header for the correct keyword

** and exit the program if the keyword doee not agree.

if( eite file header.keyword 1- eite key ) {

(void)printf("Invalid configuration file keyword\n"); rβturn(l); >

/*

** Obtain the number of βiteβ from the configuration file header and

** uβe that number to determine the βize of the configuration file

** contentβ. Reed the configuration contentβ from the configuration

** file. That provideβ a liβt of βiteβ and indicateβ which are NEXRAD

** and which are active (to be uββd in the mosaic).

V num_βiteβ - βite_file_hβader. num_βitββ; nread-fread(eite file contentβ, num βiteβ*βizeof (βite contentβ format) , 1, fpr) if ( nread— NODATA ^"

<

(void)printf("error reading configuration file contentβ, errno-%d\n",err return(1);

> βtatuβ - fcloee(fpr); if( βtatuβ 1- 0 )

{

(void)printf("error cloβing configuration file, errno-ld\n",errno);

} while(l) /* DO FOREVER */ {

/*

** Obtain the current time in βeconde βince 1/1/1970 which

•* le ueed to determine the age of exlβting productβ

** aβ they are read in from the databaee.

*/

Jan_01_1970_βtruct.tm_year - 70; Jan 0l l970~βtruct.tm mon - 0; Jan~01 1970 βtruct. tπfmday - 0;

Jan_01_1970_t - mktimβ(£Jan_01_1970_βtruct) ; if (Jan_01_1970 t — (time t)-l) {

(void)printf ( "error: mktime(l/l/70) error. \n" ) ; return(l) ; } total_ti_ne-tiroe(NULL); gmt_βtruct - *gmtimβ(fitotal_time) ; year ■ gmt_βtruct.tm year; if(year>99)

< year—100; } month - gmt_βtruct.tm mon + 1; day « gmt βtruct.tm_m3ay; hour - gn-__βtruct.tm_hour; minute - gmt_βtruct.tm_min;

(void)printf("Current (gmt) year-%d month-%d day-%d hour-%d min-%d\n", year,month,day,hour,minute) ; gmt_t « mktime(£gmt βtruct); if (gmt t — (time t)-l) {

(void)printf("error: mktime(gmt) error.\n"); return(l); elβe

< ti__e_ββc - difftime( gmt_t, Jan_01_1970_t ); cur_j_datβ - time_βec/86400; cur eeconde - time βec%86400; }

/*

** Initialize the counterβ to track βitee included in and excluded from the moβaic.

*/ num_includβd - 0; num_excluded ■ 0;

/*

** Clear the grid in which the moβaic iβ to be built.

*/ for( grid_y-MIN_Y; grid_y<GRID_Y_SIZE; grid_y++ ) for( grid x-MIN_X; grid x<GRID_X SIZE; grid_x++ )

{ grid[grid x][grid y] - NO VALUE;

} )

/*

** The following ie the main loop to read in individual NBXRAD

** producte and map the data to a large grid. Loop through

** the eitee contained in the configuration βtructurβ.

*/ for (j-MIN_-TUM; j<num βiteβ; j++)

<

/•

•* Extract the βite id to requeet from the configuration file

** data. */ βite_id - βite_file_contentβ(j].βite_id;

/*

** Process NEXRAD data here. Requeet compoβite reflectivity.

*/ prod_id - 38;

/*

** Find a path to the lateβt product of the βpecified type

** for the specified βite with the following. A nonzero

** βtatuβ code indicateβ that no product iβ available

** BO jump to the next.

*/ βtatuβ - get_lateβt( eite_id,ρrod_id,prod_narae); prod_name[strlen(prod name)-1]-'\0' ; if(status) ~

{

(void)printf(" product is unavailable.\n"); num_excluded++; tab_excluded[num_excluded-l] - βite_id; goto jump here; ~ ~

}

/*

** Open the product file uβing the path to the lateβt.

*/ fpr_prod-fopen(prod_name,"rb");

;

/*

** Read in the product header.

*/ nread-fread(buf,header size,1,fpr_prod); if( nread—NODATA )

{

(void)printf("product read error, errno-%d\n",errno); num_excluded++; tab~excluded(num excluded-1] * βite_id; printf(" path-%β\n",prod_name); goto jump here; )

/*

** Extract data from the product header. Firβt poβition the start

** index variable at the beginning of the buffer. */ et_indβx-(unβigned ehort *)buf; index-βt_index; meg_code«*index;

** Extract the operational mode and volume coverage pattern from

** the product header. Use thie information to check if the radar

** wae in a valid mode. If the radar iβ in maintenance mode

** then βkip the product and jump (goto) to the end of the loop. */ index-βt_index+16; op_m'odβ-~indβx; index-βt_indβx+17; vcp-*index; if ( (op_model-l)fifi(op_model-2)) (void)printf(" invalid op_mode-%d vcp-%d\n", ' op_mode,vcp); ~ num_excludβd++; ~ tab excluded[num_excluded-l] - βite_id; printf(" path-%β\n",prod_namβ); goto jump here; }

/*

** Extract the input product βite id number. Compare with the requeβtβd

** eite id number. If the extracted and requeeted eite id'e do not

** match report an error and jump (goto) the next βite'β product.

*/ index-βt_indβx+6; input_βoϋrcβ_id-*index; if (eite file contentβ[j].βite_.id 1- input_βourcβ id) %d) site miβmatch\n",

/*

•* Extract the date (in dayβ βince 1/1/70) and the time (eeconde)

** from the product header. Uee this to determine the total time

** in eeconde since Jan 1, 1970. Thie iβ uβed in a compariβon with the

** current time to determine the "age" of the product.

*/ index-βt index+20; j_date-"Index; index-βt_indβx+21; timβl-*indβx; indβx-βt_indβx+22; time2»*index; tot_aβc«(j_datβ*86400)+(ti__βl*65536)+timβ2;

( age eeconde - timβ_ββc - tot_ββc; agβ ninutββ - agβ_ϊβcondβ/607

>

/*

** Check the age of the product. If the product iβ too old then

** jump to the next βite.

*/ if (age minutes > max age minutes)

/*

** Determine the length of the data layer (in bytes).

*/ , index-βt_index+66; lengthl-*indβx; indβx-βt_indβx+67; length2- ^" index; data_βizβ-( long )( (lengthl*65536) + lengths );

to read");

/*

** Read in the reβt of the product with the following.

*/ nread-fread(buf,data βize,l,fpr_prod); if( nread—NODATA )

{

(void)printf("error reading product file, errno-%d\n",errno); num_βxcludβd++; t_JO_excluded(num_exeluded-l] ■ βitβ_id; print (" path-%β\n",prod_namβ); goto jump here; ~

} βtatue - fcloββ(fpr_prod); if( etatuβ 1- 0 )

{

(void)printf("error cloβing input product file, errno-%d\n",errn return(1);

}

/*

** Locate the poβition of the run length encoded raeter data

*• and extract the opcode and check it for validity.

*/ at index-(uneigned ehort *)buf; in3βχ-et index;

opcode-* index; if(opcodeI-0xba07) {

(void)printf("error: invalid raeter opcode in input product«%x \n",opcode) ; num_excluded++; tab~excluded[num excluded-1] - βite id; printf(" path-%e\n" ,prod_namβ)7 goto jump here; ~

>

/*

** Extract the number of rowβ from the product data layer.

*/ indβx+-9; nrowβ-*index;

/*

** Point to the firet row of run length encoded raeter data.

*/ index-*—2;

/*

** Check for a valid number of rowe.

*/ i ((nrowβ>MIN )fifi(nrowe<-PR ROWS)) {

/*

** Loop through the rowe of raeter data.

*/ for(row»MIN_Y;rowcnrowe;row++)

{ n bytββ-*index; index-*-*-; cp-( char *)index; bin « NO_VλLUE;

/*

** Loop through the byteβ of RLE data in the current row.

*/ for(m-1;m<«n_byteβ;m++)

/*

*• Extract the current byte and determine the run and color.

*/ byte-*( cp ); cp++; run-(( ehort )byte) / 16; col-(( ehort )byte) % 16; if(run>0) for(k-MIN Y;k<runk++)

{ if((bin<PR_BINS)£6(row<PR_ROWS)) product[bin] [row]-col; elβe

{

/*

** In the unlikely event that the bin or row numb ** exceedβ the array limit then report an error a ** jump to the next site's product. */

(void)printf("error: bin-%d row-%d\n",bin,row) ; nuα_excluded++; tab_ ^" excludedfnum_excluded-l] « sits id; printf(" path ^« % ^β \n",prod_namβ)7 goto jump here; ~

) bin++;

} > > index-(unβigned short *)cp; elβe of rowβ-%d",nrowβ);

/*

** Next, open and read in the look up table file. Firet, create

** a etring containing the look up table (lut) file name (which iβ

** iβ baaed on the βite number) and open the lut file for "read binary".

*/

(void) ^β printf(lut_name,"lut/cr4_%3.3d.lut\0",βitβ_file_cont βntβ[ j.βite_id fpr-fopen(lut name,"rb"); if(fpr—NULL) ^"

{

(void)printf("error: error opening lut file (%β) for read, errno-%d\n" lut_name,errno) ; num_excluded++; tab~excluded[num_excluded-l] »βitβ_id; goto jump here; }

/*

** Next, read in the look up table (lut) file header from which the

** number of entries iβ extracted which iβ uβed to compute the

** βize of the lut file contentβ (lut βize). Alβo, check the

** keyword and the header aize to verify that a valid "lut"

** header wae read in.

*/ nread - frβad(£lut_header, βizeof ( lut_header_ ormat ) , 1, pr) ; if ( nread— 0 ) ~ ~ ;

num entry - lut header.num_entrieβ; i (ήum entry > LUT_MAX)

{

num_entry - LUT MAX; > if( lut header.keyword 1- lut_key ) {

(void)printf("invalid lut keyword * %d\n",lut header.keyword) ; num_entry - 0; ~

} if( lut header.header βize I* βizeof(lut header format) ) { . ^{" "}

(void)printf("invalid lut header βize - %d\n",lut_header.header_βize) ; num_entry - 0; ~

} lut_βize - num_entry*LUT_WIDTH*βizeof(ehort);

/*

** Do the read to read in the lut ltβelf.

*/ nread - fread(lut,lut_aize,l,fpr) ;

/*

** Now that the lut haβ been read in, cloβe the file deβcriptor.

*/ βtatue - fcloβe(fpr); if( βtatue 1- 0 ) {

(void)printf("error: error cloβing look up table file, errno-%d\n",errno); return(1) ; }

/*

** Loop through all of the entrieβ in the look up table

** and map radar data bine to grid boxeβ.

*/ for(entry-MIN_NUM; entry<num_βntry; entry-*-*-) {

/*

** Extract the current radar product'a bin and row coordinate.

*/ bin - lut(entry][BIN_X IDX]; row - lut[βntry][BIN_Y~IDX];

/*

** Check that the bin and row are valid indicieβ in the

*• product data array (that iβ, they are valid).

*/ if((bin>-MIN X)£fi(bin<PR BINS)fi£(row»MIH_Y)fi£(row<PR ROWS))

{

/*

** Extract the product color (data level) value from the

** product array for the current bin and row.

*/ product_color - product[bin)[row);

/*

** Extract the corresponding x and y grid coordinateβ.

*/ grid x - lut[entry][GRID_X IDX]; grid^y - lut[entry] [GRID_Y~IDX];

/*

** Check the current x and y grid coordinateβ to verify that

** the valueβ are valid array indicieβ.

*/ if((grid x>*MIN X)£fi(grid x<GRID X SIZE)fifi

(gri3_y>-MIN_Y7fifi(gri3_y<GRID_Y_SIZE) )

/*

** Examine the exiβting grid box for a value.

*/ if( grid(grid x](grid_y] — NO VALUE )

( grid[grid_x][grid y] - product color; elβe

{

/*

** In the caβe with an exiβting value in the grid the

** a compariβon iβ needed between the exiβting grid

** value and the current product color value.

*• Map the product color when it iβ greater than the

** exiating grid color value.

*/ grid_color - grid[grid_xj[grid_yj; if(product color > grid color)

{ grid[grid_x][grid_y] - product_color;

} eli

{

/*

** Thie handleβ the unlikely caβe that an x,y grid coordina

** wae not a valid array index. For now, print the value

•* and go to the next lut table entry.

*/

(void)printf("entry-%d grid_x-%d grid_y-%d\n",entry,grid_x,

> elee

{

/*

** Thie handleβ the unlikely caβe that a product bin row and co

*• waβ not a valid array index. For now, print the value

•* and go to the next lut table entry.

*/

(void) rintf("entry»%d bin»%d row»%d\n",entry,bin,row);

}

jump here:; /* Jump here to βkip a product and go on to the next βite.

} /•

96/15504

- 59 -

** Now that the loop to proceββ the individual eite productβ

** into the grid iβ complete, print out the βiteβ which were included

** and the βitee which were excluded.

*/

(void)printf("Siteβ proceββing complete:\n");

(void)printf(" number excluded-%d\n",num_excluded); for(n-0; n<num excluded; n++) ~

{

(void)printf(" excluded-%d\n",tab excluded[n]); )

/*

** The producte have been mapped to the grid. Build a moβaic product hart

*/ prid - 5;

/*

** Poβition the pointer at the beginning of the buffer in

** order to be ready to aββign data to the buffer.

*/ βt_indβx-(uneigned ehort *)buf; index>et_index+0;

/*

** Build a etring containing the output filename.

*/

(void)atrcpy(prod name,out βtr home); (void)βtreat(prod~namβ,out~βtr~βite);

/*

** Continue to build the product pathname. The following buildβ

** the actual file name in year, month,day,hour,minute format.

** Concatenate the filename to complete the product pathname.

*/

(void)βprintf(prod_file,"%.2d%.2d%.2d.%.2d%.2d\0", year,month,day,hour,minute); (void)βtrcat(prod_namβ ,prod_file);

/*

** Open the product pathname for write binary with the

* ^■ following. */ fpw-fopen(prod namβ,"wb"); if(fpw—NULL)

(void)printf("fopen for write error, return(l); >

/*

** Define the paramβtβrβ which epecify how the product

** ie to be built.

*/

number_rowβ _* 896; number bine - 1280; βcalβ_x - 1; βcalβ~y - 1; etart~row - 0; βtart~bin - 0; output_βourcβ_id - 10000; lat_origin_deg - phi_center_deg; lon_origin~deg - lamda_center_deg;

/*

** Then βtart aββigning the data to the product header. Every field

** muet get a value to overwrite any exiβting data in the buffer.

*/ meg code « prid; /* Aββign the meββage code */

*in3ex - mβg_codβ; /* which iβ the product ID . */

/* λββign current date in dayβ since 1/1/70 meg_ti__β - cur_βecondβ; /* Convert seconds to 2 halfwords. * long2βhort.l_value - mβg_time; index-et_index+2; mew mβg_ti__β - long2βhort.β_value(0];

*in3ex - mβw_mβg_time; ~ /* Aββign current time (mew). */ indβx-βt_indβx+3; lew mβg_time - long2ehort.β_value(l];

*in3βx - lβw_mβg_timβ; ~ /* Aββign current time (lew). */ mβw_mβg_len_indβx-βt_indβx+ ; /*Save poβition of the */ lβw~mβg_len_indβx-βt_indβx+5; /* of the message length. */ index-βt_index+6;

•index -~output_βource_id; /* λββign βource id. •/ index-βt index+7; deβtinatIon_id - NO_VALUE;

•index - deϊtination_id; /* Aββign deβtination id (unuβed) . •/ index-et indβx+β; num bloclβ - 3;

*in3ex - num_blockβ; /* Aββign the number of blockβ. */ indβx-βt_indβx+9 ; divider - -1;

*index-divider; /* Aββign block divider. */ latitudβ_origin - (long)(lat_origin_deg*1000.); long2ehort.l_valuβ - latitude_origin; indβx-βt_index+10; raβw_latitudβ - long2βhort.β_value[0];

♦index - raβw_latitude; /* Aββign origin latitude (mew). •/ index-et indβx+11; lew latitude - long2ehort.β valuβ(l);

•in3ex - lβw_latitude; /* Aββign origin latitude (lew). */ longitude origin - (long)(lon_origin_dβg*1000.); long2ehort.l_valuβ - longitudβ_origiή; index-et indβx+12;

/15504

- 61 - mew longitude « long2ehort.β_valuβ[0];

*in3ex - mew_longitude; /* AΪβign origin longitude (mew). */ index-βt_index+13; lew longitude « long2βhort.β_valuβ(l];

*in3ex > lβw_longitude; /• Aββign origin longitude (lew). */ index-βt index+14; hβight_valuβ - NO_VALUE;

*index ^" - hβight_value; /* λββign the height value (unuβed).*/ indβx-βt_index+15;

•index -~prid; /* λββign the product ID. */

•index - op_mode; /* λββign the operational mode value. */ index-βt index+17; vcp - NO~VλLUE;

•index -~vcp; /* λββign the volume coverage pattern. */ index-βt_indβx+18; eequence~num ■ NO_VλLUB;

•index «~βequence~num; /* Aββign the sequence number (unueed). •/ indβx-et_index+19; βcan_num~- NO_VALUB;

•index - βcan_num; /* λβeign the volume βcan number (unueed). */ index-et_indβx+20; βcan_datϊ - cur_j_datβ;

*index - βcan_date; /* λeeign the βcan date. */ βcan time - cur eeconde; long_ ^" βhort.l_vaTuβ - ecan_timβ; indβx-βt_index+21; raβw_βcan~timβ - long2βhort.β_value[0);

•index -~mew_ecan_time; /* Aββign the ecan time (mew). */ index-βt_indβx+22; lew ecan~timβ - long2βhort.β_valuβ[l];

•index -~lβw_βcan_timβ; /• λeeign the βcan time (lew). */ index-st_index-*-23; generation_date - cur_j_date;

*index - genβration_datβ; /* Aββign the generation date. */ generation_time - cur_βecondβ; long2ehort l_valuβ > generation_time; indβx-βt_indβx+24; mew generation time « long2ehort.e valuβfO]; *in3βx - mew generation time; /* λββign the generation */ ~ ^~ /* time (mew).*/ index-βt index+25; lew generation time - long2ehort.β valuβ[l]; *in3βx - lew glneration time; /* λϊaign the generation */

/* time (lew).*/ index-βt indβx+26;

*index -~NO_VλLUB; /* Product Dependent (unuβed). */ index-βt indβx+27;

*index -~NO VALUE; /* Product Dependent (unueed). •/

index*et_index+28; elβvatiort_num - NO VλLUE;

•index > elevation ^' num;/* λββign the elevation number (unuβed). */ index-βt index+29; •index -~NO_VλLUE; /* Product Dependent (unuβed). •/ index-βt index+30/ /* Save poβition of data levels. */

** λββign the appropriate βet of data levels valueβ.

for( dl-MIN_VλLUE; dl<NUM_DλTλ_LEVELS; dl++ )

{

•index - data_lβvel_16[dlj; index-*--*-; ~ ~

> index-βt index+46; •index «~NO_VλLUE; /• Product Dependent (unueed). */ index-βt indβx+47; •index »~NO_VALUE; /* Product Dependent (unuβed). */ index-βt index+4β; •index -~NO_V___UB; /* Product Dependent (unuβed). •/ index-βt indβx+49; •index -~NO_VλLUB; /* Product Dependent (unuβed). */ index-βt indβx+50; •index -~NO_VλLUB; /• Product Dependent (unueed). */ index-βt indβx+51; •index -~NO_VλLUE; /* Product Dependent (unueed). •/ index-et_index+52; •index -~number_binβ; /• Number of columns (bins). •/ index-et indβx+53; •index -~NO_VλLUE; /* Unuβed. */ index-βt index-*-54; *index «~NO_VλLUB; /* Offeet to eymbology (mew) . */ indβx-βt_indβx+55; *index -~60; /* Offeet to eymbology (lew). */ index-βt index+56; ♦index -~NO_VλLUB; /* Offeet to graphic (mew) . */ index-βt indβx+57; ♦index -~NO_VλLUB; /* Offeet to graphic (lew). */ index-βt indβx+58; •index -~NO_VALUE; /* Offeet to tabular (mew). */ index-βt index+59; •index -~NO_VλLUB; /* Offeet to tabular (lew). */ indβx-βt_indβx+60; /* Beginning of eymbology block */ divider - -1; *index-divider; /• Aββign block divider. ^♦ / index - βt indβx+61;

- 63 - βym block_id « 1;

*in3ex-eym_block_id; /• λββign eymbology block ID. */ index- βt_index+62; mβw_block_len_index * index; /* Save location of block βize value*/ index- βt_index+63; lew_block_len_index » index; index- βt_index+64; num layere - 1;

•in3ex - num_layerβ; /* Aββign number of layere */ index- et index+65;

*index - 3ivider; /• λββign layer divider. */ index- βt_indβx+66; mew_layer_len_indβx « index; /* Save location of layer length. */ index- ,et_index+67; lβw_layer_len_index > index; index- βt_index+6β; raeter_opcodβ - 0xba07;

*index-- raeter_opcode; /* Aβeign Raeter data format opcode. */ index- βt_index+69; raβtβr_conβtl - 0x8000;

*index-- raβtβr_conβtl; /* λβeign Raeter data format conetant. */ index- et_indβx+70; raβter_conet2 - OxOOcO;

•index-- raβter_conβt2; /* λββign Raeter data format conetant. */ index- βt_index+71;

•index - NO_VALUB; /♦ λββign placeholder for raeter x origin. ♦/ index- et indβx+72;

♦index - NO_VALUB; /♦ λββign placeholder for raeter y origin. ♦/ index- βt indβx+73;

•index - T; /♦ λββign ecale x (integer portion). •/ index- βt_indβx+74;

♦index - NO_VλLUB; /♦ λββign scale x (fractional portion). ^♦ / index- βt index+75;

♦index - T; /♦ λββign βcale y (integer portion). ^♦ / index- βt_indβx+76;

♦index - NO_VALUE; /♦ λββign βcale y (fractional portion). ^♦ / index- βt_indβx+77;

♦index - numbβr_rowβ; /♦ λββign number of rowe. •/ index- βt_index+78; raeter conβt3 - 2;

•index"- raβter_conβt3; /• λββign packing dββcriptor. •/ indβx++; nbytβe index-index; /*Savβ poβition of the "number of */ ~ /• byteβ thie row" for the let row. •/ layβr_lβngth « 22;

/•

*• The following run length encodeβ the data in the grid.

,,, _Λ _ 15504

- 64 -

•* Loop through all of the rowe in the grid with the following. */ for(row-MIN Y;row<number rowe;row++)

{ prev col—1; run-MTN VλLUE; nbyteβ-MIN_VλLUE; index-*-*-; cp-(char *)index;

/*

** Loop through all of the bine in the current row.

*/ for(bin-MIN X;bin<number bine;bin-*--*-)

{

/*

** Extract the color value from the grid at the

** current bin and row location. Mask off the upper

** bite βo that only the 4 bit color value remainβ.

*/ col * grid[bin] [row] ;

/*

*• Now that the color haβ been determined, either start

** a new "run" or update an exiβting "run" with the following.

*/ if( ((prev coll—1)££(coll-prev col)) || (run—15) )

{

/*

** Start a new "run" at the beginning of a row, when

*• the color changeβ, and when the current run iβ 15

** (which iβ the maximum run value). The run/color

*» combination ie aββigned to the buffer at the poβition

** of the character pointer (cp).

*/ byte-(run*16)+prβv_col;

•(cp)-byte; ^" cp++; nbytββ-*"*-; run-1; prev col - col;

} else

{

/* Otherwise, update an exieting run and go on to the next

** grid box bin value.

*/ run++; prev col - col; > >

/*

** When all bine in the current row have been proceeβed

** then check if a value remainβ to be etored.

*/ if (run>l) {

/*

** Store the laβt "run color" combination in the

6/15504

- 65 -

** output buffer with the following. */ byte-(run*16)+prev_co1; •(cp)-byte; ^" cp++; nbytββ-*"*-;

}

I I*.

•• If the number of byteβ iβ odd then pad •• one more byte with a value of zero.

*/ i (nbytββ % 2)

/*

** λββign the number of byteβ for thie row to the

*• buffer in the proper location. */

•(nbytββ index)-nbytββ; layer_lβngth+-nbytββ;

/*

** Update layβr_lβngth to include 2 byteβ

*• for the nbytee value.

*/ layβr_length+-2;

/*

•* λdvance the index pointer eo that it pointβ to the

** location of the next row'β byte length.

*/ index-(uneigned ehort *)cp; nbytee index-index; >

/*

** Now that all of the rowe in the grid have been proceeeed,

** convert the layer length value to 2 "halfwordβ"

** and aβeign thoβe to the layer eize poβition in the

** product eymbology layer one header.

*/ long2ehort.1 value ■ layβr_lβngth; mew_layβr_lβn - long2ehort.β_valuβ(0]; lβw~layβr~len - long2βhort.e~valuβ(lj; *πww_layer_lβn__indβx - mβw__layβr_len; *lβw_layer~lβn ^— index - lβw-layer__lβn; block_lβngth - layβr_lβngth+16;

/*

•• Convert the block length value to 2 "halfwordβ"

•* and aββign thoβe to the block length poβition in the

** product eymbology block header.

*/

long2βhort . l_value - block_length; mβw_block_len - long2 ehort 7β_value[0J ; lβw ^" block~len - long2ehort . β~valuβ( l j ; *m_rw_block_len_index - mβw_bTock len; *lβw~block_len_index - lβw_block~len;

/*

** Convert the meββage length value to 2 "halfwordβ"

•• and aaeign thoβe to the meββage length poβition in the

** product meββage header block.

*/ total_length - block length+18+102; meeβagβ_length - total length; long2βhort.l_value - mβββagβ_length; mew_mβg_len - long2βhort.β_value(0]; lew_mβg_len - long2βho_t.e_value[lj; •mβw_mβg_len index - mew meg len; *lβw_mβg~len~indβx - lβw~mβg~len;

/*

^• • Write the product to the file and cloβe the file.

*/ nwrite - fwrite(buf,total length,1,fpw); if(nwrite—NODλTλ)

{

(void)print ( "error writing product file, errno-%d\n" ,errno) ; return ( 1 ) ; > etatuβ - ffluβh(f w); if{ βtatuβ I- 0 ) {

(void)printf("error: error fluβhing product file, errno-%d\n",errno); return(1); > etatuβ * fcloβe(fpw); if( βtatuβ 1- 0 ) {

(void)printf("error: error cloβing product file, errno-%d\n",errno) ; return(l) ; > printf("Product iβ complete. \n\n");

(void)ffluβh(βtdout);

(void)ffluβh(βtderr); } return(0);

}

/*

•* Routine to obtain lateβt vereion of the requeeted product for the requeeted β

*/ int get_latββt(short eite id, short prod id, char *β)

{

/*

*• Thie paβββe a βtring containing the complete pathname

•• to the lateet product in the data baβe of the type epecified

•• by the product id and for the βite βpecified by eite id.

•*

•• A etatuβ code iβ returned by the routine. The following lietβ the *• return valuβe: return - 0 -> product name found (OX to proceed). •• return - 1 -> no product found in directory.

** return - 2 -> invalid product id number.

•• return - 3 -> invalid aite id number.

•• return - 4 -> popen failed.

- 67 -

char *teet; char command[512); char etring[512]; conβt char *commandl - "/uβr/bin/lβ -r "; conβt char *command2 « "777777.???? 2>/dev/null|/uer/ucb/head -1 "; conat char *βtring_homβ - "/home/Productβ/radar/";

FILE *cfd; ehort j;

/*

•* Copy the firet part of the command to the command βtring.

*/

(void)βtrcpy(command,command1) ;

/*

•* Concatenate the main (home) directory βtring to the

** command βtring.

*/

(void)etrcat(command,βtring_homβ);

/•

** Search for the βpecified eite id with the following loop.

*/ for(j-0; j<NUM SITE PλTH; j++)

{ if(βitβ_path(jj.βite id — βite id )

{

/*

** When the βite id iβ found, concatenate the βite

^♦ • name path βtring to the command βtring and then

** jump (break) out of the look.

*/

(void) etrcat (command, βite pathf j ] . βite_tβxt ) ; goto jumpl_hβre; ~ y

>

/*

** The following 3 statementβ are executed only when the βite id

** iβ not found. Null the path etring, copy it to the output βtring,

•* and return from the routine with a nonzero value indicating

•• that the product wae not found.

*/ βtring[0]-'\0';

( oid)βtrcpy(β,etring); return(2); jumpl_hβrβ:; /* jump to here and continue. */

/*

*• The following βearcheβ the liet of product id'β for the

** βpecifed id.

*/ for(j-0/ j<NUK PROD PλTH; j++) {

if(prodjpath[j].prod id -- prod_id ) (

/*

** Whan the product id iβ found, concatenate the product

** name path βtring to the command etring and then

** jump (break) out of the look.

*/

(void)etrcat(command,prod_pat [j].prod_text);

/*

*• The following ia executed only when the product id iβ

** not found. Null the output βtring, copy it to the output,

** and return from the routine with a nonzero value indicating

** an invalid path.

*/ stringf0]-'\0';

(void)etrcpy(e ,etring); return(3); jump2_herβj ; /* jump to here and continue. */

/*

** Concatenete the laet part of the command on to the end of the

•* command βtring.

*/

(void)etrcat(command,command2);

/*

*• Open a pipe and execute the command under the "βh" shell.

•• The reβult iβ returned in a βtring pointed to by the pointer "cfd".

*/ cfd-popen(command,"r"); if(cfd—NULL)

{

/*

•* Thie handleβ the unlikely caee that the popβn function failed.

** Print a meββage and return with a nonzero βtatuβ code.

*/

(void)printf("gβt_latββt error -> cfd iβ NULL: errno-%d\n",errno); return(4); ~ }

/*

*• Get the etring with the following. Remember that fgetβ returns

** a βtring which include a newline character at the end (which muet

** be removed by replacing it with a NUL character.

*/ etring[0]-' \0' ;

( oid) etrcpy (β, βtring) ; tββt-f gets ( βtring, βizeof ( βtring) , cfd) ; pcloβe(cfd) ;

/*

*• Check the firet character of the βtring for a null.

*/ if(βtring[0]—'\0')

- 69 -

/*

** Thie handleβ a βtring with the firet character a null.

** It meanβ that the requeeted product waβ not found.

** The null βtring iβ paββed and the routine returns with a

** nonzero value.

*/

(void)etrcpy(β,βtring); rβturn(l);

>

/*

** When the etring iβ valid, replace the laβt character (a carriage return)

** with the null character. Thie inβureβ that juβt the path name iβ

*• returned to the calling routine. βtring[βtrlen(etring)-1]-'\0';

/*

* ^• Copy the βtring to the output etring which iβ paββed to the calling

** routine. The value of zero iβ returned with return indicating

** that a valid product waβ found.

*/

(void)etrcpy(β,βtring); return(0); >

The apparatus and method of the present invention converts a radar data bin (or the like) to a grid location or box in a resulting mosaic, rather than vice versa. If the reverse is done, occasional radar data bins could be missed.

The present invention can construct a mosaic using any data that is geographical or otherwise territorial (in its broadest sense) in nature or relating to merging into one mosaic points from a plurality of surfaces. However, it is strongly preferred that such data be in a raster format. Data from different sites can thereby be combined to produce a single mosaic.

The present invention can be practiced using prebuilt lookup tables, or by performing the reprojection calculations in real time. However, use of lookup tables is preferred, for faster operation and quicker response in real time.

Some of the many advantages of the invention should now be readily apparent. For example, apparatus and method are provided for combining of data from multiple sources into one or more mosaics such as would be suitable for

utilization by end users. Such apparatus and method i capable of accurately and efficiently constructing suc mosaics. Also, the location of any holes in such a mosaic can thereby be identified. The present invention also provides for filling such holes.

Thus there has been provided apparatus and method for combining data from a plurality of sites into a single image covering a desired geographical area. Such apparatus and method are capable of combining, in real time, radar data (such as weather radar data) or other data in digital format from multiple radar sites into a mosaic covering a geographical area. Such apparatus and method are furthermore capable of combining data from a plurality of sources into a single image that can be run on a general purpose computer such as a personal computer. Such apparatus and method are capable of producing a mosaic of weather radar data that can provide location and intensity of precipitation and other meteorological phenomena throughout a large geographical area. When lookup tables are used, the present invention can accurately combine a substantial amount of data from a plurality of desired sites fast enough to provide a current depiction of the weather or other phenomena.

Obviously, many modifications and variations of the present invention are possible in light of the above teachings. It is therefore to be understood that the foregoing embodiments are presented by way of example only and that, within the scope of the appended claims and equivalents thereto, the invention may be practiced otherwise than as specifically described.