Changeset 5465


Ignore:
Timestamp:
Jul 4, 2008, 1:56:15 PM (17 years ago)
Author:
ole
Message:

updated boxingday scenario to use ordered sts file

Location:
anuga_work/development/boxingday08
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/boxingday08/boxingday_boundary_order.txt

    r5464 r5465  
    1 index,longitude,latitude
    2 90,472194.058708,193430.03714
    3 91,468910.350862,417814.660921
    4 0,336294.174659,659735.865848
    5 1,335940.023759,664344.009834
    6 2,335585.169534,668952.072994
    7 3,333458.007771,671787.84005
    8 4,333103.527881,676395.93421
    9 5,335584.913055,679586.176909
    10 6,338066.17797,683130.901664
    11 7,337712.337412,687384.526034
    12 8,340193.365548,691283.679243
    13 9,342674.942465,695182.864565
    14 10,343028.996544,698727.563095
    15 11,343383.181996,702626.692455
    16 12,345510.456522,706171.398037
    17 13,347636.763496,710070.561254
    18 14,348345.993026,713969.747598
    19 15,347991.999887,717868.90241
    20 16,348700.260229,721413.598337
    21 17,348700.74816,725312.746857
    22 18,348346.261066,729566.410886
    23 19,347637.463051,733465.539648
    24 20,346219.280415,737010.261839
    25 21,345864.407267,741263.85975
    26 22,345510.331768,745163.03914
    27 23,342674.316521,747289.866983
    28 24,340547.547708,750480.065429
    29 25,339484.479114,755088.19808
    30 26,341965.31107,757923.969427
    31 27,343028.978631,761468.662515
    32 28,338420.592193,761823.097012
    33 29,335585.22943,764304.383975
    34 30,334167.626591,768203.581251
    35 31,335584.735827,771393.786837
    36 32,337003.09138,775292.952764
    37 33,335939.319181,779546.605803
    38 34,333103.507345,782027.857453
    39 35,329204.500134,783800.237744
    40 36,327432.550064,787699.371861
    41 37,328141.129947,790889.621245
    42 38,329204.642377,795143.228566
    43 39,330268.050113,799396.879865
    44 40,331331.840794,803296.055516
    45 41,332394.68402,807195.228032
    46 42,332749.712168,811803.336354
    47 43,332749.085392,815347.988474
    48 44,331686.000498,819247.16197
    49 45,330977.236649,823500.800482
    50 46,330622.34502,827754.458718
    51 47,330268.46938,831299.140929
    52 48,329913.427618,835198.318404
    53 49,328141.18488,838742.994622
    54 50,325305.544465,841933.245518
    55 51,323533.487091,845123.425076
    56 52,321051.760662,848668.14935
    57 53,321051.898764,852921.802336
    58 54,320343.173374,857529.873001
    59 55,316443.694676,860011.188558
    60 56,312898.80714,861429.03829
    61 57,310772.31805,864973.73467
    62 58,309354.462261,868872.930518
    63 59,305454.953286,870999.733764
    64 60,301555.903121,872772.053637
    65 61,299429.071562,876671.238128
    66 62,299075.077259,880924.887741
    67 63,298365.530409,884824.028562
    68 64,297302.823004,889077.699828
    69 65,294820.924212,892267.85487
    70 66,290922.375994,894394.741886
    71 67,287377.363778,896167.086731
    72 68,284895.736516,899711.774988
    73 69,283123.376077,902901.93998
    74 70,280997.234969,905737.733813
    75 71,278515.811642,908927.944896
    76 72,277097.388615,911763.685916
    77 73,276389.008202,915308.431193
    78 74,277452.296071,919562.067608
    79 75,277097.344968,923461.192048
    80 76,274616.125999,927360.398292
    81 77,272489.492172,930196.118259
    82 78,270008.294828,932677.417175
    83 79,266817.763524,935867.66411
    84 80,262919.01753,937639.955008
    85 81,258665.020804,936576.567479
    86 82,254411.760125,936576.586726
    87 83,257956.617175,988187.298388
    88 85,431895.346948,1020339.55099
    89 84,456154.607059,1099910.74796
    90 87,485676.850637,892038.719914
    91 86,436057.013086,865548.249097
    92 93,422179.921754,873660.847571
    93 92,417478.780955,873745.865662
    94 88,533099.268408,829034.042352
    95 89,571839.718716,740632.637627
     1index, longitude, latitude
     290,98.75,1.75
     391,98.7200012207,3.77999997139
     40,97.5210189819,5.96663379669
     51,97.5177078247,6.00829792023
     62,97.5143890381,6.04996109009
     73,97.4951019287,6.07555246353
     84,97.4917831421,6.11721515656
     95,97.5141220093,6.14612770081
     106,97.5364532471,6.17824554443
     117,97.5331497192,6.21670341492
     128,97.5554733276,6.25202655792
     139,97.5778045654,6.28734970093
     1410,97.5809173584,6.31941461563
     1511,97.584022522,6.35468482971
     1612,97.6031646729,6.38679361343
     1713,97.6222915649,6.42210769653
     1814,97.6286087036,6.45738744736
     1915,97.6253128052,6.49264097214
     2016,97.6316299438,6.52471494675
     2117,97.6315383911,6.5599770546
     2218,97.6282272339,6.59843635559
     2319,97.6217193604,6.63368034363
     2420,97.6088027954,6.66570091248
     2521,97.6054840088,6.70415878296
     2622,97.602180481,6.73941135406
     2723,97.5764694214,6.75857067108
     2824,97.557144165,6.78736352921
     2925,97.5474014282,6.82900667191
     3026,97.569770813,6.85471820831
     3127,97.5792999268,6.8868021965
     3228,97.5375900269,6.88988161087
     3329,97.5118637085,6.91224050522
     3430,97.4989242554,6.94746017456
     3531,97.5116577148,6.97634935379
     3632,97.5243835449,7.01164960861
     3733,97.5146331787,7.05008459091
     3834,97.4888916016,7.07243967056
     3935,97.4535446167,7.08835077286
     4036,97.4373855591,7.12355518341
     4137,97.4437026978,7.15242481232
     4238,97.4532012939,7.19092082977
     4339,97.4626998901,7.22941732407
     4440,97.4722137451,7.26470851898
     4541,97.4817199707,7.2999997139
     4642,97.4847946167,7.3416800499
     4743,97.4846801758,7.37373304367
     4844,97.4749298096,7.4089589119
     4945,97.4683761597,7.44740056992
     5046,97.4650268555,7.48585319519
     5147,97.4617080688,7.51789474487
     5248,97.458366394,7.55314159393
     5349,97.4421920776,7.58513689041
     5450,97.4163894653,7.61389112473
     5551,97.4002227783,7.6426782608
     5652,97.3776092529,7.67464590073
     5753,97.3774642944,7.71310758591
     5854,97.370880127,7.7547492981
     5955,97.3354415894,7.77704811096
     6056,97.3032531738,7.78974056244
     6157,97.2838439941,7.8217124939
     6258,97.2708435059,7.85691452026
     6359,97.2354049683,7.87599658966
     6460,97.1999816895,7.89187002182
     6561,97.1805419922,7.92703914642
     6662,97.1771621704,7.96548223495
     6763,97.1705703735,8.00070571899
     6864,97.1607589722,8.03911972046
     6965,97.1381149292,8.06785964966
     7066,97.1026611328,8.08692550659
     7167,97.0704269409,8.10279750824
     7268,97.0477600098,8.13473510742
     7369,97.031539917,8.16349697113
     7470,97.0121231079,8.18903827667
     7571,96.9894638062,8.21776580811
     7672,96.9764633179,8.24333572388
     7773,96.969871521,8.27534675598
     7874,96.9793243408,8.31384754181
     7975,96.9759216309,8.34907817841
     8076,96.9532165527,8.38420963287
     8177,96.9337768555,8.40974235535
     8278,96.9111328125,8.4320526123
     8379,96.8820114136,8.46073436737
     8480,96.8465270996,8.47656059265
     8581,96.8079605103,8.4667339325
     8682,96.7693481445,8.46651554108
     8783,96.7988052368,8.93318271637
     8885,98.3799972534,9.22999954224
     8984,98.5999984741,9.94999980927
     9087,98.8700027466,8.06999969482
     9186,98.4199981689,7.82999992371
     9293,98.2940063477,7.90318584442
     9392,98.2513580322,7.90388059616
     9488,99.3000030518,7.5
     9589,99.6500015259,6.69999980927

  • anuga_work/development/boxingday08/project.py

    r5458 r5465  
    22from Scientific.IO.NetCDF import NetCDFFile
    33from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
    4     compress,zeros,fabs,allclose
     4    compress,zeros,fabs,allclose,ones
    55from anuga.utilities.polygon import inside_polygon,read_polygon
    66from os import sep
    77from time import localtime, strftime, gmtime
    8 
    9 def create_sts_boundary_order_file(order_filename,stsfilename,verbose=False):
     8from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
     9import os
     10
     11def create_sts_boundary_order_file(urs_filename,order_filename,sts_filename,verbose=False):
    1012    """
    1113    Note This is not a robust algorithm. Be Careful when applying that
    … …  
    1315    follows the path it is supposed too.
    1416    """
    15     fid = NetCDFFile(stsfilename+'.sts', 'r')    #Open existing file for read
    16     x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
    17     y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
    18     coordinates=transpose(asarray([x.tolist(),y.tolist()]))
     17    from anuga.shallow_water.data_manager import read_mux2_py
     18    times, y, x, elevation, quantity, starttime = read_mux2_py([urs_filename],weights=ones(1,Float))
    1919   
    2020    #find indices of the southernmost gauges
    … …  
    6060             
    6161    #So far i cannot tell if the wrong point is closer.
     62   
     63    boundary=ensure_numeric(boundary)
     64    from anuga.coordinate_transforms import redfearn,convert_from_latlon_to_utm
     65    boundary_utm,zone=convert_from_latlon_to_utm(longitudes=list(boundary[:,0]),latitudes=list(boundary[:,1]))
     66
    6267    d=","
    6368    fid=open(order_filename,'w')
    64     header="index,longitude,latitude\n"
     69    header='index, longitude, latitude\n'
    6570    fid.write(header)
    6671    line=str(west_index)+d+str(x[west_index])+d+\
    … …  
    7277            fid.write(line)
    7378    fid.close()
    74     return boundary
    75            
    76 
    77 #read in polyline boundary
    78 from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
    79 urs_name='polyline'
    80 #base_name=urs_name
    81 tide=0.35
    82 base_name='tide_polyline'
    83 #base_name='most_polyline'
    84 import os
    85 if os.path.exists(base_name+'.sts'):
    86     print 'sts boundary file already created.'
    87 else:
    88     print 'creatin sts file'
    89     urs2sts(urs_name,base_name,mean_stage=tide,verbose=False)
    90 
    91 order_filename='sts_order.txt'
    92 boundary=create_sts_boundary_order_file(order_filename,base_name,verbose=False)
    93 urs_boundary=create_sts_boundary(order_filename,base_name,lat_long=False)
    94 assert allclose(boundary.tolist(),urs_boundary)
    95 
    96 #Read in large domain boundary
    97 dir='mesh_polygons'
    98 extent = 'extent.csv'
    99 bp = read_polygon(dir+sep+extent)
    100 
    101 #Now order of boundary is correct clip sts boundary to small region for
    102 #this particular scenario
    103 indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False)
    104 bounding_polygon=[]
    105 for i in indices[:-2]:
    106     #disregard last two entries because they are not on boundary
    107     bounding_polygon.append(urs_boundary[i])
    108 
    109 #Append clipped sts boundary to large domain boundary
    110 bounding_polygon.extend(bp[2:len(bp)])
    111 
    112 ###############################
    113 # Domain definitions
    114 ###############################
     79    return boundary,boundary_utm
     80
     81######################
     82# Domain definitions #
     83######################
    11584
    11685# Bathymetry and topography filenames
    … …  
    12998extent = 'extent.csv'
    13099patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')
     100   
     101tide=0.35
     102
     103######################################
     104# Create urs boundary from mux2files #
     105######################################
     106urs_filename='polyline'
     107base_name='tide_polyline'
     108order_filename='boxingday_boundary_order.txt'
     109
     110# Create ordering file
     111# Usually this would be provided before hand
     112boundary,boundary_utm=create_sts_boundary_order_file(urs_filename+'_waveheight-z-mux2',order_filename,base_name,verbose=False)
     113
     114# Create ordered sts file
     115if not os.path.exists(base_name+'.sts'):
     116    print 'creating sts file'
     117    urs2sts(urs_filename,basename_out=base_name,
     118    ordering_filename=order_filename,
     119    mean_stage=tide,
     120    verbose=False)
     121else:
     122    print 'sts file exists'
     123
     124# Read in boundary from ordered sts file
     125urs_boundary=create_sts_boundary(base_name)
     126
     127# Checks that the boundary read from the sts file is the same as the one
     128# defined by the function create_sts_boundary_order_file
     129assert allclose(boundary_utm,urs_boundary)
     130
     131#############################################################
     132# Add urs boundary segment to user defined boundary polygon #
     133#############################################################
     134
     135# Read in specific user defined boundary
     136dir='mesh_polygons'
     137extent = 'extent.csv'
     138user_boundary_polygon = read_polygon(dir+sep+extent)
     139
     140# The following is specific to this scenario
     141indices=inside_polygon(urs_boundary, user_boundary_polygon, closed=True, verbose=False)
     142bounding_polygon=[]
     143for i in indices[:-2]:
     144    #disregard last two entries because they are not on boundary
     145    bounding_polygon.append(urs_boundary[i])
     146
     147# Append clipped sts boundary to user defined domain boundary
     148bounding_polygon.extend(user_boundary_polygon[2:len(user_boundary_polygon)])
    131149
    132150###############################
    133 # Interior region definitions
     151# Interior region definitions #
    134152###############################
    135153
    … …  
    141159patong = read_polygon(dir+sep+'patong.csv')
    142160bay = read_polygon(dir+sep+'bay.csv')
     161
     162##########################################
     163# Plot urs boundary and internal regions #
     164##########################################
    143165
    144166plot=False
    … …  
    156178    show()
    157179
     180   
     181#######################
     182# For MOST SIMULATION #
     183#######################
    158184east  = 97.7
    159185west  = 96.7
Note: See TracChangeset for help on using the changeset viewer.