[4297] | 1 | """Script for running tsunami inundation scenario for boxing day Phuket, |
---|
| 2 | Thailand. Source data such as elevation and boundary data is assumed to be |
---|
| 3 | available in directories specified by project.py The output sww file is stored |
---|
| 4 | in project.output_time_dir |
---|
| 5 | |
---|
| 6 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
| 7 | the elevation data and a simulated submarine landslide. |
---|
| 8 | |
---|
| 9 | John Jakeman - ANU 2007 |
---|
| 10 | """ |
---|
| 11 | |
---|
| 12 | #use _read_asc() in data_manager |
---|
| 13 | |
---|
| 14 | |
---|
| 15 | #------------------------------------------------------------------------------ |
---|
| 16 | # Import necessary modules |
---|
| 17 | #------------------------------------------------------------------------------ |
---|
| 18 | |
---|
| 19 | # Standard modules |
---|
| 20 | from os import sep |
---|
| 21 | from os.path import dirname, basename |
---|
| 22 | from os import mkdir, access, F_OK |
---|
| 23 | from shutil import copy |
---|
| 24 | import time |
---|
| 25 | import sys |
---|
| 26 | from time import localtime, strftime |
---|
| 27 | |
---|
| 28 | # Related major packages |
---|
| 29 | from anuga.shallow_water import Domain |
---|
| 30 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts, ferret2sww |
---|
| 31 | from anuga.geospatial_data.geospatial_data import * |
---|
| 32 | from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files |
---|
| 33 | from anuga.utilities.polygon import polygon_area |
---|
| 34 | |
---|
| 35 | # Application specific imports |
---|
| 36 | import project # Definition of file names and polygons |
---|
| 37 | |
---|
| 38 | def convert_arcgis_latlon_list_to_utm(points): |
---|
| 39 | #Used because arc gis produced csv files put lat lon in |
---|
| 40 | #reverse order to those accpeted by convert_latlon_to_utm() |
---|
| 41 | |
---|
| 42 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
| 43 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
| 44 | old_geo = Geo_reference() |
---|
| 45 | utm_points = [] |
---|
| 46 | for point in points: |
---|
| 47 | zone, easting, northing = redfearn(float(point[1]),float(point[0])) |
---|
| 48 | new_geo = Geo_reference(zone) |
---|
| 49 | old_geo.reconcile_zones(new_geo) |
---|
| 50 | utm_points.append([easting,northing]) |
---|
| 51 | |
---|
| 52 | return utm_points, old_geo.get_zone() |
---|
| 53 | """ |
---|
| 54 | #------------------------------------------------------------------------------ |
---|
| 55 | # Preparation of topographic data |
---|
| 56 | # |
---|
| 57 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 58 | # Do for coarse and fine data |
---|
| 59 | # Fine pts file to be clipped to area of interest |
---|
| 60 | #------------------------------------------------------------------------------ |
---|
| 61 | print"project.bounding_polygon_latlon",project.bounding_polygon_latlon |
---|
| 62 | |
---|
| 63 | # topography and bathymetry directory filenames |
---|
| 64 | bathymetry_file = project.bathymetry |
---|
| 65 | |
---|
| 66 | #gets time for new dir |
---|
| 67 | time = strftime('%Y%m%d_%H%M%S',localtime()) |
---|
| 68 | |
---|
| 69 | # creates DEM from asc data |
---|
| 70 | print "creates DEMs from asc data" |
---|
| 71 | convert_dem_from_ascii2netcdf(bathymetry_file, basename_out=project.bathymetry_points, use_cache=True, verbose=True) |
---|
| 72 | |
---|
| 73 | #creates pts file for onshore DEM |
---|
| 74 | print "creates pts file for onshore DEM" |
---|
| 75 | dem2pts(project.bathymetry_points, |
---|
| 76 | use_cache=True, |
---|
| 77 | verbose=True) |
---|
| 78 | |
---|
| 79 | print'create Geospatial data objects from topographies' |
---|
| 80 | G = Geospatial_data(file_name = project.bathymetry_points + '.pts') |
---|
| 81 | |
---|
| 82 | bounding_polygon, zone = convert_arcgis_latlon_list_to_utm(project.bounding_polygon_latlon) |
---|
| 83 | print bounding_polygon |
---|
| 84 | print 'zone of bounding polygon', zone |
---|
| 85 | print 'Area of bounding polygon', polygon_area(bounding_polygon)/1000000.0 |
---|
| 86 | |
---|
| 87 | print'clip combined geospatial object by bounding polygon' |
---|
| 88 | G_clipped = G.clip(bounding_polygon) |
---|
| 89 | |
---|
| 90 | print'export combined DEM file' |
---|
| 91 | G_clipped.export_points_file(project.combined_dir_name + '.pts') |
---|
| 92 | """ |
---|
| 93 | |
---|