Optimization of Bathymetry Database for Coastal Areas

Problem statement: A database containing gridded bathymetries with dif ferent coverages and resolutions has been established. Approach: The interpolation techniques of the gridded bathymetries have been developed to define the sphe rical, rectangular and curvilinear coordinate grid systems. Results: Subsequently, merging of the bathymetry results h as also been developed to define the optimal bathymetry. Furthermore, the bathymetri es were established in order to apply for the Gulf of Thailand (GoT) and also for the storm wave and s urge models of the Coasts of Thailand (CoT). Conclusion: The bathymetry covering the coastal waters was prod uce to extend the applicability of the operational 2D and 3D ocean models. In addition , he data sets used in this study were integrated and collected in order to be the database and post pr cessing tools to convert the output to the speci fic formats required for various operational models at the CoT.


INTRODUCTION
A high quality fine resolution bathymetry data set has long been a priority for the operational ocean modelers. Currently ten bathymetry data sets covering different parts of the Coasts of Thailand (CoT) are available from GEODAS (Amante and Eakins, 2008;Edwards, 1989;Sloss, 2006). These data sets have different coverages, resolutions and accuracies in different areas. Ongoing the mathematical oceanic models (Blumberg and Mellor, 1987;Booij et al., 1999;Gunther et al., 1992;Hasselmann et al., 1988;Komen et al., 1994) with the bathymetry data (Amante and Eakins, 2008;Edwards, 1989;Sloss, 2006) have been carried out in many organizations and in the global, regional and coastal scales (Wittmann and Farrar, 1997). In recent years of the operational ocean prediction system, the data have faced a frequent demand of predictions of water levels (Vongvisessomjai et al., 2008), waves, surges (Wannawong et al., 2010d), tides (Wannawong et al., 2010b), currents and circulations at a very high resolution in many locations (Wannawong et al., 2010c). This leads to increase the high resolution of the regional and coastal models continuously. A quick setup of these models in new regions is a great asset in the operational ocean prediction systems such as navigation, search, rescue and oil spill combating. An important premise of the new model setup is the presence of the bathymetry in the required resolution. The database containing data sets covering the CoT from 99-111°E in longitude and from 2-14°N in latitude has therefore been established based on existing bathymetry data at the CoT (Fig. 1). A software package has been developed to interpolate the data sets in order to define the spherical and rectangular coordinate grid systems (Wannawong et al., 2008;2010a). The software has also been applied to merge the results of the interpolation data sets into the best single predefined bathymetry.

MATERIALS AND METHODS
Data sets: Ten data sets were assembled and integrated into the database. These were listed as shown in Table 1 and Fig. 2. It is anticipated that more data sets will be incorporated into the database in the future. The data sets were stored in either spherical or Universal Transverse Mercator (UTM) coordinates (Hager et al., 1989). The ETOPO5 is a global topography which has been extracted for the entire areas of interests (Edwards, 1989). The original version of the ETOPO5 (1993) on a 5 min latitude/longitude grid (1 min of latitude = 1 nautical mile, or 1.853 km) was updated in June 2005 for the deep water. While the ETOPO5 data which are based on a satellite altimetry showed poor quality with shallow water.
Recently, other global bathymetries: ETOPO2 and ETOPO1 became available with higher resolutions (Amante and Eakins, 2008;Sloss, 2006;Smith and Sandwell, 1997). These three bathymetries have been denoted by the WAve Model Cycle 4 (WAMC4) (Gunther et al., 1992;Hasselmann et al., 1988;Komen et al., 1994;Monbaliu et al., 2000) and the Simulating Waves Nearshore (SWAN) models (Booij et al., 1999). The ETOPO2 original version (1997) on a 2 min latitude/longitude grid was updated in December 2, 2006 for the shallow water. Several studies were reported that the old version of ETOPO5 has been used to operate wave models from the eye storm generated at the Pacific Ocean entering into the South China Sea (SCS) and the Gulf of Thailand (GoT) (Vongvisessomjai, 2007;2009) and then the new version of ETOPO5 has been improved and used instead of the old version by Wannawong et al. (2010a;2010c). The ETOPO1 has also been utilized to operate the wave models at the CoT (Amante and Eakins, 2008;Wannawong, et al., 2010b;2010d;2010e). The latest version of the ETOPO1 (on July 28, 2008) (Amante and Eakins, 2008) on a 1 min latitude/longitude grid with the high resolution was available on Ice Surface (top of Antarctic and Greenland ice sheets) and Bedrock (base of ice sheets) versions. The ETOPO1 (Bedrock version) was chosen to apply inside domains (Wannawong, et al., 2010b;2010d;2010e). The details of the combined bathymetries and nested grid domains to drive the Princeton Ocean Model (POM) (Monbaliu et al., 2000) were described by Xia et al. (2008).
The bathymetry used in the operational Thai ocean general circulation model has also been incorporated into the database. The CoT was covered by the curvilinear and rectangular coordinates with the coarse resolutions in a range of 2-55 km, by the horizontal coordinates with 19.3×9.2 km 2 and by the vertical coordinate with 10 sigma levels (Wannawong et al., 2008). While other Thai nearshore areas were covered by the fine grid resolution of 0.1°×0.1° of the horizontal coordinates and by 21 sigma levels of the vertical coordinate (Wannawong et al., 2010b;2010c;2010d). The Dynamics of Connected Seas (DCS) provided the high resolution bathymetry for the larger area covering the CoT. The wave and storm surge models at the CoT have been applied by the DCS bathymetry between the GoT and CoT (Wannawong et al., 2010a;2010b;2010c;2010d;2010e). At the GoT and CoT, the bathymetries have been successfully improved. The bathymetry data of the GoT was interpolated to obtain the high resolution at the CoT. The results of bathymetry were denoted by the POM and MIKE21 models. The quality of POM, MIKE21 and SWAN bathymetry models were clearly higher than the DCS bathymetry.
Format file: Five data formats were defined as shown in Table 2. The first column showed the data sets used in this study. The second and third columns presented the formats of the longitude and latitude of the data files. The fourth and fifth columns expressed the sizes of the grid cell in the form files. The depths of sea water were written in the sixth column in meters. The negative values of the depths denoted the grid cell points below sea levels while the positive values denoted the land grid cells. The height of topography might not exceed 10 m which was defined as a true land. The last column was a flag column which was presented not to use and can be left empty. All data sets shown in Table 2 were merged in order to apply in the UTM-47 coordinates. The modification of the bathymetry data in the UTM-47 coordinates is shown in Table 3.

