... | ... | @@ -131,6 +131,40 @@ Exit properly |
|
|
```python
|
|
|
slimPre.exit(0)
|
|
|
```
|
|
|
### Define a Manning coefficient based on polygons from shapefiles
|
|
|
This example requires ogr from the python-gdal package.
|
|
|
```
|
|
|
import ogr
|
|
|
import slimPre
|
|
|
import sys
|
|
|
import numpy as np
|
|
|
|
|
|
def get_index_of_points_in_polygons(points,polygonshp) :
|
|
|
ds = ogr.Open(polygonshp)
|
|
|
lyr_in = ds.GetLayer(0)
|
|
|
point_ref=ogr.osr.SpatialReference()
|
|
|
point_ref.ImportFromProj4("+proj=latlong +ellps=WGW84")
|
|
|
print(lyr_in.GetSpatialRef().ExportToProj4())
|
|
|
ctran=ogr.osr.CoordinateTransformation(point_ref,lyr_in.GetSpatialRef())
|
|
|
def check(x,y):
|
|
|
[lon,lat,z]=ctran.TransformPoint(x,y)
|
|
|
pt = ogr.Geometry(ogr.wkbPoint)
|
|
|
pt.SetPoint_2D(0, lon,lat)
|
|
|
lyr_in.SetSpatialFilter(pt)
|
|
|
return len(lyr_in) > 0
|
|
|
mask = [i for i in range(points.shape[0]) if check(points[i,0],points[i,1])]
|
|
|
return mask
|
|
|
|
|
|
mesh = slimPre.Mesh("torres.msh",'+proj=utm +ellps=WGS84 +zone=54 +south')
|
|
|
region = slimPre.Region(mesh)
|
|
|
points = slimPre.Coordinate_system(region,"+proj=longlat +ellps=WGS84").coordinates*180/np.pi
|
|
|
|
|
|
data = np.zeros([points.shape[0]])
|
|
|
data[:] = 0.023
|
|
|
data[get_index_of_points_in_polygons(points,"SourceLocationsSLIM/Original_Seagrass_Polygon_UTM54s.shp")] = 0.055
|
|
|
data[get_index_of_points_in_polygons(points,"SourceLocationsSLIM/ReefsDry.shp")] = 0.072
|
|
|
slimPre.write_file("manning.nc",region=region,data=[("manning",data)])
|
|
|
```
|
|
|
|
|
|
## Command line
|
|
|
|
... | ... | |