| 124 | A user may have access to a range of different resolution DEMs and raw data points (such |
| 125 | as beach profiles, spot heights, single or multi-beam data etc) and will need |
| 126 | to combine them to create an overall elevation data set. |
| 127 | |
| 128 | If there are multiple DEMs, say of 10m and 25m resolution, then the technique is similar to |
| 129 | that defined in the Cairns example described earlier, that is: |
| 130 | |
| 131 | {{{ |
| 132 | convert_dem_from_ascii2netcdf(10m_dem_name, use_cache=True, verbose=True) |
| 133 | convert_dem_from_ascii2netcdf(25m_dem_name, use_cache=True, verbose=True) |
| 134 | }}} |
| 135 | followed by: |
| 136 | {{{ |
| 137 | dem2pts(10m_dem_name, use_cache=True, verbose=True) |
| 138 | dem2pts(25m_dem_name, use_cache=True, verbose=True) |
| 139 | }}} |
| 140 | These data sets can now be combined by |
| 141 | {{{ |
| 142 | from anuga.geospatial_data.geospatial_data import * |
| 143 | G1 = Geospatial_data(file_name = 10m_dem_name + '.pts') |
| 144 | G2 = Geospatial_data(file_name = 25m_dem_name + '.pts') |
| 145 | G = G1 + G2 |
| 146 | G.export_points_file(combined_dem_name + �.pts�) |
| 147 | }}} |
| 148 | This is the basic way of combining data sets, however, the user will need to |
| 149 | assess the boundaries of each data set and whether they overlap. For example, consider |
| 150 | if the 10m DEM is describing by {{{poly1}}} and the 25m DEM is described by {{{poly2}}} |
| 151 | with {{{poly1}}} completely enclosed in {{{poly2}}} as shown here: |
| 152 | \begin{figure}[hbt] |
| 153 | \centerline{\includegraphics{graphics/polyanddata.jpg}} |
| 154 | \caption{Polygons describing the extent of the 10m and 25m DEM.} |
| 155 | \label{fig:polydata} |
| 156 | \end{figure} |
| 157 | To combine the data sets, the geospatial addition is updated to |
| 158 | {{{ |
| 159 | G = G1 + G2.clip_outside(Geospatial_data(poly1)) |
| 160 | }}} |
| 161 | For this example, we assume that {{{poly2}}} is the domain, otherwise an additional dataset |
| 162 | would be required for the remainder of the domain. |
| 163 | |
| 164 | This technique can be expanded to handle point data sets as well. In the case |
| 165 | of a bathymetry data set available in text format in an {{{.csv}}} file, then |
| 166 | the geospatial addition is updated to |
| 167 | {{{ |
| 168 | G3 = Geospatial_data(file_name = bathy_data_name + '.csv') |
| 169 | G = G1 + G2.clip_outside(Geospatial_data(poly1)) + G3 |
| 170 | }}} |
| 171 | The {{{.csv}}} file has the data stored as {{{x,y,elevation}} with the text {{{elevation}}} |
| 172 | on the first line. |
| 173 | |
| 174 | The coastline could be included |
| 175 | as part of the clipping polygon to separate the offshore and onshore datasets if required. |
| 176 | Assume that {{{poly1}}} crosses the coastline |
| 177 | In this case, two new polygons could be created out of {{{poly1}}} which uses the coastline |
| 178 | as the divider. As shown in Figure \ref{fig:polycoast}, {{{poly3}}} describes the |
| 179 | onshore data and {{{poly4}}} describes the offshore data. |
| 180 | \begin{figure}[hbt] |
| 181 | \centerline{\includegraphics{graphics/polyanddata2.jpg}} |
| 182 | \caption{Inclusion of new polygons separating the 10m DEM area into an |
| 183 | onshore (poly3) and offshore (poly4) data set.} |
| 184 | \label{fig:polycoast} |
| 185 | \end{figure} |
| 186 | Let's include the bathymetry |
| 187 | data described above, so to combine the datasets in this case, |
| 188 | {{{ |
| 189 | G = G1.clip(Geospatial_data(poly3)) + G2.clip_outside(Geospatial_data(poly1)) + G3 |
| 190 | }}} |
| 191 | |
| 192 | Finally, to fit the elevation data to the mesh, the script is adjusted in this way |
| 193 | {{{ |
| 194 | domain.set_quantity('elevation', |
| 195 | filename = combined_dem_name + '.pts', |
| 196 | use_cache = True, |
| 197 | verbose = True) |
| 198 | }}} |