Parameter file:
In principle, all necessary information about the data files can be obtained from the data files themselves. In practice, however, it is very time consuming to extract an amount of this information and it is therefore included in a parameter file. In the parameter file, some additional information needed for the interpolation of new bathymetries is also included. A simple example of the parameter file is displayed below:  (8)  999 format (i6)  888 format(49i6) open(9,file='bathymetry_worachat.asc',form='formatte d') write(9,888)((deep(i,j),i=1,nx),j=1,ny) close (9)  The example showed that these reduced bathymetry database consisted of three data sets named as ETOPO5, bsh_worachat and DCS. The parameter values of each data set were separated by an empty line. The first line of each data set contained the name of the input data file (The extension was used to identify which coordinate system is used in the data file). The information of the data set was written in the second line while the following lines contained the information determined where to use the bathymetry in the interpolation. The first two values of the second line defined the coordinates of the origin of the data set (lower-left corner) while the third and fourth values held the grid size. Depending on the coordinate systems used in the data sets either degrees (spherical coordinates) or meters (UTM-47 coordinates) were used as units for the data fields. The fifth and sixth values contained the number of data points in each of the coordinate directions. The last field was an identifier for each data set. This value was used in the flag column of the interpolation data set and should be set between 1 and 99 (0 was used for the merged data set). The grid of the data set was calculated based on these parameters. It was therefore necessary to specify the grid sizes with a high accuracy as a round-off error in one of these parameters which sum up over the grid. The magnitude of the maximum error in a coordinate direction can be calculated by multiplying the number of grid points in the direction by the error of a single grid cell. In the above example, the data sets specified the longitudinal grid size of the ETOPO5, ETOPO2 and ETOPO1 as 0.0833333, 0.03333333 and 0.01666667 respectively with similar format (i6). For the sample of ETOPO5 data, the error of the single grid cell was 10-8. The maximum error then became the order of 721×10-8 ~ 10-5° longitude. The parameter file therefore also contained the information to determine where and which data set is used in the merged bathymetry. For each data set in the parameter file zero or more bounding boxes were specified (in spherical coordinates). Inside each box of the interpolation data, the specified data set was used for the merged data set. The western limit of the bounding box was firstly specified, followed by the southern, eastern and northern limits (see the above example).
Bathymetry creation: For a new bathymetry generation, the following steps have to be made: • Choosing existing data sets which will be used to produce the new data set, • Interpolation and merging the selected bathymetry data sets onto the required area, resolution and quality • Checking the quality of the new generated bathymetry Data set selections for the interpolation: Since the existing bathymetry data have different coverages (availability), resolutions and qualities, it is recommended to specify the first input data set as a base map covering the entire area of interests and let the following data sets to be nested inside this area based on the quality and resolution. It is important to check whether the data set covers the entire bounding box. If this is not the case, spurious output might occur. The data set that is specified lastly in the parameter file takes precedence when bounding boxes overlap. If the base map is specified as the first data set in the parameter file, all other data sets would thus be used although they always overlap the base map.

