Index: Scripts/e2o-getfromwci.py
===================================================================
diff -u -rfbd077ecd824c16e9e207b542e5ee2f9f9ff92dd -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- Scripts/e2o-getfromwci.py	(.../e2o-getfromwci.py)	(revision fbd077ecd824c16e9e207b542e5ee2f9f9ff92dd)
+++ Scripts/e2o-getfromwci.py	(.../e2o-getfromwci.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -7,6 +7,7 @@
 
 # Uncomment these if needed
 from numpy import *
+import numpy as np
 from scipy import *
 from matplotlib import *
 from pylab import *
@@ -136,7 +137,7 @@
         if len(tot) == 0:
             tot = p_mean.copy()
         else:
-            tot = hstack((tot, p_mean))
+            tot = np.hstack((tot, p_mean))
 
         arcnt = 0
         for a in timeObj:
@@ -152,7 +153,7 @@
                 "PCRaster",
                 lon,
                 lat[::-1],
-                flipud(p_select[arcnt, :, :]),
+                np.flipud(p_select[arcnt, :, :]),
                 -999.0,
             )
             arcnt = arcnt + 1
Index: Scripts/pcr2netcdf.py
===================================================================
diff -u -rfbd077ecd824c16e9e207b542e5ee2f9f9ff92dd -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- Scripts/pcr2netcdf.py	(.../pcr2netcdf.py)	(revision fbd077ecd824c16e9e207b542e5ee2f9f9ff92dd)
+++ Scripts/pcr2netcdf.py	(.../pcr2netcdf.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -70,6 +70,7 @@
 import getopt
 import sys
 from numpy import *
+import numpy as np
 import netCDF4 as nc4
 import cftime
 import osgeo.gdal as gdal
@@ -126,7 +127,7 @@
     # get rasterband entry
     TempBand = TempDataset.GetRasterBand(1)
     # fill rasterband with array
-    TempBand.WriteArray(data.astype(float32), 0, 0)
+    TempBand.WriteArray(data.astype(np.float32), 0, 0)
     TempBand.FlushCache()
     TempBand.SetNoDataValue(FillVal)
 
@@ -357,7 +358,7 @@
         latlen = len(nc_trg.variables["lat"])
         lonlen = len(nc_trg.variables["lon"])
 
-    timestepbuffer = zeros((bufsize, latlen, lonlen))
+    timestepbuffer = np.zeros((bufsize, latlen, lonlen))
 
     # now loop over all time steps, check the date and write valid dates to a list, write time series to PCRaster maps
     for nn, curTime in enumerate(timeList):
Index: Scripts/wflow_sbm_rtc.py
===================================================================
diff -u -rfbd077ecd824c16e9e207b542e5ee2f9f9ff92dd -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- Scripts/wflow_sbm_rtc.py	(.../wflow_sbm_rtc.py)	(revision fbd077ecd824c16e9e207b542e5ee2f9f9ff92dd)
+++ Scripts/wflow_sbm_rtc.py	(.../wflow_sbm_rtc.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -11,6 +11,7 @@
 
 from wflow.wf_DynamicFramework import *
 
+import pcraster as pcr
 import wflow.wflow_bmi as bmi
 import logging
 import wflow.wflow_adapt as adapter
@@ -181,22 +182,22 @@
 
 
 # In[]:  Read and map reservoir inflow and outflow locations
-Reservoir_inflow = pcr2numpy(
-    scalar(pcraster.readmap(os.path.abspath(inflow_map))), np.NaN
+Reservoir_inflow = pcr.pcr2numpy(
+    pcr.scalar(pcr.readmap(os.path.abspath(inflow_map))), np.NaN
 )
-Reservoir_outflow = pcr2numpy(
-    scalar(pcraster.readmap(os.path.abspath(outflow_map))), np.NaN
+Reservoir_outflow = pcr.pcr2numpy(
+    pcr.scalar(pcr.readmap(os.path.abspath(outflow_map))), np.NaN
 )
 
 inflow_list = list(np.unique(Reservoir_inflow[~np.isnan(Reservoir_inflow)]))
 outflow_list = list(np.unique(Reservoir_outflow[~np.isnan(Reservoir_outflow)]))
 
 
 # In[]:  Overwrite TopoLdd with modified version
-ldd = pcraster.pcr2numpy(
-    pcraster.readmap(os.path.join(dir_wflow, ldd_map)), np.NaN
+ldd = pcr.pcr2numpy(
+    pcr.readmap(os.path.join(dir_wflow, ldd_map)), np.NaN
 ).astype(np.float32)
-LA_model.set_value("TopoLdd", flipud(ldd).copy())
+LA_model.set_value("TopoLdd", np.flipud(ldd).copy())
 
 
 ########################################################################
@@ -212,15 +213,15 @@
         toread = gettimestepfname(
             thisstack, os.path.join(dir_wflow, "inmaps"), timecounter + 1
         )
-        inmstackbuf[thisstack] = flipud(
-            pcr2numpy(scalar(pcraster.readmap(os.path.abspath(toread))), -999.0)
+        inmstackbuf[thisstack] = np.flipud(
+            pcr.pcr2numpy(pcr.scalar(pcr.readmap(os.path.abspath(toread))), -999.0)
         ).copy()
         LA_model.set_value(thisstack, inmstackbuf[thisstack])
 
     print("calculation timestep = " + str(timecounter))
 
     # Get the inflow from the wflow model runoff map and map
-    inflowQ = flipud(LA_model.get_value("SurfaceRunoff")).copy()
+    inflowQ = np.flipud(LA_model.get_value("SurfaceRunoff")).copy()
 
     # Map the sum of WFlow Inflow to RTC
     for idx, wflow_id in enumerate(id_in_wflow):
@@ -239,10 +240,10 @@
         rtc_id = id_out_rtc[id_out_wflow.index(str(wflow_id))]
         Qout = RTC_model.get_var(rtc_id)
         Qout = 300.0
-        if isfinite(Qout):  # no nan's into wflow
+        if np.isfinite(Qout):  # no nan's into wflow
             inflowfield[Reservoir_outflow == int(wflow_id)] += Qout
 
-    LA_model.set_value("IF", flipud(inflowfield).copy())
+    LA_model.set_value("IF", np.flipud(inflowfield).copy())
     # This not not bmi but needed to update the kinematic wave reservoir
     LA_model.update()
     t += LA_dt
Index: doc/pcraster-for-doc-only/pcraster/framework/aggregationfunctions.py
===================================================================
diff -u -r9d329e51cdc321edb86cc2518eee59bbc4f94159 -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- doc/pcraster-for-doc-only/pcraster/framework/aggregationfunctions.py	(.../aggregationfunctions.py)	(revision 9d329e51cdc321edb86cc2518eee59bbc4f94159)
+++ doc/pcraster-for-doc-only/pcraster/framework/aggregationfunctions.py	(.../aggregationfunctions.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -5,6 +5,7 @@
 import string
 import numpy
 import numpy.ma
+import numpy as np
 from pcraster import *
 from frameworkBase import generateNameS, generateNameT, generateNameST
 import generalfunctions, regression
Index: doc/pcraster-for-doc-only/pcraster/framework/regression.py
===================================================================
diff -u -r9d329e51cdc321edb86cc2518eee59bbc4f94159 -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- doc/pcraster-for-doc-only/pcraster/framework/regression.py	(.../regression.py)	(revision 9d329e51cdc321edb86cc2518eee59bbc4f94159)
+++ doc/pcraster-for-doc-only/pcraster/framework/regression.py	(.../regression.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -3,6 +3,7 @@
 
 
 from numpy import *
+import numpy as np
 import sys
 
 
Index: wflow/wflow_delwaq.py
===================================================================
diff -u -rc427e9b948a7fbad40a1c828b68355cb277dc5e0 -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- wflow/wflow_delwaq.py	(.../wflow_delwaq.py)	(revision c427e9b948a7fbad40a1c828b68355cb277dc5e0)
+++ wflow/wflow_delwaq.py	(.../wflow_delwaq.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -173,30 +173,30 @@
     nr_inflowtypes = len(inflowtypes)
 
     # for i in range(nr_inflowtypes-1):
-    #    totareas = vstack((totareas,areas))
+    #    totareas = np.vstack((totareas,areas))
     totareas = areas
 
     arid = 0
     for i in range(len(pointer)):
         if pointer[i, 1] < 0:
             print(
                 "'BD_"
-                + str(absolute(pointer[i, 1]))
+                + str(np.absolute(pointer[i, 1]))
                 + "'  '"
-                + str(absolute(pointer[i, 1]))
+                + str(np.absolute(pointer[i, 1]))
                 + "'"
                 + " 'Outflow'",
                 file=exfile,
             )
         elif pointer[i, 0] < 0:
-            # ar = int(absolute(totareas[arid]))
+            # ar = int(np.absolute(totareas[arid]))
             ar = totareas[arid]
             print(
                 "'BD_"
-                + str(absolute(pointer[i, 0]))
+                + str(np.absolute(pointer[i, 0]))
                 + "' "
                 + "'"
-                + str(absolute(pointer[i, 0]))
+                + str(np.absolute(pointer[i, 0]))
                 + "'"
                 + " 'Area_"
                 + str(ar)
@@ -244,10 +244,10 @@
     # First convert the array to a 32 bit float
     totareas = datablock
     for i in range(boundids - 1):
-        totareas = vstack((totareas, datablock))
+        totareas = np.vstack((totareas, datablock))
 
-    artow = np.array(totareas, dtype=float32).copy()
-    timear = np.array(ttime, dtype=int32)
+    artow = np.array(totareas, dtype=np.float32).copy()
+    timear = np.array(ttime, dtype=np.int32)
     if os.path.isfile(fname):  # append to existing file
         fp = open(fname, "ab")
         tstr = timear.tostring() + artow.tostring()
@@ -327,79 +327,79 @@
 
     # remove all non-active cells
     np_catchid = np_catchid[np_catchid > 0.0]
-    np_upbound = np_upbound[isfinite(np_upbound)]
-    np_flowto = np_flowto[isfinite(np_flowto)]
-    np_ptid = np_ptid[isfinite(np_ptid)]
+    np_upbound = np_upbound[np.isfinite(np_upbound)]
+    np_flowto = np_flowto[np.isfinite(np_flowto)]
+    np_ptid = np_ptid[np.isfinite(np_ptid)]
     np_flowto = np_flowto.reshape(len(np_flowto), 1)
     np_ptid = np_ptid.reshape(len(np_ptid), 1)
     np_catchid = np_catchid.reshape(len(np_catchid), 1)
     # Now make catchid a list
     np_catchid = np_catchid.flatten()
-    np_catchid = np.array(int_(np_catchid), dtype="|S").tolist()
+    np_catchid = np.array(np.int_(np_catchid), dtype="|S").tolist()
     # find all downstream segments (flowto == ptid)
     # now set the flowto points (outflows, usually just one) also to negative
-    lowerck = absolute(np_ptid) == absolute(np_flowto)
+    lowerck = np.absolute(np_ptid) == np.absolute(np_flowto)
     # mak epointer matrix and add to zero zolumns
-    orgpointer = hstack(
-        (np_ptid, np_flowto, zeros((len(np_flowto), 1)), zeros((len(np_flowto), 1)))
+    orgpointer = np.hstack(
+        (np_ptid, np_flowto, np.zeros((len(np_flowto), 1)), np.zeros((len(np_flowto), 1)))
     )
     pointer = orgpointer.copy()
     # Pointer labels:
     #    negative: outflow boundary
     #    zero    : internal flow
     #    positive: inflow boundary
-    pointer_labels = zeros((len(np_flowto)), dtype=numpy.int)
+    pointer_labels = np.zeros((len(np_flowto)), dtype=numpy.int)
     extraboun = []
     # Add the inflow boundaries here.
     cells = pointer[:, 0]
     cells = cells.reshape(len(cells), 1)
     bounid = cells.copy()
-    zzerocol = zeros((len(np_flowto), 1), dtype=numpy.int)
+    zzerocol = np.zeros((len(np_flowto), 1), dtype=numpy.int)
 
     # outflow to pointer
     # point -> - point
     lopt = np_ptid[lowerck]
     lopt = lopt.reshape(len(lopt), 1)
-    zerocol = zeros((len(lopt), 1))
+    zerocol = np.zeros((len(lopt), 1))
     lowerids = np.arange(1, len(lopt) + 1) * -1
-    # of = hstack((lopt,lopt * -1.0,zerocol,zerocol))
+    # of = np.hstack((lopt,lopt * -1.0,zerocol,zerocol))
     lowerids = lowerids.reshape(len(lowerids), 1)
-    of = hstack((lopt, lowerids, zerocol, zerocol))
+    of = np.hstack((lopt, lowerids, zerocol, zerocol))
 
     # Now remove double pointer to itself and replace by lower boundary
     lowerck = pointer[:, 0] == pointer[:, 1]
     pointer[lowerck, :] = of
     pointer_labels[lowerck] = -1
-    start = absolute(lowerids.min()) + 1
+    start = np.absolute(lowerids.min()) + 1
     bouns = 1
     for idd in range(1, difboun + 1):
         bounid = np.arange(start, (len(cells) + start)).reshape((len(cells), 1)) * -1.0
         if bouns == 1:
-            extraboun = hstack((bounid, cells, zzerocol, zzerocol))
+            extraboun = np.hstack((bounid, cells, zzerocol, zzerocol))
         else:
-            extraboun = vstack((extraboun, hstack((bounid, cells, zzerocol, zzerocol))))
-        pointer_labels = hstack((pointer_labels, zzerocol[:, 0] + bouns))
+            extraboun = np.vstack((extraboun, np.hstack((bounid, cells, zzerocol, zzerocol))))
+        pointer_labels = np.hstack((pointer_labels, zzerocol[:, 0] + bouns))
         bouns = bouns + 1
         start = start + len(cells)
 
     res = []
     for idd in range(1, difboun + 1):
         ct = list(np_catchid)
         print("ct: ")
-        print(unique(ct))
+        print(np.unique(ct))
         for i in range(0, len(np_catchid)):
             ct[i] = np_catchid[i] + "_" + str(idd)
         res = res + ct
-    print(unique(res))
+    print(np.unique(res))
     np_catchid = res
-    # pointer = vstack((pointer,extraboun))
+    # pointer = np.vstack((pointer,extraboun))
     # now catchment id's
-    # zerocol = zeros((len(np_catchid),1))
-    # extraboun= hstack((np_catchid,cells,zerocol,zerocol))
+    # zerocol = np.zeros((len(np_catchid),1))
+    # extraboun= np.hstack((np_catchid,cells,zerocol,zerocol))
     # print np_catchid
 
     if len(extraboun) > 0:
-        pointer = vstack((pointer, extraboun))
+        pointer = np.vstack((pointer, extraboun))
 
     return ptid, pointer, pointer_labels, np_ptid.flatten(), np_catchid
 
@@ -410,7 +410,7 @@
     missing values are removed. Used for generating delwaq data
     """
     ttar = pcr.pcr2numpy(pcrmap, np.nan).flatten()
-    ttar = ttar[isfinite(ttar)]
+    ttar = ttar[np.isfinite(ttar)]
 
     return ttar
 
@@ -568,13 +568,13 @@
     np_gauges = pcr.pcr2numpy(gauges, np.nan).flatten()
     np_segs = pcr.pcr2numpy(segs, np.nan).flatten()
 
-    np_gauges = np_gauges[isfinite(np_gauges)]
-    np_segs = np_segs[isfinite(np_segs)]
+    np_gauges = np_gauges[np.isfinite(np_gauges)]
+    np_segs = np_segs[np.isfinite(np_segs)]
 
     if len(np_segs) != len(np_gauges):
         logger.error("Gauges and segments do not match!")
 
-    pts = size(np_segs)
+    pts = np.size(np_segs)
     exfile = open(fname, "w")
     print("%d ; nr of locations" % pts, file=exfile)
     print("; 'outlocname' numberofsegments segment list", file=exfile)
@@ -1363,7 +1363,7 @@
 
     thecells = pcr.pcr2numpy(modelmap, np.nan).flatten()
     nrcells = len(thecells)
-    nractcells = len(thecells[isfinite(thecells)])
+    nractcells = len(thecells[np.isfinite(thecells)])
 
     logger.info("Total number gridcells (including inactive): " + str(nrcells))
     logger.info("Total number of used gridcells: " + str(nractcells))
@@ -1374,7 +1374,7 @@
     upar = pcr.pcr2numpy(pcr.scalar(upbound), np.nan).flatten()
     logger.info(
         "Number of upstream cells (without upstream connection): "
-        + str(len(upar[isfinite(upar)]))
+        + str(len(upar[np.isfinite(upar)]))
     )
     pcr.report(upbound, dwdir + "/debug/upbound.map")
 
@@ -1425,15 +1425,15 @@
     internalflowlength = pcr.readmap(caseId + "/" + runId + "/outsum/DCL.map")
     surface_map = internalflowwidth * internalflowlength
     surface_block = dw_pcrToDataBlock(surface_map)
-    logger.info("Writing surface.dat. Nr of points: " + str(size(surface_block)))
+    logger.info("Writing surface.dat. Nr of points: " + str(np.size(surface_block)))
     dw_WriteSegmentOrExchangeData(
         0, dwdir + "/includes_flow/surface.dat", surface_block, 1, WriteAscii
     )
 
     # create dummy length file
-    length_block = zeros(pointer.shape[0] * 2) + 0.5
+    length_block = np.zeros(pointer.shape[0] * 2) + 0.5
     # write  length file
-    logger.info("Writing length.dat. Nr of points: " + str(size(length_block)))
+    logger.info("Writing length.dat. Nr of points: " + str(np.size(length_block)))
     dw_WriteSegmentOrExchangeData(
         0, dwdir + "/includes_flow/length.dat", length_block, 1, WriteAscii
     )
@@ -1474,7 +1474,7 @@
 
             # volume for each timestep and number of segments
 
-            logger.info("Writing volumes.dat. Nr of points: " + str(size(volume_block)))
+            logger.info("Writing volumes.dat. Nr of points: " + str(np.size(volume_block)))
             dw_WriteSegmentOrExchangeData(
                 i, dwdir + "/includes_flow/volume.dat", volume_block, 1, WriteAscii
             )
@@ -1498,14 +1498,14 @@
                 thesource = read_timestep(nc, source, ts, logger, caseId, runId)
                 thesource = zero_map + thesource
                 flow_block_IN = dw_pcrToDataBlock(thesource)
-                flowblock = hstack((flowblock, flow_block_IN))
-                area_block = hstack((area_block, surface_block))
+                flowblock = np.hstack((flowblock, flow_block_IN))
+                area_block = np.hstack((area_block, surface_block))
 
-            logger.info("Writing flow.dat. Nr of points: " + str(size(flowblock)))
+            logger.info("Writing flow.dat. Nr of points: " + str(np.size(flowblock)))
             dw_WriteSegmentOrExchangeData(
                 i, dwdir + "/includes_flow/flow.dat", flowblock, 1, WriteAscii
             )
-            logger.info("Writing area.dat. Nr of points: " + str(size(area_block)))
+            logger.info("Writing area.dat. Nr of points: " + str(np.size(area_block)))
             dw_WriteSegmentOrExchangeData(
                 i, dwdir + "/includes_flow/area.dat", area_block, 1, WriteAscii
             )
@@ -1528,22 +1528,22 @@
         i = i + timestepsecs
         logger.info("Writing last step..")
 
-        logger.info("Writing area.dat. Nr of points: " + str(size(area_block)))
+        logger.info("Writing area.dat. Nr of points: " + str(np.size(area_block)))
         dw_WriteSegmentOrExchangeData(
             i, dwdir + "/includes_flow/area.dat", area_block, 1, WriteAscii
         )
 
-        # logger.info("Writing surface.dat. Nr of points: " + str(size(surface_block)))
+        # logger.info("Writing surface.dat. Nr of points: " + str(np.size(surface_block)))
         # dw_WriteSegmentOrExchangeData(i,dwdir + '/includes_flow/surface.dat',surface_block,1,WriteAscii)
 
-        logger.info("Writing flow.dat. Nr of points: " + str(size(flowblock)))
+        logger.info("Writing flow.dat. Nr of points: " + str(np.size(flowblock)))
         dw_WriteSegmentOrExchangeData(
             i, dwdir + "/includes_flow/flow.dat", flowblock, 1, WriteAscii
         )
 
         volume_map = read_timestep(nc, "voln", ts, logger, caseId, runId)
         volume_block = dw_pcrToDataBlock(volume_map)
-        logger.info("Writing volumes.dat. Nr of points: " + str(size(volume_block)))
+        logger.info("Writing volumes.dat. Nr of points: " + str(np.size(volume_block)))
         dw_WriteSegmentOrExchangeData(
             i, dwdir + "/includes_flow/volume.dat", volume_block, 1, WriteAscii
         )
Index: wflow/wflow_hbv.py
===================================================================
diff -u -ra9498adee6baab0a0abaa331041be8948510167b -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- wflow/wflow_hbv.py	(.../wflow_hbv.py)	(revision a9498adee6baab0a0abaa331041be8948510167b)
+++ wflow/wflow_hbv.py	(.../wflow_hbv.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -584,10 +584,10 @@
             self.sh = {}
             res_ids = pcr.ifthen(self.ResStorFunc == 2, self.ReserVoirComplexLocs)
             np_res_ids = pcr.pcr2numpy(res_ids, 0)
-            np_res_ids_u = np.unique(np_res_ids[nonzero(np_res_ids)])
+            np_res_ids_u = np.unique(np_res_ids[np.nonzero(np_res_ids)])
             if np.size(np_res_ids_u) > 0:
-                for item in nditer(np_res_ids_u):
-                    self.sh[int(item)] = loadtxt(
+                for item in np.nditer(np_res_ids_u):
+                    self.sh[int(item)] = np.loadtxt(
                         self.Dir
                         + "/"
                         + self.intbl
@@ -598,10 +598,10 @@
             self.hq = {}
             res_ids = pcr.ifthen(self.ResOutflowFunc == 1, self.ReserVoirComplexLocs)
             np_res_ids = pcr.pcr2numpy(res_ids, 0)
-            np_res_ids_u = np.unique(np_res_ids[nonzero(np_res_ids)])
-            if size(np_res_ids_u) > 0:
-                for item in nditer(np_res_ids_u):
-                    self.hq[int(item)] = loadtxt(
+            np_res_ids_u = np.unique(np_res_ids[np.nonzero(np_res_ids)])
+            if np.size(np_res_ids_u) > 0:
+                for item in np.nditer(np_res_ids_u):
+                    self.hq[int(item)] = np.loadtxt(
                         self.Dir
                         + "/"
                         + self.intbl
Index: wflow/wflow_sbm.py
===================================================================
diff -u -rd7827e63ba2a822d7107471121efec49ad5a36d4 -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4
--- wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision d7827e63ba2a822d7107471121efec49ad5a36d4)
+++ wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4)
@@ -1284,10 +1284,10 @@
             self.sh = {}
             res_ids = pcr.ifthen(self.ResStorFunc == 2, self.ReserVoirComplexLocs)
             np_res_ids = pcr.pcr2numpy(res_ids, 0)
-            np_res_ids_u = np.unique(np_res_ids[nonzero(np_res_ids)])
+            np_res_ids_u = np.unique(np_res_ids[np.nonzero(np_res_ids)])
             if np.size(np_res_ids_u) > 0:
-                for item in nditer(np_res_ids_u):
-                    self.sh[int(item)] = loadtxt(
+                for item in np.nditer(np_res_ids_u):
+                    self.sh[int(item)] = np.loadtxt(
                         self.Dir
                         + "/"
                         + self.intbl
@@ -1298,10 +1298,10 @@
             self.hq = {}
             res_ids = pcr.ifthen(self.ResOutflowFunc == 1, self.ReserVoirComplexLocs)
             np_res_ids = pcr.pcr2numpy(res_ids, 0)
-            np_res_ids_u = np.unique(np_res_ids[nonzero(np_res_ids)])
-            if size(np_res_ids_u) > 0:
-                for item in nditer(np_res_ids_u):
-                    self.hq[int(item)] = loadtxt(
+            np_res_ids_u = np.unique(np_res_ids[np.nonzero(np_res_ids)])
+            if np.size(np_res_ids_u) > 0:
+                for item in np.nditer(np_res_ids_u):
+                    self.hq[int(item)] = np.loadtxt(
                         self.Dir
                         + "/"
                         + self.intbl