# Convert ASC 2 DEM 2 PTS using source data and store result in source data
# Do for coarse and fine data
# Fine pts file to be clipped to area of interest
#-------------------------------------------------------------------------------
print"project.combined_dir_name",project.combined_dir_name

# topography directory filenames
onshore_dir_name = project.onshore_dir_name
onshore_dir_name1 = project.onshore_dir_name
#island_in_dir_name = project_urs.island_in_dir_name
Multibeam_dir_name = project.Multibeam_dir_name
added_data_dir_name =project.added_data_dir_name

print'create Geospatial data1 objects from topographies',onshore_dir_name + '.txt'
G1 = Geospatial_data(file_name = onshore_dir_name + '.txt')

print'create Geospatial data1 objects from topographies',onshore_dir_name1 + '.txt'
G2 = Geospatial_data(file_name = onshore_dir_name1 + '.txt')

print'create Geospatial data2 objects from multibeam', Multibeam_dir_name + '.txt'
G_off1 = Geospatial_data(file_name = Multibeam_dir_name + '.txt')

print'create Geospatial data objects from added data',added_data_dir_name + '.txt'
G_off2 = Geospatial_data(file_name = added_data_dir_name + '.txt')

print'add all geospatial objects'
Gt1 = G1.clip(Geospatial_data(project.poly_west_reef))+G1.clip(Geospatial_data(project.poly_west))+G1.clip(Geospatial_data(project.poly_centerKolonia))+G1.clip(Geospatial_data(project.poly_east))
Gt2= G2.clip_outside(Geospatial_data(project.poly_centerKolonia)).clip_outside(Geospatial_data(project.poly_east))
Gt3= G_off1.clip_outside(Geospatial_data(project.poly_west_reef)).clip_outside(Geospatial_data(project.poly_west)).clip_outside(Geospatial_data(project.poly_centerKolonia)).clip_outside(Geospatial_data(project.poly_east))
                         
G=Gt1+Gt2+Gt3+G_off2
                         
#print'clip combined geospatial object by bounding polygon'
#G_clipped = G.clip(project_urs.poly_all)
#FIXME: add a clip function to pts
#print'shape of clipped data', G_clipped.get_data_points().shape

print'export combined DEM file'
if access(project.topographies_dir,F_OK) == 0:
    mkdir (project.topographies_dir)
print'export',project.combined_dir_name+ '.txt'
G.export_points_file(project.combined_dir_name+ '.txt')



'''
print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
print'split'
G_all_1, G_all_2 = G_all.split(.10)
print'export 1'
G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
print'export 2'
G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')


#-------------------------------------------------------------------------
# Convert URS to SWW file for boundary conditions 
#-------------------------------------------------------------------------
print 'starting to create boundary conditions'
boundaries_in_dir_name = project.boundaries_in_dir_name

from anuga.shallow_water.data_manager import urs2sww

print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary

#import sys; sys.exit()

#if access(project.boundaries_dir,F_OK) == 0:
#    mkdir (project.boundaries_dir)

from caching import cache
cache(urs2sww,
      (boundaries_in_dir_name,
       project.boundaries_dir_name1), 
      {'verbose': True,
       'minlat': project.south_boundary,
       'maxlat': project.north_boundary,
       'minlon': project.west_boundary,
       'maxlon': project.east_boundary,
#       'minlat': project.south,
#       'maxlat': project.north,
#       'minlon': project.west,
#       'maxlon': project.east,
       'mint': 0, 'maxt': 40000,
#       'origin': domain.geo_reference.get_origin(),
       'mean_stage': project.tide,
#       'zscale': 1,                 #Enhance tsunami
       'fail_on_NaN': False},
       verbose = True,
       )
#       dependencies = source_dir + project.boundary_basename + '.sww')

'''