Skip to content

Commit

Permalink
Added fileintype to fort14togdf
Browse files Browse the repository at this point in the history
  • Loading branch information
jmpmcmanus committed Aug 1, 2023
1 parent 9d85155 commit 451be23
Showing 1 changed file with 24 additions and 1 deletion.
25 changes: 24 additions & 1 deletion kalpana/export.py
Original file line number Diff line number Diff line change
Expand Up @@ -878,7 +878,7 @@ def nc2shp(ncFile, var, levels, conType, pathOut, epsgOut, vUnitOut='ft', vUnitI
logger.info(f'Ready with exporting code after: {(time.time() - t00)/60:0.3f} min') # Changed
return gdf

def fort14togdf(filein, epsgIn, epsgOut):
def fort14togdf(filein, epsgIn, epsgOut, fileintype='netcdf'):
''' Write adcirc mesh from fort.14 file as GeoDataFrame and extract centroid of each element.
Used in the downscaling process
Parameters:
Expand All @@ -888,12 +888,15 @@ def fort14togdf(filein, epsgIn, epsgOut):
coordinate system of the adcirc input
epsgOut: int
coordinate system of the output shapefile
fileInType: str
input file type, which are netcdf (default) or fort.14
Returns
gdf: GeoDataFrame
GeoDataFrame with polygons as geometry and more info such as: area, representative
element size, centroids coordinates, and vertices
'''
## read only the two first lines of the file to get the number of elements and nodes
'''
with open(filein) as fin:
head = list(islice(fin, 2))
data = [int(x) for x in head[1].split()]
Expand All @@ -904,6 +907,26 @@ def fort14togdf(filein, epsgIn, epsgOut):
x = nodes[:, 0]
y = nodes[:, 1]
z = nodes[:, 2]
'''
if fileintype == 'fort.14':
## read only the two first lines of the file to get the number of elements and nodes
with open(filein) as fin:
head = list(islice(fin, 2))
data = [int(x) for x in head[1].split()]
## read nodes
nodes = np.loadtxt(filein, skiprows = 2, max_rows = data[1], usecols = (1, 2, 3))
## read elements
elem = np.loadtxt(filein, skiprows = 2 + data[1], max_rows = data[0], usecols = (2, 3, 4)) - 1
x = nodes[:, 0]
y = nodes[:, 1]
z = nodes[:, 2]
elif fileintype == 'netcdf':
with netcdf.Dataset(filein) as nc:
x = nc['x'][:].data
y = nc['y'][:].data
z = nc['depth'][:].data
elem = nc['element'][:].data - 1

## matplotlib triangulation
tri = mpl.tri.Triangulation(x, y, elem)
## select the coordinate of each vertex
Expand Down

0 comments on commit 451be23

Please sign in to comment.