Changeset 5012 for anuga_work/production

Timestamp:

Feb 8, 2008, 5:13:08 PM (17 years ago)

Author:

nick

Message:

update geraldton

Location:

anuga_work/production/geraldton

Files:

• build_geraldton.py (modified) (6 diffs)
• project.py (modified) (8 diffs)

anuga_work/production/geraldton/build_geraldton.py

@@ -25 +25 @@
 # Related major packages
 from anuga.shallow_water import Domain
-from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts,start_screen_catcher, copy_code_files
 from anuga.geospatial_data.geospatial_data import *
-from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
 
 # Application specific imports
@@ -37 +36 @@
 #------------------------------------------------------------------------------
 
-start_screen_catcher(project.output_build_time_dir)
 
-print 'time stamp: ',project.gtime
+print 'time stamp: ',project.build_time
 print 'USER: ', project.user
 
@@ -46 +44 @@
                dirname(project.__file__)+sep+ project.__name__+'.py' )
 
+start_screen_catcher(project.output_build_time_dir)
 
 
@@ -59 +58 @@
 onshore_in_dir_name = project.onshore_in_dir_name
 coast_in_dir_name = project.coast_in_dir_name
-#island_in_dir_name = project.island_in_dir_name
+island_in_dir_name = project.island_in_dir_name
 offshore_in_dir_name = project.offshore_in_dir_name
 
 onshore_dir_name = project.onshore_dir_name
 coast_dir_name = project.coast_dir_name
-#island_dir_name = project.island_dir_name
+island_dir_name = project.island_dir_name
 offshore_dir_name = project.offshore_dir_name
 
@@ -70 +69 @@
 print "creates DEMs from ascii data"
 convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
-#convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
 
 #creates pts file for onshore DEM
 print "creates pts file for onshore DEM"
 dem2pts(onshore_dir_name,
-#        easting_min=project.eastingmin,
-#        easting_max=project.eastingmax,
-#        northing_min=project.northingmin,
-#        northing_max= project.northingmax,
         use_cache=True, 
         verbose=True)
 
 #creates pts file for island DEM
-#dem2pts(island_dir_name, use_cache=True, verbose=True)
+dem2pts(island_dir_name, use_cache=True, verbose=True)
 
 print'create Geospatial data1 objects from topographies'
 G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
-G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
-#G3 = Geospatial_data(file_name = island_dir_name + '.pts')
-G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
+print'finished reading in %s.pts' %onshore_dir_name
+G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt')
+G3 = Geospatial_data(file_name = island_dir_name + '.pts')
+G4 = Geospatial_data(file_name = offshore_in_dir_name + '.txt')
+print'finished reading in files'
+
+G_onshore_clip = G1.clip(project.poly_all, verbose=True)
+print 'finished clip'
+G_onshore_clip_clipout=G_onshore_clip.clip_outside(project.poly_cbd, verbose=True)
+print 'finished clipout'
+# reduce resolution by 9 time (1/9) of the original file. This was determined
+# as the resolution for this region is going to be 20000 and the asc grid cellsize is
+# 30m there (900). The sqrt of 20000 is about 140 or 150, a 150m x150m grid would contain
+# 36 points at a 30m grid spacing and we only need 4. So 4/36 is 1/9 or 0.11111
+G_onshore_clip_clipout_reduced=G_onshore_clip_clipout.split(0.111, verbose=True)
+
+G_cbd = G1.clip(projec.poly_cbd,verbose=True)
+G_cbd_reduced = G_cbd.split(0.25, verbose=True)
+
+G_bathy = G4.clip(project.poly_all)
+G_island = G3.clip(project.poly_all)
+G_coast = G2.clip(project.poly_all)
 
 print'add all geospatial objects'
-G = G1 + G2 + G_off 
+G = G_onshore_clip_clipout_reduced + G_cbd_reduced + G_bathy + G_island + G_coast 
 
 print'clip combined geospatial object by bounding polygon'
@@ -102 +116 @@
 if access(project.topographies_dir,F_OK) == 0:
     mkdir (project.topographies_dir)
-G.export_points_file(project.combined_dir_name + '.xya')
+G.export_points_file(project.combined_dir_name + '.txt')
 #G_clipped.export_points_file(project.combined_dir_name + '.xya')
