# ---------- Qgridder utils tutorial ---------- # A. Pryet - alexandre.pryet@ensegid.fr # ---------------------------------------------- # ------------------------------------------------------------ # Case 1 : Steady-state simulation with Modflow through Flopy # ------------------------------------------------------------ # ------------ Preprocessing in Qgis 2.0 ------------ # Build a grid with the Qgridder plugin # Add and fill "IBOUND" and "hstart" columns to the grid layer. # Option : prepare a "wells" layer with points # Option : prepare a "drains" layer with polylines # Option : prepare a "ghb" layer with polylines # ************************************************************ # FOLLOWING COMMANDS MUST BE EXECUTED FROM QGIS PYTHON CONSOLE # ************************************************************ # -------- Load python packages and set path ------------ import numpy as np import os import sys import pickle from matplotlib import pyplot as plt # path to flopy source code home = os.path.expanduser('~') sys.path.append(home+'/Programmes/flopy/svn/trunk/') # add path to Qgridder source code sys.path.append(home+'/Programmes/PyQgis/Qgridder/') # load flopy and qgridder from flopy.modflow import * from flopy.utils import * import qgridder_utils import ftools_utils # set current directory os.chdir(home + '/Programmes/flopy/tutorial') # ------------------------------------------------------------ # Pre-processing # ------------------------------------------------------------ # ------------ Grid data ------------ # select grid layer gridLayer = ftools_utils.getVectorLayerByName('grid') # Number of rows (along height) and columns (along width) nrow, ncol = qgridder_utils.get_rgrid_nrow_ncol(gridLayer) # delc, height or rows (along columns) and delrow width of cells of columns (along rows) delr, delc = qgridder_utils.get_rgrid_delr_delc(gridLayer) # Set ibound and hstart 2D array ibound = qgridder_utils.get_param(gridLayer, output_type = 'array', fieldName = 'IBOUND') hstart = qgridder_utils.get_param(gridLayer, output_type = 'array', fieldName = 'H0') # Get coordinates of centroids cx = qgridder_utils.get_param(gridLayer, output_type = 'array', fieldName = 'CX') cy = qgridder_utils.get_param(gridLayer, output_type = 'array', fieldName = 'CY') # ------------ Well data ------------ # select vector layer, add wells wellsLayer = ftools_utils.getVectorLayerByName('wells') wellsRowCol = qgridder_utils.get_ptset_centroids(wellsLayer, gridLayer, idFieldName = 'ID') wellsQvalues = qgridder_utils.get_feat_param(wellsLayer, idFieldName = 'ID', FieldName = 'Q') # init wells lay = 1 wellsData = [] for wellID in wellsRowCol.keys(): # fill wellData for wellID wellData = [lay] # lay wellData.append(wellsRowCol[wellID][0][0]) # row wellData.append(wellsRowCol[wellID][0][1]) # col wellData.append(wellsQvalues[wellID]) # Q # fill wellsData wellsData.append(wellData) # ------------ Drain data ------------ # select drain layer, add drains drainsLayer = ftools_utils.getVectorLayerByName('drains') drainsRowCol = qgridder_utils.get_pline_centroids(drainsLayer, gridLayer, idFieldName = 'ID') drainsEvalues = qgridder_utils.get_feat_param(drainsLayer, idFieldName = 'ID', FieldName = 'E') drainsCvalues = qgridder_utils.get_feat_param(drainsLayer, idFieldName = 'ID', FieldName = 'C') # init drains lay = 1 # drainsData [ [lay,row,col,e,c], drainsData = [] for drainID in drainsRowCol.keys(): for cellRowCol in drainsRowCol[drainID]: cellDrainData = [lay] # lay cellDrainData.append( cellRowCol[0] ) # row cellDrainData.append( cellRowCol[1] ) # col cellDrainData.append( drainsEvalues[drainID] ) # Elevation cellDrainData.append( drainsCvalues[drainID] ) # Conductance # fill drainData drainsData.append(cellDrainData) # ------------ Ghb data ------------ # select General-head boundary condition ghbLayer = ftools_utils.getVectorLayerByName('ghb') ghbRowCol = qgridder_utils.get_pline_centroids(ghbLayer, gridLayer, idFieldName = 'ID') ghbHvalues = qgridder_utils.get_feat_param(ghbLayer, idFieldName = 'ID', FieldName = 'H') ghbCvalues = qgridder_utils.get_feat_param(ghbLayer, idFieldName = 'ID', FieldName = 'C') # init ghb lay = 1 # ghbData [ [lay,row,col,e,c], ghbData = [] for drainID in ghbRowCol.keys(): for cellRowCol in ghbRowCol[drainID]: ghbDrainData = [lay] # lay ghbDrainData.append( cellRowCol[0] ) # row ghbDrainData.append( cellRowCol[1] ) # col ghbDrainData.append( ghbHvalues[drainID] ) # Elevation ghbDrainData.append( ghbCvalues[drainID] ) # Conductance # fill drainData ghbData.append(ghbDrainData) # ------------ Model settings ------------ # -- Space discretization # Model geometry # height (vertically) H = 100.0 # elevation of top top= 0 # elevation of bottom botm = -100 # Number of layers nlay = 1 # -- Aquifer properties, options for BCF package laycon = 0 # 0:confined 1:confined # Hydraulic conductivity k = 0.005 # Transmissivity T = k*H / nlay # if T field is not uniform: # Storage coefficient S = 0.01 # if S field is not uniform: # -- Time discretization steady = True # -- Boundary condition # For constant fixed-head (Dirichlet type boundary condition) # Set ibound to -1, initial values of h0 will be maintained. # fill layer_row_column_shead_ehead # ------------------------------------------------------------ # Modflow simulation with Flopy # ------------------------------------------------------------ name = 'demo1' # New model ml = Modflow(modelname=name, exe_name= home + '/Programmes/modflow/mf2005/src/mf2005', version='mf2005',silent=0) # DIS (Discretization) package discret = ModflowDis(ml,nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, top=top, botm=botm, steady=steady) # BAS package bas = ModflowBas(ml,ibound=ibound,strt=hstart) # BCF (Block-centered Flow) package. ibcf = -1 to save budget. bcf = ModflowBcf(ml, ibcfcb = 53, laycon=laycon, tran=T, sf1=S) # Well package. Note that wells data must be a list of list of layer_row_column_Q mfwel = ModflowWel(ml, layer_row_column_Q=[wellsData]) # Drain package mfdrn = ModflowDrn(ml, layer_row_column_elevation_cond=[drainsData]) # General head boundary condition package mfghb = ModflowGhb(ml, layer_row_column_head_cond=[ghbData]) # Output control package oc = ModflowOc(ml, words = ['head','budget'], save_head_every=1) # PCG package pcg = ModflowPcg(ml) # Write input files ml.write_input() # Run model ml.run_model3() # -- Read results mread_heads = ModflowHdsRead(ml,compiler='l') t, h, m = mread_heads.read_all(name+'.hds') h = np.array(h) h = np.swapaxes(h,2,3) h = np.swapaxes(h,1,2) mread_budget = ModflowCbcRead(ml,compiler='l') t, b, m = mread_budget.read_all(name+'.cbc') b = np.array(b) b = np.swapaxes(b,1,4) # ------------------------------------------------------------ # Post-processing # ------------------------------------------------------------ # Contour plot of heads with matplotlib plt.contour(cx[0,:],cy[:,0],h[0,0,:,:]) plt.show() # plot flow through right face plt.imshow(b[0,0,:,:,1]) # Compute budget for one cell cellWell = [1, 36, 28 ] # [lay,row,col] budgetCellWell = qgridder_utils.cells_budget(ml, b, [ cellWell ] ) # send simulation results back to Qgis qgridder_utils.data_to_grid(h[0,0,:,:], gridLayer, fieldName = 'H')