Create boundary and initial file from seagrid and GLBa0.08

[bangdt]

Hi everyone, please help me to create boundary and initial file from seagrid and GLBa0.08 - Hycom data.

I created seagrid nc file for my domain and have Hycom data. I also install pyroms_toolbox. My problem is I have no idea about how to make boundary and initial file. Could anybody help me a example by using pyroms_toolbox to do that please? especially how to fill in the file to run

My questions are

1) in seagrid how can I chose the number of vertical cell? (if I want my grid size is x=30; y=30; z=15)

2) for Grid_HYCOM, get_nc_Grid_HYCOM and make_remap_grid_file which information I should use in?
------------------------------------------------------------------------------------
class Grid_HYCOM(object):
"""
Grid object for HYCOM
"""

def __init__(self, lon_t, lat_t, lon_vert, lat_vert, mask_t, z_t, h, angle, name):

self.name = name

self.lon_t = lon_t
self.lat_t = lat_t

self.lon_vert = lon_vert
self.lat_vert = lat_vert

self.mask_t = mask_t

self.z_t = z_t

self.h = h

self.angle = angle
-----------------------------------------------------------------------------
3) should I run only one of flood, flood_fast and flood_fast_weighted or all of them?

Regard.

[kate]

I have successfully made ROMS IC and BC files from HYCOM using the scripts in pyroms/examples/Arctic_HYCOM/. First you run make_remap_weights_file, then make_ic_file/make_bdry_file.
A more equatorial example is in Palau_HYCOM.

For the vertical, set the number of vertical points in the gridid.txt file. My ARCTIC4 chunk looks like:

Code: Select all

id      = ARCTIC4
name    = ARCTIC4
grdfile = /import/AKWATERS/kshedstrom/gridpak/Arctic2/grid_Arctic_4.nc
N       = 50
grdtype = roms
Vtrans  = 4
theta_s = 7
theta_b = 2
Tcline  = 250

[bangdt]

. First you run make_remap_weights_file, then make_ic_file/make_bdry_file.
A more equatorial example is in Palau_HYCOM.

Dear Kate

I download Hycom data with empty z (to import all z to the netcdf file)
then try to run make_remap_weights_file and got problem with this
------------------------------------------------------------------

Code: Select all

>>> srcgrd = pyroms_toolbox.Grid_HYCOM.get_nc_Grid_HYCOM('/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/Hycomdata/1906201200.nc')

Traceback (most recent call last):
  File "<pyshell#6>", line 1, in <module>
    srcgrd = pyroms_toolbox.Grid_HYCOM.get_nc_Grid_HYCOM('/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/Hycomdata/1906201200.nc')
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/Grid_HYCOM/get_nc_Grid_HYCOM.py", line 19, in get_nc_Grid_HYCOM
    depth = nc.variables['z'][:]
KeyError: 'z'

------------------------------------------------------------------
by open the hycom netcdf file data by IDV and I have all z data can be visualized in 3d model so I do not understand why nc.variables['z'][:] have that error?

by running ncdump -h /work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/Hycomdata/1906201200.nc > my_file.cdl I got 40 depth levels like this.
-----------------------------------------------------------------

Code: Select all

