testshp.py 954 Bytes
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
3
4
import msea
import numpy as np
from osgeo import osr 

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
6
7
def mesh_size(x,projection) :
    s_coast = np.clip((dist_coast_field(x,projection)-400)*0.5,200,5000)
    s_porquerolles = np.clip((dist_porquerolles_field(x,projection)-200)*0.5,50,5000)
Jonathan Lambrechts's avatar
fixes    
Jonathan Lambrechts committed
8
    return np.minimum(s_coast,s_porquerolles)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
9
10

domain_srs = osr.SpatialReference()
Jonathan Lambrechts's avatar
fixes    
Jonathan Lambrechts committed
11
12
13

domain_srs.ImportFromEPSG(32631)
#domain_srs.ImportFromProj4("+proj=utm +ellps=WGS84 +zone=31")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
14
15

domain = msea.Domain(domain_srs)
16
domain.add_shapefile("test/data_no_duplicate.shp","physical",curve_type=msea.POLYLINE)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
#domain.add_shapefile("test/interior.shp",None,True,curve_type=msea.POLYLINE)
18
#domain.add_interior_points("test/interior.shp")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
19
bath_field = msea.RasterField("medit.tiff")
Jonathan Lambrechts's avatar
fixes    
Jonathan Lambrechts committed
20
21
22
dist_coast_field = msea.DistanceField(domain,100,["coast","island"])
dist_porquerolles_field = msea.DistanceField(domain,20,["porquerolles"])
domain.unrefine_boundaries(8e5,4.68e6,20,mesh_size)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
23
msea.mesh_gmsh(domain,"test.msh",mesh_size)