Changes between Version 8 and Version 9 of ModellingQuestions

Oct 20, 2008, 1:54:45 PM (16 years ago)



  • ModellingQuestions

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