netcdf \1906201200 {
dimensions:
	time = 5 ;
	lat = 626 ;
	lon = 313 ;
	depth = 40 ;
variables:
...
float salinity(time, depth, lat, lon) ;
		salinity:_CoordinateAxes = "time_run time depth lat lon" ;
		salinity:units = "psu" ;
		salinity:long_name = "Salinity" ;
		salinity:standard_name = "sea_water_salinity" ;
		salinity:NAVO_code = 16 ;
		salinity:coordinates = "time_run time depth lat lon " ;
	double depth(depth) ;
		depth:units = "m" ;
		depth:long_name = "Depth" ;
		depth:standard_name = "depth" ;
		depth:positive = "down" ;
		depth:axis = "Z" ;
		depth:NAVO_code = 5 ;
		depth:_CoordinateAxisType = "Height" ;
		depth:_CoordinateZisPositive = "down" ;
...

----------------------------------------------------------------
my file /get_nc_Grid_HYCOM.py", line 19 got

Code: Select all

    depth = nc.variables['Z'][:]
#   depth = nc.variables['z'][:]

I tried booth Z and z with same error
-----------------------------------------------------------------
Thank you so much for your allway very quick and helpful replies.

Regard.

[kate]

Looks like it should be 'depth', not 'z'.

[bangdt]

Looks like it should be 'depth', not 'z'.

Dear Kate
you are right I changed it to depth and pass, and have this error after that

Code: Select all

>>> srcgrd = pyroms_toolbox.Grid_HYCOM.get_nc_Grid_HYCOM('/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/Hycomdata/1906201200.nc')

Traceback (most recent call last):
  File "<pyshell#4>", line 1, in <module>
    srcgrd = pyroms_toolbox.Grid_HYCOM.get_nc_Grid_HYCOM('/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/Hycomdata/1906201200.nc')
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/Grid_HYCOM/get_nc_Grid_HYCOM.py", line 21, in get_nc_Grid_HYCOM
    var = nc.variables['sea_water_temperature'][0,:,1:-1,1:-1]
KeyError: 'sea_water_temperature'

my cdl file have this information

Code: Select all

float water_temp(time, depth, lat, lon) ;
		water_temp:_CoordinateAxes = "time_run time depth lat lon" ;
		water_temp:units = "degC" ;
		water_temp:long_name = "Water Temperature" ;
		water_temp:standard_name = "sea_water_temperature" ;
		water_temp:NAVO_code = 15 ;
		water_temp:comment = "in-situ temperature" ;
		water_temp:coordinates = "time_run time depth lat lon " ;

So what should I need to replace for "sea_water_temperature" ?
Change it to water_temp and passed

then I got this

Code: Select all

>>> srcgrd = pyroms_toolbox.Grid_HYCOM.get_nc_Grid_HYCOM('/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/Hycomdata/1906201200.nc')

Traceback (most recent call last):
  File "<pyshell#4>", line 1, in <module>
    srcgrd = pyroms_toolbox.Grid_HYCOM.get_nc_Grid_HYCOM('/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/Hycomdata/1906201200.nc')
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/Grid_HYCOM/get_nc_Grid_HYCOM.py", line 24, in get_nc_Grid_HYCOM
    lon_t = lon[1:-1,1:-1]
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/numpy/ma/core.py", line 3174, in __getitem__
    dout = self.data[indx]
IndexError: too many indices for array

My installed pyroms is for python 2.7 with included Hycom folder with scripts so I hope these scripts are wrote for python 2.7 but it seem that "IndexError: too many indices for array" related to version conflict doesn't it?


Thank you
Regard.

[kate]

No, I don't think you have a version problem, unless it's a HYCOM version problem.

Code: Select all

lon_t = lon[1:-1,1:-1]

Longitude is often a 1-D array, not a 2-D array. What do you have?

[bangdt]

No, I don't think you have a version problem, unless it's a HYCOM version problem.

Code: Select all

lon_t = lon[1:-1,1:-1]

Longitude is often a 1-D array, not a 2-D array. What do you have?

Dear Kate
I am sorry, I still do not understand my problem.
I would like to make IC BC from this data https://ncss.hycom.org/thredds/ncss/gri ... taset.html
and I have get_nc_Grid_HYCOM.py code after install pyroms like this

Code: Select all

import numpy as np
import pyroms
import netCDF4
from mpl_toolkits.basemap import pyproj
from pyroms_toolbox.Grid_HYCOM import Grid_HYCOM


def get_nc_Grid_HYCOM(grdfile, name='GLBa0.08_NEP'):

    """
    grd = get_nc_Grid_HYCOM(grdfile)

    Load grid object for HYCOM_GLBa0.08_NEP
    """

    nc = netCDF4.Dataset(grdfile)
    lon = nc.variables['lon'][:]
    lat = nc.variables['lat'][:]
    depth = nc.variables['z'][:]
    var = nc.variables['temp'][0,:,1:-1,1:-1]
    nc.close()

    lon_t = lon[1:-1,1:-1]
    lat_t = lat[1:-1,1:-1]

    lon_vert = 0.5 * (lon[:,1:] + lon[:,:-1])
    lon_vert = 0.5 * (lon_vert[1:,:] + lon_vert[:-1,:])

    lat_vert = 0.5 * (lat[1:,:] + lat[:-1,:])
    lat_vert = 0.5 * (lat_vert[:,1:] + lat_vert[:,:-1])

    mask_t = np.array(~var[:].mask, dtype='int')

    z_t = np.tile(depth,(mask_t.shape[2],mask_t.shape[1],1)).T

    depth_bnds = np.zeros(len(depth)+1)
    for i in range(1,len(depth)):
        depth_bnds[i] = 0.5 * (depth[i-1] + depth[i])
    depth_bnds[-1] = 5750

    bottom = pyroms.utility.get_bottom(var[::-1,:,:], mask_t[0], spval=var.fill_value)
    nlev = len(depth)
    bottom = (nlev-1) - bottom
    h = np.zeros(mask_t[0,:].shape)
    for i in range(mask_t[0,:].shape[1]):
        for j in range(mask_t[0,:].shape[0]):
            if mask_t[0,j,i] == 1:
                h[j,i] = depth_bnds[int(bottom[j,i])+1]


    geod = pyproj.Geod(ellps='WGS84')
    az_forward, az_back, dx = geod.inv(lon_vert[:,:-1], lat_vert[:,:-1], lon_vert[:,1:], lat_vert[:,1:])
    angle = 0.5 * (az_forward[1:,:] + az_forward[:-1,:])
    angle = (90 - angle) * np.pi/180.


    return Grid_HYCOM(lon_t, lat_t, lon_vert, lat_vert, mask_t, z_t, h, angle, name)

Thank you
Regard

[kate]

Are you sure you are downloading on the GLBa0.08 HYCOM grid? They have moved that to the "inactive" status. It was on a tripole grid and needed lat,lon to be functions if both i and j. The newer HYCOM products are on a simpler grid, where lat is a function of j only and lon is a function of i only. I'm afraid I haven't tried to run these scripts with the new HYCOM products. It's on my list of things to do, just not necessarily this week.

[bangdt]

Are you sure you are downloading on the GLBa0.08 HYCOM grid? They have moved that to the "inactive" status. It was on a tripole grid and needed lat,lon to be functions if both i and j. The newer HYCOM products are on a simpler grid, where lat is a function of j only and lon is a function of i only. I'm afraid I haven't tried to run these scripts with the new HYCOM products. It's on my list of things to do, just not necessarily this week.

Dear Kate
Thank you for your suggestion. I double check and found out I have downloaded wrong data for GLBa0.08. Actually I would like to have that data for predict/ forecasting works but at the moment, for just practice using ROMS, It is better for me if I keep up with GLBa0.08 examples to understand how pyroms work.

I followed instruction in this topic viewtopic.php?f=23&t=4819
mdessert wrote

1-I've well downloaded files from HYCOM.org.
2-I've formatted them through get_hycom_GLBa0.08_'varname'_2014.py.
3-I've built my source grid through get_hycom_GLBa0.08_grid.py.
4-I've computed my weigth files from make_remap_weights_file.py.
5-And then, I've 'concatenated' all my variable files to obtain one big file HYCOM_GLBa0.08_2014_095.nc.

at 1- I am not sure how he downloaded these files, I try to run 2- then get the *.nc file for salt, ssh, temp, u and v

then I run 3- once time and get the GLBa0.08 grid.nc file
when I run 4- for the 1st time I stuck as mention bellow.
I do not know how to run step 5- to have all in one big file HYCOM_GLBa0.08_2014_095.nc


I stuck at running make_remap_weights_file.py :

Code: Select all

Traceback (most recent call last):
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/make_remap_weights_file.py", line 7, in <module>
    srcgrd = pyroms_toolbox.Grid_HYCOM.get_nc_Grid_HYCOM('/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/HYCOM_GLBa0.08_PALAU_grid.nc')
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/Grid_HYCOM/get_nc_Grid_HYCOM.py", line 52, in get_nc_Grid_HYCOM
    h[j,i] = depth_bnds[bottom[j,i]+1]
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

I am sorry because I am not a code/programmer then It should take me time to pass basic questions like this.
Regard and have a nice weekend.

[kate]

My latest version of that troublesome line is:

Code: Select all

h[j,i] = depth_bnds[int(bottom[j,i])+1]

Sometimes it pays to "git pull".

Step 5 happens with this script:

Code: Select all

#!/bin/ksh

print -n "Enter year to process: "; read year; print ""

leap=`echo $(($year % 4))`

if [ $leap == 0 ] ; then
    nday=366
else
    nday=365
fi
#nday=31

#set -A days {1..331}
#set -A days {2..4}
#set -A days {33..$nday}
set -A days {1..$nday}

for day in ${days[@]} ; do
    day=`echo $day | awk '{printf "%03d", $1}'`
    ncks -a -O HYCOM_GLBa0.08_ssh_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_temp_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_salt_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_u_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_v_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    rm -f HYCOM_GLBa0.08_ssh_${year}_${day}.nc HYCOM_GLBa0.08_temp_${year}_${day}.nc \ HYCOM_GLBa0.08_salt_${year}_${day}.nc HYCOM_GLBa0.08_u_${year}_${day}.nc \ 
HYCOM_GLBa0.08_v_${year}_${day}.nc
done

Gosh, should convert to bash.

[bangdt]

My latest version of that troublesome line is:

Code: Select all

h[j,i] = depth_bnds[int(bottom[j,i])+1]

Sometimes it pays to "git pull".

Dear Kate
by looking in get_nc_Grid_HYCOM.py I have

Code: Select all

for i in range(1,len(depth)):
        depth_bnds[i] = 0.5 * (depth[i-1] + depth[i])
    depth_bnds[-1] = 5750

It suggested that after that, depth_bnds become not integer anymore (0.5*....), then could it be leading to the "IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices" error?
Regard.

[kate]

Sometimes Python evolves so that something that used to run no longer runs. Sometimes when that happens I can figure out how to change it and get it pushed to github. You need to update your code through either "git pull" or just editing it in place.

[bangdt]

Sometimes Python evolves so that something that used to run no longer runs. Sometimes when that happens I can figure out how to change it and get it pushed to github. You need to update your code through either "git pull" or just editing it in place.

Dear Kate
Thank you so much, I do git pull and updated the code, then pyroms go to this

Code: Select all

Traceback (most recent call last):
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/make_remap_weights_file.py", line 8, in <module>
    dstgrd = pyroms.grid.get_ROMS_grid('PALAU1')
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms/grid.py", line 467, in get_ROMS_grid
    gridinfo = ROMS_gridinfo(gridid,hist_file=hist_file,grid_file=grid_file)
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms/grid.py", line 71, in __init__
    self._get_grid_info(grid_file,hist_file)
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms/grid.py", line 82, in _get_grid_info
    data = open(gridid_file,'r')
TypeError: coercing to Unicode: need string or buffer, NoneType found

I think it is related to gridid.txt. My gridid.txt has no infor about PALAU1 (I used script in PALAU_HYCOM in example folder. So how can I find information to fill in gridid.txt to run?

Regard.

[kate]

I can give you the PALAU1 chunk of gridid.txt, but then it points to a grid file which you don't have either. Why don't you put in a chunk for your grid, since that is what you want anyway?

Code: Select all

id      = PALAU1
name    = PALAU1
grdfile = /import/AKWATERS/kshedstrom/gridpak/Palau/grid_Palau_1.nc
N       = 50
grdtype = roms
Vtrans  = 4
theta_s = 7
theta_b = 2
Tcline  = 250

This is the PALAU chunk which has the vertical coordinates we use for most domains these days.

[bangdt]

My latest version of that troublesome line is:

Code: Select all

h[j,i] = depth_bnds[int(bottom[j,i])+1]

Sometimes it pays to "git pull".

Step 5 happens with this script:

Code: Select all

#!/bin/ksh

print -n "Enter year to process: "; read year; print ""

leap=`echo $(($year % 4))`

if [ $leap == 0 ] ; then
    nday=366
else
    nday=365
fi
#nday=31

#set -A days {1..331}
#set -A days {2..4}
#set -A days {33..$nday}
set -A days {1..$nday}

for day in ${days[@]} ; do
    day=`echo $day | awk '{printf "%03d", $1}'`
    ncks -a -O HYCOM_GLBa0.08_ssh_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_temp_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_salt_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_u_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    ncks -a -A HYCOM_GLBa0.08_v_${year}_${day}.nc HYCOM_GLBa0.08_${year}_${day}.nc
    rm -f HYCOM_GLBa0.08_ssh_${year}_${day}.nc HYCOM_GLBa0.08_temp_${year}_${day}.nc \ HYCOM_GLBa0.08_salt_${year}_${day}.nc HYCOM_GLBa0.08_u_${year}_${day}.nc \ 
HYCOM_GLBa0.08_v_${year}_${day}.nc
done

Gosh, should convert to bash.

Dear Kate

Thank you so much for sharing your works. I finished the 5th step and going to run make ic and bc file and got these error.

bc file

Code: Select all

Traceback (most recent call last):
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/make_bdry_file.py", line 57, in <module>
    year = int(sys.argv[1])
IndexError: list index out of range
>>> 

ic file (Ic problem has been solved).

Code: Select all

Traceback (most recent call last):
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/make_ic_file.py", line 30, in <module>
    remap(file, 'temp', src_grd, dst_grd, dst_dir=dst_dir)
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/remap.py", line 127, in remap
    dxy=dxy, cdepth=cdepth, kk=kk)
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/Grid_HYCOM/flood_fast.py", line 107, in flood_fast
    varz[bottom[j,i]:,j,i] = varz[bottom[j,i],j,i]
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices
>>> 

Thank you.
Regard.

[kate]

year = int(sys.argv[1])

Did you give it a year as a command-line argument? It is expecting one. There's also a commented out line:

Code: Select all

#lst_year = sys.argv[1:]

that can be enabled, allowing you to give it a list of years. Or you can edit the script to provide the year(s) internally.

For the other, I thought I had fixed all those "bottom" problems. Indeed, my copy of that line contains:

Code: Select all

              varz[int(bottom[j,i]):,j,i] = varz[int(bottom[j,i]),j,i]

[bangdt]

year = int(sys.argv[1])

Did you give it a year as a command-line argument? It is expecting one. There's also a commented out line:

Code: Select all

#lst_year = sys.argv[1:]

that can be enabled, allowing you to give it a list of years. Or you can edit the script to provide the year(s) internally.

For the other, I thought I had fixed all those "bottom" problems. Indeed, my copy of that line contains:

Code: Select all

              varz[int(bottom[j,i]):,j,i] = varz[int(bottom[j,i]),j,i]

Dear Kate
thank you for your answers, I put the year in then passed (when not python asked me for the year, I input but it still give me error)
the bc script run for a while then a message box appear to say that "Your program is still running! Do you want to kill it?" I said "No" for 4 time then the script seem to be stopped. May be I have some problems with these lines

Code: Select all

processes = 4
p = Pool(processes)
# Trick to pass more than one arg
partial_do_file = partial(do_file, src_grd=src_grd, dst_grd=dst_grd)
results = p.map(partial_do_file, lst_file)

I have been run ROMS for multiple cores for several times. I have no idea about this.

Regard.

[kate]

Each system is different. I run on one with front-end nodes containing 8 cores, while the compute nodes contain 24-28 cores. I run these scripts on the front-end nodes because they have access to the filesystems containing these files, but I only ask for 4 cores in my "pool" so other people can get things done too. The limits on these two types of jobs are different, and you might have run into some system limits on cputime, etc. Do you have a sys admin to talk to about where and how to run jobs? These things can take time and my current ROMS-to-ROMS boundary extraction is taking about two days per year.

[bangdt]

The limits on these two types of jobs are different, and you might have run into some system limits on cputime, etc. Do you have a sys admin to talk to about where and how to run jobs? .

Dear Kate
Thank you for help. I am trying to run all in my laptop before asking admin about my cputime. I run the script and have this message box



I notice that the cpu usage all are very low so may be the script was actually stopped.
please forgive me if I am a newbie in linux base systems so some time it makes me hard to understand the problems.
Have a nice week.
Regard

[kate]

Maybe you should try a pool size of 1, not running in parallel. I'm not familiar with that exact display, but the swap seems high, like it doesn't all fit in memory and none of them are making any progress.

[bangdt]

Maybe you should try a pool size of 1, not running in parallel. I'm not familiar with that exact display, but the swap seems high, like it doesn't all fit in memory and none of them are making any progress.

Dear Kate
Thanks for help. I put some flag into the code then

so I have error with p = Pool(processes), I also try with processes = 1 but the error stills. Could you please help me other way to skip using Pool?

Thank you.
Regard.

[kate]

Does it object to processes = 1 the same way as processes = 4, or was there a different error? What if you don't let it kill the job?

[bangdt]

Does it object to processes = 1 the same way as processes = 4, or was there a different error? What if you don't let it kill the job?

I have same err for 4 and 1 processes. If i dont kill it then the script hang like picture above with no err message and activity.

[kate]

Can you try on a system with more memory?

[bangdt]

Can you try on a system with more memory?

I restart and try to run it and the err stills remain. My system is 8gb of ram with nearly 6gb free. May be problem is related to system open mpi libraries can not linked to anaconda packages? I have succeced running ROMS test cases with mpi run.
I try to run with this

Code: Select all

print('start')
processes = 4
print('stop1')
p = Pool(processes)
# Trick to pass more than one arg
print('stop2')
partial_do_file = partial(do_file, src_grd=src_grd, dst_grd=dst_grd)
print('stop3')
results = p.map(partial_do_file, lst_file)
print('stop4')

when running with these stop1-4 flags, the message boxes appear 4 times, when stop1 printed I have 1st message box, I picked "No" 3 times then the stop2 printed out together with 4th message box, I continue pick "No" on 4th message box then the stop3 printed out then script hang and seem to be never printout stop4, at this time my system still have nearly 2gb ram free.

I changed the code from

Code: Select all

partial_do_file = partial(do_file, src_grd=src_grd, dst_grd=dst_grd)
results = p.map(partial_do_file, lst_file)

to

Code: Select all

results = p.map(partial(do_file(lst_file, src_grd=src_grd, dst_grd=dst_grd)))

then the script passed then I have this

Code: Select all

Traceback (most recent call last):
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/make_bdry_file.py", line 101, in <module>
    results = p.map(partial(do_file(file=lst_file, src_grd=src_grd, dst_grd=dst_grd)))
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/make_bdry_file.py", line 30, in do_file
    zeta = remap_bdry(file, 'ssh', src_grd, dst_grd, dst_dir=dst_dir)
  File "/work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/remap_bdry.py", line 28, in remap_bdry
    dst_file = src_file.rsplit('/')[-1]
AttributeError: 'list' object has no attribute 'rsplit'

Regard.

[kate]

You can read up on multiprocessing online. Something wasn't right about your last attempt.

I would check your system limits. With a bash shell:

Code: Select all

(snowdrifts) [kshedstrom@chinook02 Tidal_bay]$ ulimit -a
core file size          (blocks, -c) 0
data seg size           (kbytes, -d) unlimited
scheduling priority             (-e) 0
file size               (blocks, -f) unlimited
pending signals                 (-i) 257144
max locked memory       (kbytes, -l) unlimited
max memory size         (kbytes, -m) unlimited
open files                      (-n) 1024
pipe size            (512 bytes, -p) 8
POSIX message queues     (bytes, -q) 819200
real-time priority              (-r) 0
stack size              (kbytes, -s) 10240
cpu time               (seconds, -t) unlimited
max user processes              (-u) 1024
virtual memory          (kbytes, -v) unlimited
file locks                      (-x) unlimited

and with the csh family:

Code: Select all

chinook03.rcs.alaska.edu 318% limit
cputime      unlimited
filesize     unlimited
datasize     unlimited
stacksize    10240 kbytes
coredumpsize 0 kbytes
memoryuse    unlimited
vmemoryuse   unlimited
descriptors  1024 
memorylocked unlimited
maxproc      1024 

Honestly, what we really want is to have this using xarray and dask instead of multiprocessing. I don't (yet) know enough to do that rewrite.

[bangdt]

You can read up on multiprocessing online. Something wasn't right about your last attempt.

I would check your system limits. With a bash shell:

My information

Code: Select all

(pr2) bangdt@C2:~$ ulimit -a
core file size          (blocks, -c) 0
data seg size           (kbytes, -d) unlimited
scheduling priority             (-e) 0
file size               (blocks, -f) unlimited
pending signals                 (-i) 31039
max locked memory       (kbytes, -l) 16384
max memory size         (kbytes, -m) unlimited
open files                      (-n) 1024
pipe size            (512 bytes, -p) 8
POSIX message queues     (bytes, -q) 819200
real-time priority              (-r) 0
stack size              (kbytes, -s) 8192
cpu time               (seconds, -t) unlimited
max user processes              (-u) 31039
virtual memory          (kbytes, -v) unlimited
file locks                      (-x) unlimited

and

Code: Select all

# limit
cputime 	unlimited
filesize 	unlimited
datasize 	unlimited
stacksize 	8192 kbytes
coredumpsize 	0 kbytes
memoryuse 	unlimited
memorylocked 	16384 kbytes
maxproc 	31039 
openfiles 	1024 

[kate]

Maybe this is biting you:

Code: Select all

memorylocked    16384 kbytes

[bangdt]

Maybe this is biting you:

Code: Select all

memorylocked    16384 kbytes

Dear Kate

by running this code

Code: Select all

(pr2) root@C2:~# ulimit -l
16384
(pr2) root@C2:~# ulimit -l unlimited
(pr2) root@C2:~# ulimit -l
unlimited

then my system has

Code: Select all

(pr2) root@C2:~# ulimit -a
core file size          (blocks, -c) 0
data seg size           (kbytes, -d) unlimited
scheduling priority             (-e) 0
file size               (blocks, -f) unlimited
pending signals                 (-i) 31040
max locked memory       (kbytes, -l) unlimited
max memory size         (kbytes, -m) unlimited
open files                      (-n) 1024
pipe size            (512 bytes, -p) 8
POSIX message queues     (bytes, -q) 819200
real-time priority              (-r) 0
stack size              (kbytes, -s) 8192
cpu time               (seconds, -t) unlimited
max user processes              (-u) 31040
virtual memory          (kbytes, -v) unlimited
file locks                      (-x) unlimited

and

Code: Select all

(pr2) root@C2:~# csh
# limit
cputime 	unlimited
filesize 	unlimited
datasize 	unlimited
stacksize 	8192 kbytes
coredumpsize 	0 kbytes
memoryuse 	unlimited
memorylocked 	unlimited
maxproc 	31040 
openfiles 	1024 

I assume that by unlock the memories size I could load the full lst_file in to the memory and run,
but when I try to run the script the problem stills remain.
may be stuck at this line and can not load the lst_file to call do_file.

Code: Select all

results = p.map(partial_do_file, lst_file)

Regard

[kate]

What you want is to call do_file for each file in the list. Did it print the list of files to operate on? Is it what you expected? Then you can comment out the chunk with processes = 4 to the end. Add a loop going through the list of files, calling do_file(file, src_grd, dst_grd), something like:

Code: Select all

for file in lst_file:
   do_file(file, src_grd, dst_grd)

[bangdt]

What you want is to call do_file for each file in the list. Did it print the list of files to operate on? Is it what you expected? Then you can comment out the chunk with processes = 4 to the end. Add a loop going through the list of files, calling do_file(file, src_grd, dst_grd), something like:

Code: Select all

for file in lst_file:
   do_file(file, src_grd, dst_grd)

Dear Kate
Thank you so much for your suggestion. Now I have all ic and bc created. But I still have problem with them.

1- All my ic and bc file are 2d grid file, not 3d comparing to the source data
2- My bc file doesnt have temp, salt, u and v data, only have 2d grid data in mask on v and u, and all the mask on u, v are the same value (=1), the bathymetry data in the bc file are ok.

my grid id is

Code: Select all

# gridid definition file
#
id      = EASTSEA
name    = EASTSEA
grdfile = /work/apps/anaconda3/envs/pr2/lib/python2.7/site-packages/pyroms_toolbox/works/East-Sea_HYCOM/East_sea_400.nc
N       = 20
grdtype = roms
Vtrans  = 1
theta_s = 5
theta_b = 0.4
Tcline  = 3

I am not sure if my ic and bc file are ok or how to check it (I use IDV to view them)
https://drive.google.com/drive/folders/ ... sp=sharing
Regard.

[kate]

I like to use ncview to get a quick and dirty view of what's in a netCDF file. Looking at 'ncdump -h' on a file can also be useful. If you didn't get temp and salt, then something didn't go right in the boundary extraction.

If you like the idea of using these Python scripts, perhaps it's time to learn what you can about Python/numpy and the rest. Check out pdb, the Python debugger.

[bangdt]

Dear Kate

It is long time no see, I hope you are doing well.
thank you so much for help us go through this forum smoothly.
several months ago I have install pyroms 2 and successfully made boundary and initial file. Today after install pyroms 3 (python 3) it return an error.

when I run this code

Code: Select all

import matplotlib
matplotlib.use('Agg')

import numpy as np
import netCDF4
from datetime import datetime
import pyroms
import pyroms_toolbox
#from pyroms_toolbox import date2jday
import sys



def create_HYCOM_file(name, time, lon, lat, var):

    #create netCDF file
    nc = netCDF4.Dataset(name, 'w', format='NETCDF3_64BIT')
    nc.Author = sys._getframe().f_code.co_name
    nc.Created = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
    nc.title = 'HYCOM + NCODA Global 1/12 Analysis (GLBa0.08)'

    #create dimensions
    Mp, Lp = lon.shape
    nc.createDimension('lon', Lp)
    nc.createDimension('lat', Mp)
    nc.createDimension('ocean_time', None)

    #create variables        
    nc.createVariable('lon', 'f', ('lat', 'lon'))
    nc.variables['lon'].long_name = 'longitude'
    nc.variables['lon'].units = 'degrees_east'
    nc.variables['lon'][:] = lon

    nc.createVariable('lat', 'f', ('lat', 'lon'))
    nc.variables['lat'].long_name = 'latitude'
    nc.variables['lat'].units = 'degrees_north'
    nc.variables['lat'][:] = lat

    nc.createVariable('ocean_time', 'f', ('ocean_time'))
    nc.variables['ocean_time'].units = 'days since 1900-01-01 00:00:00'
    nc.variables['ocean_time'].calendar = 'LEAP'
    nc.variables['ocean_time'][0] = time

    nc.createVariable(outvarname, 'f', ('ocean_time', 'lat', 'lon'), fill_value=spval)
    nc.variables[outvarname].long_name = long_name
    nc.variables[outvarname].units = units
    nc.variables[outvarname].coordinates = 'lon lat'
    nc.variables[outvarname][0] = var

    nc.close()

    print('Done with file %s' %name)




# get HYCOM Northeast Pacific data from 2007 to 2011

year = 2015
retry='True'

invarname = 'ssh'
outvarname = 'ssh'

#read grid and variable attributes from the first file
url='http://tds.hycom.org/thredds/dodsC/datasets/GLBa0.08/expt_91.1/2015/2d/archv.2015_001_00_2d.nc'
dataset = netCDF4.Dataset(url)
lon = dataset.variables['Longitude'][1550-9:1950,300:700]
lat = dataset.variables['Latitude'][1550-9:1950,300:700]
#spval = dataset.variables[invarname]._FillValue
units = dataset.variables[invarname].units
long_name = dataset.variables[invarname].long_name
dataset.close()


retry_day = []

# loop over daily files
if year%4 == 0:
    daysinyear = 366
else:
#    daysinyear = 365
    daysinyear = 32
for day in range(1,daysinyear+1):
    print('Processing file for day %03d, year %04d' %(day, year))
    url='http://tds.hycom.org/thredds/dodsC/datasets/GLBa0.08/expt_91.1/2015/2d/archv.%04d_%03d_00_2d.nc' %(year,day)
    #get data from server
    try:
        dataset = netCDF4.Dataset(url)
        var = dataset.variables[invarname][0,1550-9:1950,300:700]
        spval = var.get_fill_value()
        dataset.close()
    except:
        print('No file on the server... We skip this day.')
        retry_day.append(day)
        continue

    #create netCDF file
    outfile = 'data/HYCOM_GLBa0.08_%s_%04d_%03d.nc' %(outvarname,year,day)
    jday = pyroms_toolbox.date2jday(datetime(year, 1, 1)) + day - 1
    create_HYCOM_file(outfile, jday, lon, lat, var)


if retry == 'True':
    if len(retry_day) != 0:
        print("Some files have not been downloded... Let's try again")
    while len(retry_day) != 0:
        for day in retry_day:
            print('Retry file for day %03d, year %04d' %(day, year))
            url='http://tds.hycom.org/thredds/dodsC/datasets/GLBa0.08/expt_91.1/2015/2d/archv.%04d_%03d_00_2d.nc' %(year,day)
            #get data from server
            try:
                dataset = netCDF4.Dataset(url)
                var = dataset.variables[invarname][0,1550-9:1950,300:700]
                spval = var.get_fill_value()
                dataset.close()
            except:
                print('No file on the server... We skip this day.')
            continue

            #create netCDF file
            outfile = 'data/HYCOM_GLBa0.08_%s_%04d_%03d.nc' %(outvarname,year,day)
            jday = pyroms_toolbox.date2jday(datetime(year, 1, 1)) + day - 1
            create_HYCOM_file(outfile, jday, lon, lat, var)

            retry_day.remove(day)

it returned this error

Code: Select all

>>> 
= RESTART: /home/bangdt/hymenet/pyroms/pyroms_toolbox/pyroms_toolbox/Grid_HYCOM/1-get_hycom_GLBa0.08_ssh_2015.py
Processing file for day 001, year 2015
Traceback (most recent call last):
  File "/home/bangdt/hymenet/pyroms/pyroms_toolbox/pyroms_toolbox/Grid_HYCOM/1-get_hycom_GLBa0.08_ssh_2015.py", line 100, in <module>
    jday = pyroms_toolbox.date2jday(datetime(year, 1, 1)) + day - 1
TypeError: 'module' object is not callable
>>> 

Best
Bang

[Dan_chan]

Hi,

I test pyroms example Palau_HYCOM directly and make grid nc file. But when I run 'python make_remap_weight_file.py' or 'python make_bdry_file.py 2015', it return an error on spheical. I find out the error command in .py script is dstgrd = pyroms.grid.get_ROMS_grid('PALAU1')

Code: Select all

>>> dstgrd = pyroms.grid.get_ROMS_grid('PALAU1')

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/dan_chen/ROMS/pyroms/pyroms/pyroms/grid.py", line 487, in get_ROMS_grid
    hgrd = get_ROMS_hgrid(gridid)
  File "/Users/dan_chen/ROMS/pyroms/pyroms/pyroms/grid.py", line 228, in get_ROMS_hgrid
    spherical = nc.variables['spherical'][0]
KeyError: 'spherical'

I don't find spherical variable in my grid.nc, so should I add this var into grid.nc ?

Thanks for any help in advance.

[andrearuju]

Hi Dan,
I got the same error. Did you figure out how to solve it?
Apparently the nc file created by get_hycom_GLBa0.08_Palau_grid.py has no spherical variable. Even if I comment something in the code and put spherical = 'T' in the grid.py I get the following error:
NetCDF file must contain lon_vert and lat_vert or lon_rho, lat_rho, lon_psi, lat_psi for a geographical grid.
That is reasonable since the grid nc file has lon lat as variables, not lon_vert, lat_vert
Any suggestion?
Andrea

[kate]

There are two options at least: (a) change pyroms to accept the lat, lon variables as they are in the file or (b) change the file to include the fields that pyroms is looking for.