RESULTS
An area weighting interpolation algorithm is used. The depths of the output grid cells were calculated by averaging the depths of the input cells weighed by the fraction of the output cell area that is covered by each input cell. This process is applied both when the output grid is finer than the input grid and when the opposite is true (Fig. 3). Land areas are included in the calculations using a maximum value of +10 m. If the land (in the input grid cells) is more covered than a half of the output grid cell area, the depth is set to +10 m. For the spherical to spherical interpolations, the convergence of the meridians is taken into an account when calculating areas. The only degradation of data therefore stems from the implicit smoothing in the interpolation process. It is more complicated to interpolation data from UTM-47 coordinates to the spherical coordinates. As illustrated in Fig. 4 the UTM-47 input grid is generally rotated compared to the spherical output grid. The angle of rotation is a function of both longitude and latitude. In the interpolation coordinates of the corners of the spherical and rectangular coordinates, output grid cell are transformed to UTM-47 coordinates. Neglecting the curvature of a single output grid cell, these corner coordinates define a rectangular box in the UTM-47 grid. Each input grid cell depth is then weighed by its area inside the rectangular box.

DISCUSSION
The database containing gridded bathymetries of different coverages and resolutions has been established. A tool for the gridded bathymetry interpolations to a user defined the spherical, rectangular and curvilinear coordinate grid systems and subsequently merging the bathymetry results to a user defined the optimal bathymetry has been developed. The example of the interpolation was presented in this report (Appendix A). Furthermore, the bathymetry has been produced for the GoT by the storm wave and surge models. For the CoT, the bathymetry covering the coastal waters was produced to extend the applicability of the operational 2D and 3D ocean models. In the future, more data sets will be integrated into the database and post processing tools for converting the output to specific formats required for various operational models at the CoT will be developed. An extension of the interpolation algorithm is also able to interpolate the UTM coordinates which might be useful in the future. The data sets with higher quality and resolution will be assembled and integrated into thedatabase as they become available. Presently, the highest resolution of the data set was 0.33 nm in very limited areas. There was a clear need for higher resolution in large areas. More measurements are needed in certain areas as the offshore area and the shelf break. It will be necessary to share the bathymetry data between institutes and to improve the data exchange networks in order to produce these high resolution large area bathymetries.

CONCLUSION
It is self-evident that the best result is obtained by using the high quality of input bathymetries and with the finest possible resolution. The interpolation algorithms produce the acceptable results as long as the input bathymetry is comparable to the output bathymetry or the resolution of output bathymetry is finer than that of the input bathymetry. If the output bathymetry is much finer than the input bathymetry, staircase coastlines and isobaths will probably emerge. The problem can be understood by considering on the left hand side of Fig. 3. It is clearly considered that if the output grid cell is much smaller than the input grid cell then the depths of most of the output grid cells will just be the value of the input grid cell inside which they are fully contained. Only at the borders the depths will be calculated using more than one input grid cell. The problem could easily be elevated by using an algorithm using the squared inverse distance to the nearest four points (disregarding the actual sizes of the grid boxes). This algorithm would be poor for fine to coarse grid interpolations and the choice of the interpolation algorithm would thus have to depend on the resolutions of the grids. For the interpolations from the UTM-47 to spherical and rectangular coordinates, this would result in a very severe complication as a single interpolation could contain both coarse to fine and fine to coarse interpolations (due to the convergence of meridians). It has therefore been decided to keep the single interpolation algorithm. If the algorithm is used for this type of interpolation, smoothing the bathymetry results might reduce the problem.

Merging bathymetries:
The merged data set was produced by specifying where each data set is used. The depths in the nearest grid points to the division lines between two data sets were calculated as the average of these data sets. This smoothed the transition between the data sets somewhat. It is possible that more advanced merging/smoothing algorithm will be implemented in the future.

Appendix A:
Interpolation code of bathymetry: In this appendix, the simple example illustrating how to interpolate a bathymetry was reported. The bathymetry in the spherical coordinates was with the resolution of 2×2 min or 1×1 min (1 min of latitude = 1 nautical mile, or 1.853 km) for water surrounding the GoT. More precisely, the target domain was focused on the bathymetry which covered the area from 99-101° east and 12-14° north (Fig. 6). The most appropriate input bathymetry for this resolution area was the DCS bathymetry with a resolution of 1 nm ( Fig. 5 and 6 and accompanying with Table 3). All data sets specified in the parameter file were then interpolated to the grid and merged to produce the best bathymetry (Fig. 7a-c). A new interpolation bathymetry is created by specifying the coordinates of the desired grid as arguments for the interpolation program. More specifically, the (south-western) origin of the grid (lon. = 99°E, lat. = 12°N), resolution (dlon. = (1/60)°, dlat. = (1/60)°) for the ETOPO1 data and the number of points in each direction (nlon. = 49, nlat. = 49) must be given as command line arguments.