# Change the following if necessary: Path = "C:/GISIranData/" ; Number_per_dot = 75 # Initializing variables Points = 0 ; Crops = 0 ; Outlets = 0 ; Adjus = 0; Toler = 0 TotArea = 0 ; TotRev = 0 ; IterNo = 0 TotDemand = [] # Dimensions - Crops Id = [] ; Area = [] ; Frac = [] # Dimensions - Points Trans = [] # Dimensions - Outlets Suit = [] ; Capac = [] ; Prod = [] ; ZoneRev = [] ; Select = [] # - Points X Crops Price = [] ; NewPrice = [] # - Outlets X Crops Dist = [] # - Points X Outlets Rev = [] # - Crops X Outlets # ---------------------- # 'Read' values from a .csv table def read(Name): global Crops,Points,Outlets,Adjus,Toler,TotDemand,Trans global Price,NewPrice, Path pTemp = [] ; import csv i = -1; stg = Path + Name + ".csv" with open(stg) as csvDataFile: csvReader = csv.reader(csvDataFile) for row in csvReader: i = i+1 if i == 1: if Name == "Parameters": Crops = int(row[0]) ; Points = int(row[1]) Outlets = int(row[2]) ; Toler = int(row[3]) Adjus = float(row[4]) if Name == "TotalDemand": TotDemand.append(int(row[0])) TotDemand.append(int(row[1])) if Name == "Outlets": for row in csvReader: Trans.append(int(row[3])) ; jj = 3 ; pTemp = [] for j in range(Crops): jj = jj + 1 ; pTemp.append(int(row[jj])) Price.append(pTemp) ; jj = 3 + Crops ; pTemp = [] for j in range(Crops): jj = jj + 1 ; pTemp.append(int(row[jj])) NewPrice.append(pTemp) if Name == "Suitabilities": for row in csvReader: stemp = [] ; l = 2 for j in range(0,Crops): l = l + 1 ; stemp.append(int(row[l])) Suit.append(stemp) #-------------- # 'Observe' some variables in a layer def observe(Name): global Area,TotArea,Frac,Dist,Capac,Id if Name == "Pattern": layer = QgsProject.instance().mapLayersByName(Name)[0] TotArea = 0 for i in range(Points): feat = layer.getFeature(i) TotArea = TotArea + int(feat['Area']) for i in range(Points): feat = layer.getFeature(i) idtemp = [feat['Id']]; Id.append(idtemp) Area = int(feat['Area']) tempno = float((Area/TotArea) * Crops) Frac.append(tempno) ; DistTemp = [] ; tempstg = "Hello" for j in range(Outlets): tempstg = 'Dist_' + str(j+1) ; dtemp = int(feat[tempstg]) DistTemp.append(dtemp) Dist.append(DistTemp) ; CapacTemp = [] for j in range(Crops): ctemp = int(Frac[i] * TotDemand[j]) ctemp = int((ctemp * (1 + ((Suit[i][j]-5)* 0.2)))) CapacTemp.append(ctemp) Capac.append(CapacTemp) #------------ # Calculate optimal crop's production in each zone def Find_soln(): global Rev,TotRev,Prod,TotDemand,Adjus,Price,NewPrice,Toler global Crops,Outlets,Capac,Trans,Select Prod =[]; Select = [] ; TotRev = 0 ; Return = [] for i in range(0,Points): Stemp = [] for j in range(0,Crops): Stemp.append(0) Select.append(Stemp) for i in range(0,Points): Rev = [] for j in range(0,Crops): OTemp = [] for k in range(0,Outlets): rtemp = int(Capac[i][j] * NewPrice[k][j]) rtemp = rtemp - int(Capac[i][j] * Trans[k]) OTemp.append(rtemp) Rev.append(OTemp) M = - 100000 ; jj = 0 for j in range(0,Crops): for k in range(Outlets): if Rev[j][k] > M: M = Rev[j][k] ties = 0 ; Return = [] ; MM = int(M - ((Toler/100) * M)) for j in range(0,Crops): tempReturn = 0 for k in range(Outlets): if Rev[j][k] > MM: ties = ties + 1 ; Select[i][j] = 1 tempReturn = Rev[j][k] ; break Return.append(tempReturn) ProdTemp = [] for j in range(0,Crops): ptemp = 0 if Select[i][j] == 1: TotRev = TotRev + int((Return[j]/ties)) ptemp = int(Capac[i][j]/ties) ProdTemp.append(ptemp) Prod.append(ProdTemp) # Find consequent price changes for j in range(0,Crops): TotProd = 0 for i in range(0,Points): TotProd = TotProd + Prod[i][j] print("Commodity",j+1," Demand =",TotDemand[j]," Prodn. =",TotProd) percentdiff = int(((TotProd - TotDemand[j]) * 100) / TotDemand[j]) Miss = percentdiff * Adjus for i in range(0,Points): change = int((Miss / 100) * Capac[i][j]) Capac[i][j] = int(Capac[i][j] - change) for k in range(0,Outlets): change = int((Miss / 100) * Price[k][j]) NewPrice[k][j] = int(NewPrice[k][j] - change) #------------ # Insert a new field (Name) def New_fields(Layer_name,Name,Iter): global Path #tempstr = "C:\GISData\\" tempstr = Path + Layer_name + ".shp" layer = QgsVectorLayer(tempstr,"park","ogr") features=layer.getFeatures() from PyQt5.QtCore import QVariant layer_provider=layer.dataProvider() for k in range(Crops): tempstg = Name + str(k+1) + Iter layer_provider.addAttributes([QgsField(tempstg,QVariant.Int)]) layer.updateFields() with edit(layer): for feature in layer.getFeatures(): layer.updateFeature(feature) layer.commitChanges() #------------ # Change values in a field (Name) to those in the (Variable) array def Change_fields_vals(Layer_name,Name,Variable,Iter): global Prod, Path ; tempstr = Path + Layer_name + ".shp" layer = QgsVectorLayer(tempstr,"park","ogr") features=layer.getFeatures() from PyQt5.QtCore import QVariant layer_provider=layer.dataProvider() with edit(layer): i = -1 for feature in layer.getFeatures(): i = i + 1 for k in range(Crops): ii = k + 1 ; b = Name + str(ii) + Iter tnber = int(Variable[i][k]) ; feature[b] = tnber layer.updateFeature(feature) # ---------- # Generates a new layer in which the optimal production of each # crop in each zone is shown def Record_optimal(Layer_name,Iter): global Prod, Path New_fields(Layer_name,"Prodn_",Iter) Change_fields_vals(Layer_name,"Prodn_",Prod,Iter) IterF = Iterations + 1 if Iter == "_"+ str(IterF): tempstr = Path + Layer_name + ".shp" iface.addVectorLayer(tempstr,"Optimal", "ogr") # ---------- # Make dot maps def Make_dot_maps(): import random ; global Prod,Crops,Path,Number_per_dot src = Path + "Pattern.shp" tractLyr = QgsVectorLayer(src, 'Pattern', "ogr") for i in range(0,Crops): xx = i + 1 ; tempstr = "Optimal_" + str(xx) popLyr = QgsVectorLayer('Point?crs=epsg:4326', tempstr , "memory") features = tractLyr.getFeatures() vpr = popLyr.dataProvider() j = -1 for feature in features: j = j + 1 ; dotFeatures = [] density = Prod[j][i] / Number_per_dot found = 0 ; dots = [] ; g = feature.geometry() minx = g.boundingBox().xMinimum() miny = g.boundingBox().yMinimum() maxx = g.boundingBox().xMaximum() maxy = g.boundingBox().yMaximum() while found < density: x = random.uniform(minx,maxx) y = random.uniform(miny,maxy) pnt = QgsPointXY(x,y) if g.contains(pnt): dots.append(pnt) ; found += 1 geom = QgsGeometry.fromMultiPointXY(dots) f = QgsFeature() ; f.setGeometry(geom) dotFeatures.append(f) vpr.addFeatures(dotFeatures) ; popLyr.updateExtents() QgsProject.instance().addMapLayers([popLyr]) # -------------------------------------------------------------- # Run script Iterations = 25 ; read("Parameters") ; read("TotalDemand") read("Outlets") ; read("Suitabilities") ; observe("Pattern") Find_soln() ; print("Total Revenue =",TotRev) ; ii = 1 for i in range(Iterations): IterNo = ii ; ii = ii + 1 ; print("Iteration",IterNo,":") Find_soln() ; print("Total Revenue =",TotRev) It = "_" + str(ii) if i == Iterations - 1: Record_optimal("Pattern",It) Make_dot_maps() # --------------------------------------------------------------