+G_small, = G.split(0.1, verbose)
+G_small.export_points_file(project.combined_dir_name_small + '.txt')
 
 '''

anuga_work/production/geraldton/project.py

@@ -31 +31 @@
 #Making assumptions about the location of scenario data
 state = 'western_australia'
-scenario_name = 'gerladton'
-scenario = 'busselton_tsunami_scenario_2006'
+scenario_name = 'geraldton'
+scenario = 'geraldton_tsunami_scenario'
 
 #Maybe will try to make project a class to allow these parameters to be passed in.
@@ -41 +41 @@
 finaltime=25000
 export_cellsize=50
-setup='final'
+setup='trial'
 source='other'
 
@@ -65 +65 @@
 
 # onshore data provided by WA DLI
-onshore_name = 'DLI_orthophoto_DEM' # original
-#islands
+onshore_name = 'landgate_dem_clip' # original
+island_name = 'dted_srtm_islands'
 
 # AHO + DPI data
-coast_name = '100k_coastline'
-offshore_name = 'Bathymetry'
+coast_name = 'coastline'
+offshore_name = 'geraldton_bathy'
 
 #final topo name
-combined_name ='busselton_combined_elevation.pts'
-combined_smaller_name = 'busselton_combined_elevation_smaller'
+combined_name ='geraldton_combined_elevation.pts'
+combined_smaller_name = 'geraldton_combined_elevation_smaller'
 
 anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
@@ -83 +83 @@
 #input topo file location
 onshore_in_dir_name = topographies_in_dir + onshore_name
-#island_in_dir_name = topographies_in_dir + island_name
+island_in_dir_name = topographies_in_dir + island_name
 
 coast_in_dir_name = topographies_in_dir + coast_name
@@ -89 +89 @@
 
 onshore_dir_name = topographies_dir + onshore_name
-#island_dir_name = topographies_dir + island_name
+island_dir_name = topographies_dir + island_name
 coast_dir_name = topographies_dir + coast_name
 offshore_dir_name = topographies_dir + offshore_name
@@ -128 +128 @@
 output_dir = anuga_dir+'outputs'+sep
 output_build_time_dir = output_dir+build_time+'_'+source+sep
-#output_run_time_dir = output_dir +run_time+dir_comment+sep
-output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
+output_run_time_dir = output_dir +run_time+dir_comment+sep
+#output_run_time_dir = output_dir+sep+run_time+sep
 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
 
@@ -150 +150 @@
 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
 
-# bounding polygon for study area
-poly_all = read_polygon(polygons_dir+'poly_all.csv')
-res_poly_all = 100000*res_factor
 
 ###################################################################
@@ -169 +166 @@
 
 #digitized polygons
-poly_busselton1 = read_polygon(polygons_dir+'neg20_pos10_polygon.csv')
-res_busselton1 = 10000*res_factor
-poly_busselton2 = read_polygon(polygons_dir+'neg5_pos5_poly_.csv')
-res_busselton2 = 500*res_factor
+# bounding polygon for study area
+poly_all = read_polygon(polygons_dir+'extent_new_points.csv')
+res_poly_all = 100000*res_factor
+
+poly_neg20_pos20 = read_polygon(polygons_dir+'neg20_pos20_points.csv')
+res_neg20_pos20 = 25000*res_factor
+
+poly_neg10_pos10= read_polygon(polygons_dir+'neg10_pos10_points.csv')
+res_neg10_pos10 = 5000*res_factor
+
+poly_cbd = read_polygon(polygons_dir+'cbd_points.csv')
+res_cbd = 500*res_factor
+
+poly_wallabi = read_polygon(polygons_dir+'island_wallabi_poly_points.csv')
+res_wallabi = 25000*res_factor
+poly_dingiville = read_polygon(polygons_dir+'island_dingiville_poly_points.csv')
+res_dingiville = 25000*res_factor
+poly_pelsaert = read_polygon(polygons_dir+'island_pelsaert_poly_points.csv')
+res_pelsaert = 25000*res_factor
+
 
 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],figname='boundingpoly2',verbose=False)
 
-interior_regions = [[poly_busselton1,res_busselton1],[poly_busselton2,res_busselton2]]
-
-boundary_tags={'back': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 14],
-               'side': [10,13], 'ocean': [11, 12]}
+interior_regions = [[poly_neg20_pos20,res_neg20_pos20],[poly_neg10_pos10,res_neg10_pos10],
+                    [poly_cbd,res_cbd],[poly_wallabi,res_wallabi],
+                    [poly_dingiville,res_dingiville],[poly_pelsaert, res_pelsaert]]
+
+boundary_tags={'back': [2, 3, 4, 5],
+               'side': [1, 6], 'ocean': [0, 7, 8, 9, 10]}
 
 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
 
-poly_mainland=read_polygon(polygons_dir+'initial_condition_south.csv')
+poly_mainland=read_polygon(polygons_dir+'initial_condition.csv')
 
 print 'min number triangles', trigs_min