issue calculating transport

[LAWTIDEBOY]

• pyfesom2 version: 0.2.0
• Python version: 3.7.4
• Operating System: HLRN

Description

I am trying to reproduce the tutorial for transport calculation. It looks like there is an array dimension issue, but I cant find the problem in transport.py

What I Did

import warnings
warnings.filterwarnings('ignore')

%autosave 5
#%matplotlib widget
import sys
import os
import numpy as np
from netCDF4 import Dataset
import pandas as pd
import xarray as xr
sys.path.append('/home/hbkoziel/pyfesom2/codes/py_f2recom/modules')
#sys.path.append('/home/hbkoziel/pyfesom2/pyfesom2')
import pyfesom2 as pf
import numba as nb
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from matplotlib.dates import DateFormatter
sys.path.append('/home/hbkoziel/pyfesom2/codes/py_f2recom/modules/cmocean-master/')
import cmocean as cmo

resultpath = '/scratch/projects/hbk00083/model_outputs/fesom2.1_recom/'
meshpath = '/scratch/usr/hbkoziel/mesh/farc'
dfpaths = '/home/hbkoziel/pyfesom2/codes/py_f2recom_develop/Erosion/Dataframes/'
runid='fesom'
meshdiag=meshpath+'/'+runid+'.mesh.diag.nc'

mesh = pf.load_mesh(meshpath)
diag = pf.get_meshdiag(mesh,meshdiag=meshdiag, runid=runid)
nod_area = diag.nod_area
mesh.path = meshpath

transport_85N, section_85N = pf.cross_section_transport([-179.8, 180, 85, 85],                                                   # select a section from the presets or [lon1, lon2, lat1, lat2]
                                                mesh=mesh,            # mesh
                                                data_path=resultpath+'EXP1/',   # directory of the u, v, (T, S, uice, vice, m_ice, a_ice) files
                                                mesh_diag=diag,                                            # mesh_diag_file
                                                years=np.arange(2000,2001),                                     # years to compute
                                                use_great_circle=False,                                  # compute the section as a great circle
                                                how='ori',                                               # 'ori' do not apply mean, 'mean' apply time mean
                                                add_TS=False,                                            # add temperature and salinity to the section
                                                add_extent=.075                                             # the extent to look for gridcells nearby the section, choose large for low resolutions
                                               )

OUTPUTS:

Javascript Error: IPython is not defined
Autosaving every 5 seconds
osgeo is not installed, conversion to Geo formats like Geotiff (fesom2GeoFormat) will not work.
/scratch/usr/hbkoziel/mesh/farc/pickle_mesh_py3_fesom2
The usepickle == True)
The pickle file for FESOM2 exists.
The mesh will be loaded from /scratch/usr/hbkoziel/mesh/farc/pickle_mesh_py3_fesom2

Starting computation...

Your section:  not specified : Start:  -179.8 °E  85 °N  End:  180 °E  85 °N

Converting grid cells to Polygons... (If this takes very long try to reduce the add_extent parameter)
HBox(children=(IntProgress(value=0, max=3554), HTML(value='')))

Looking for intersected grid cells...
HBox(children=(IntProgress(value=0, max=3554), HTML(value='')))

ERROR MESSAGE:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-15-41fbc66c174c> in <module>
      7                                                 how='ori',                                               # 'ori' do not apply mean, 'mean' apply time mean
      8                                                 add_TS=False,                                            # add temperature and salinity to the section
----> 9                                                 add_extent=.075                                             # the extent to look for gridcells nearby the section, choose large for low resolutions
     10                                                )

~/pyfesom2/codes/py_f2recom/modules/pyfesom2/transport.py in cross_section_transport(section, mesh, data_path, years, mesh_diag, how, add_extent, abg, add_TS, chunks, use_great_circle, n_points)
   1087     elem_box_nods, elem_box_indices = _ReduceMeshElementNumber(section_waypoints, mesh, section, add_extent)
   1088 
-> 1089     elem_box_nods, elem_box_indices, cell_intersections, line_section = _LinePolygonIntersections(mesh, section_waypoints, elem_box_nods, elem_box_indices)
   1090 
   1091     intersected_edge, midpoints_edge, elem_centers, elem_box_indices, elem_box_nods, cell_intersections = _FindIntersectedEdges(mesh, elem_box_nods, elem_box_indices, line_section, cell_intersections)

~/pyfesom2/codes/py_f2recom/modules/pyfesom2/transport.py in _LinePolygonIntersections(mesh, section_waypoints, elem_box_nods, elem_box_indices)
    403     # remove indices of elements that are not intersected
    404     elem_box_nods = elem_box_nods[intersection_bool]
--> 405     elem_box_indices = elem_box_indices[intersection_bool]
    406 
    407     return elem_box_nods, elem_box_indices, cell_intersections, line_section

IndexError: boolean index did not match indexed array along dimension 0; dimension is 3555 but corresponding boolean dimension is 3554

[LAWTIDEBOY]

I would like to add the following information: I am using the fArc mesh :)

Thank in advance for the support !!

[koldunovn]

@LAWTIDEBOY My first suspicion is transposed output. In the output data that you use, what is the order of dimensions [time, elem, depth] or [time, depth, elem]?

[LAWTIDEBOY]

Hello Nikolay, thank !

ncdump -h u.fesom.2000.nc

gives this:

float u(time, elem, nz1) ; u:description = "horizontal velocity" ; u:long_name = "horizontal velocity" ; u:units = "m/s" ;

so: [time, elem, depth] is that wrong ?

[koldunovn]

Both are right, but produced by different versions of the model :) Your variant is the one that the transport code was initially developed for. However the problem was handled in this PR https://github.com/FESOM/pyfesom2/pull/176 . Can you send me your data (subset) and notebook that you trying to run?

[koldunovn]

@LAWTIDEBOY can you please try #192 and see if this fixes your problem. Sorry for the delay :(

[LAWTIDEBOY]

Thank Nikolay, I checked and it works fine !!! :))))