What is the cause for the degradation of environment? Capitalism, corruption, consuming society? - OVERPOPULATION!
Please, save the Planet - kill yourself...
Sometimes in Dango you need to compose query with several `Q(query) & Q(other_query) &...` and the number of such queries changes dynamically. Here is a quick example how to solve it:
import operator
ids = [1, 2, 3]
queries = [Q(some_m2m_relation__pk=pk) for pk in id]
result = SomeModel.objects.filter(
reduce(operator.and_, queries)
)
I wanted to have UUID fields in my Peewee-bsead models for the SQLite database. I quickly found ready-to-use code, but it lacked one important thing - automatic uuid generation. Here is the solution:
import uuid
from peewee import Field
class UIDField(Field):
db_field = 'uid'def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.default = uuid.uuid4
def db_value(self, value):
return str(value) # convert UUID to strdef python_value(self, value):
return uuid.UUID(value) # convert str to UUID
Recently I did a QGIS scripting job and here is the feedback from an extremly satisfied customer:
It was fantastic to work with Yury. He provided an excellent product, that went far beyond what I was expecting. I will certainly be contacting Yury in the future for any jobs that relate to GIS and python scripting. Communication was excellent and his ability to understand the job requirement was very impressive. A+++
Guys, if you are in need of geoprocessing tool for your project - don't hesitate to contact me ;-)
This is strange, but I was unable to find instruction about importing QGIS layers into PostGIS database with PyQGIS API. The PyQGIS cookbook has the example of exporting layers as .shp-files via QgsVectorFileWriter.writeAsVectorFormat() function and says that the other other OGR-supported formats are available to use in this function as well. PostGIS is supported by OGR so people got confused and try to use this function to import data to the PostGIS with no success and have to write generic import functions.
After of couple of hours of searching the internet for the solution I gave up and decided to find the answer the hard way: started to search through the source code for the DB Manager plugin that has this nice "import layer" feature. It took about 20 minuets to trace down the function that was in charge of the import and it was QgsVectorLayerImport.importLayer(). But the quest wasn't over yet! The documentation says nothing about provider names that are accepted by this function. "PostgreSQL" would be the obvious name for the provider as it is the name for the PostgeSQL provider in OGR, but it is not the case. I had to go through the source code of DB Manager again and luckily in comments (and there are quite a few of them in DB Manager: I didn't find a single doc-string there) the author wrote that it is "postgres" for the PostgreSQL. UPD: the function to list available providers is QgsProviderRegistry.instance().providerList().
Now here is the very basic example code of importing QGIS layer into PostGIS:
For my own project I needed to create a graph based on a Delauney triangulation using NetworkX python library. And a special condition was that all the edges must be unique. Points for triangulation stored in a .shp-file.
I was lucky enough to find this tread on Delauney triangulation using NetworkX graphs. I made a nice function out of it for point NetworkX graphs processing. This function preserves nodes attributes (which are lost after triangulation) and calculates lengths of the edges. It can be further improved (and most likely will be) but even at the current state it is very handy.
Example of use:
import networkx as nx
import scipy.spatial
import matplotlib.pyplot as plt
path = '/directory/'
f_path = path + 'filename.shp'
G = nx.read_shp(f_path)
GD = createTINgraph(G, show = True)
Code for the function:
import networkx as nx
import scipy.spatial
import matplotlib.pyplot as plt
def createTINgraph(point_graph, show = False, calculate_distance = True):
'''
Creates a graph based on Delaney triangulation
@param point_graph: either a graph made by read_shp() from another NetworkX's point graph
@param show: whether or not resulting graph should be shown, boolean
@param calculate_distance: whether length of edges should be calculated
@return - a graph made from a Delauney triangulation
@Copyright notice: this code is an improved (by Yury V. Ryabov, 2014, riabovvv@gmail.com) version of
Tom's code taken from this discussion
https://groups.google.com/forum/#!topic/networkx-discuss/D7fMmuzVBAw
'''
TIN = scipy.spatial.Delaunay(point_graph)
edges = set()
# for each Delaunay triangle
for n in xrange(TIN.nsimplex):
# for each edge of the triangle
# sort the vertices
# (sorting avoids duplicated edges being added to the set)
# and add to the edges set
edge = sorted([TIN.vertices[n,0], TIN.vertices[n,1]])
edges.add((edge[0], edge[1]))
edge = sorted([TIN.vertices[n,0], TIN.vertices[n,2]])
edges.add((edge[0], edge[1]))
edge = sorted([TIN.vertices[n,1], TIN.vertices[n,2]])
edges.add((edge[0], edge[1]))
# make a graph based on the Delaunay triangulation edges
graph = nx.Graph(list(edges))
#add nodes attributes to the TIN graph from the original points
original_nodes = point_graph.nodes(data = True)
for n in xrange(len(original_nodes)):
XY = original_nodes[n][0] # X and Y tuple - coordinates of the original points
graph.node[n]['XY'] = XY
# add other attributes
original_attributes = original_nodes[n][1]
for i in original_attributes.iteritems(): # for tuple i = (key, value)
graph.node[n][i[0]] = i[1]
# calculate Euclidian length of edges and write it as edges attribute
if calculate_distance:
edges = graph.edges()
for i in xrange(len(edges)):
edge = edges[i]
node_1 = edge[0]
node_2 = edge[1]
x1, y1 = graph.node[node_1]['XY']
x2, y2 = graph.node[node_2]['XY']
dist = sqrt( pow( (x2 - x1), 2 ) + pow( (y2 - y1), 2 ) )
dist = round(dist, 2)
graph.edge[node_1][node_2]['distance'] = dist
# plot graph
if show:
pointIDXY = dict(zip(range(len(point_graph)), point_graph))
nx.draw(graph, pointIDXY)
plt.show()
return graph
Here is an interesting post about 6.81 patch influence on heroes win rate in Dota 2. One point caught my attention: a drastic change in Axe win rate: from 52% to 48%; see image below.
Axe's win rate before and after 6.81 patch (29-th of April)
The only change Axe get in this patch was implementation of the preudo random approach for the determination of Counter Helix skill triggering instead of random. As was guessed in the article:
"...where you manage to Berserker's Call several enemies, you only have 3.2 seconds of being unconditionally hit by the enemy. In this time frame, the amount of successful Counter Helixes that use PRD is probably lower than with a true random. Hence the decrease in damage output and Win Rate."
Ok, nice guess! Lets test it!
Create a simulation to test the hypothesis
Import needed modules
import random
import numpy
from tabulate import tabulate
import collections
Lets create functions to calculate number of Counter Helix PROCs for given number of hits using random and pseudo random approaches. First one will be using random approach (notice that random.randrange() actually generates pseudo random set of values)
def randHelixProc(hits):
'''
checks if Counter Helix PROCs (chance = 17%) using random approach
@param hits: number of times Axe was hit
@return procs: number of times Counter Helix PROCs
'''
procs = 0
for i in xrange(hits):
chance = random.randrange(1, 101)
if chance < 18:
procs += 1
else:
continue
return procs
I don't know how exactly pseudo random PROCs are designed for Axe or for Dota 2 in general, but they say that the probability of the pseudo random event in Dota 2 in it initial state has lower chance to proc than the stated probability for the event and each time the event is not triggered when it could the chance for this event to occur is increased. Lets say some event triggers with the 50% chance. Suggested pseudo random approach can be realised with 3 consecutive checks that have different chances: 0%, 50% and 100% (average is 50%). When event finally triggers it reset chance for the event to initial value (0% in this case) and it all starts over again.
def pseudoRandHelixProc(hits):
'''
checks if Counter Helix PROCs (chance = 17%) using pseudo random approach
@param hits: number of times Axe was hit
@return procs: number of times Counter Helix PROCs
'''
treshold = 0
prob_list = [2, 5, 9, 14, 23, 47] # ~17% on average
procs = 0
for i in xrange(hits):
chance = random.randrange(1, 101)
try:
success = prob_list[treshold]
except:
success = prob_list[5]
if chance >= success:
treshold += 1
else:
procs += 1
treshold = 0
return procs
Check if the chances for PROCs are the same for the both functions. We should have 1700 procs of Counter Helix for 10000 attacs on Axe. Launch 100 simulations of 10000 hits for each function and compute average procs for random and pseudo random approaches.
rhelix_list = []
for i in xrange(100):
value = randHelixProc(10000)
rhelix_list.append(value)
numpy.mean(rhelix_list)
>>>1702.79
p_rhelix_list = []
for i in xrange(100):
value = pseudoRandHelixProc(10000)
p_rhelix_list.append(value)
numpy.mean(p_rhelix_list)
>>>1702.3
Output difference is negligible for random and pseudo random implementations.
Now its time to create a simulation function.
def simulation(times, hits):
'''
Computes average number of PROCs for given hits for random and pseudo random
approaches and difference between them.
@param times: number of times simulation will run
@param hits: number of hits to simulate
@return: average difference in pocs between random and pseudo random approaches;
table of simul results
'''
# create lists of results
diff_list = []
rand_mean_list = []
p_rand_mean_list = []
# run simulation
for hit in xrange(1, hits + 1):
rand_list = []
pseudo_rand_list = []
for t in xrange(times):
rand = randHelixProc(hit)
p_rand = pseudoRandHelixProc(hit)
rand_list.append(rand)
pseudo_rand_list.append(p_rand)
# compute statistics and populate lists of results
rand_mean = numpy.mean(rand_list)
rand_mean_list.append(rand_mean)
p_rand_mean = numpy.mean(pseudo_rand_list)
p_rand_mean_list.append(p_rand_mean)
diff = rand_mean - p_rand_mean
diff = round(diff, 2)
diff_list.append(diff)
# print average difference in PROCs
total_diff = sum(diff_list)
l = len(diff_list)
print 'average difference:', total_diff/l
print '#######################################################################################################'
# create table output for simulation results
out_dict = {}
out_dict['(1) cumulative hits'] = range(1, l + 1)
out_dict['(2) random mean PROCs'] = rand_mean_list
out_dict['(3) pseudo random mean PROCs'] = p_rand_mean_list
out_dict['(4) difference'] = diff_list
out_dict = collections.OrderedDict(sorted(out_dict.items()))
print tabulate(out_dict, headers = 'keys', tablefmt="orgtbl")
Lets run 100 simulations of 100 consecutive hits. Of course it is possible to run 10000 hits but it is unnecessary since every time the Counter Helix is triggered we jump back to the first row from the table below (result of the simulation). As was already mentioned:
"[if] ...you manage to Berserker's Call several enemies, you only have 3.2 seconds of being unconditionally hit by the enemy"
If Axe was able to catch 3 enemies, together they will hit him about 9-12 times. This means that with the pseudo random approach he will spin 0.3-0.4 times less than with random approach. It will cause him to deal ~(205*0.35)*3 = 215,25 (or 71.75 per hero) less damage in engagement.
So the hypothesis was true - pseudo random approach caused the decrease in damage output on the short time interval even if total number of PROCs during the game stayed the same.
Script for unifying extent and resolution is now able to work with multi-band rasters. Also I fixed an error that caused output rasters to have different resolutions in some cases.
The script was designed to be used within Processing module of QGIS. This script will make two or more rasters of your choice to have the same spatial extent and pixel resolution so you will be able to use them simultaneously in raster calculator. No interpolation will be made - new pixels will get predefined value. Here is a simple illustration of what it does:
Modifications to rasters A and B
To use my humble gift simply download this archive and unpack files from 'scripts' folder to your .../.qgis2/processing/scripts folder (or whatever folder you configured in Processing settings). At the next start of QGIS you will find a 'Unify extent and resolution' script in 'Processing Toolbox' in 'Scripts' under 'Raster processing' category:
If you launch it you will see this dialogue:
Main window
Note that 'Help' is available:
Help tab
Lets describe parameters. 'rasters' are rasters that you want to unify. They must have the same CRS. Note that both output rasters will have the same pixel resolution as the first raster in the input.
Raster selection window
'replace No Data values with' will provide value to pixels that will be added to rasters and replace original No Data values with the value provided here. Note that in output rasters this value will not be saved as No Data value, but as a regular one. This is done to ease feature calculations that would include both of this rasters, but I'm open to suggestions and if you think that No Data attribute should be assigned I can add corresponding functionality.
Finally you mast provide a path to the folder where the output rasters will be stored in 'output directory' field. A '_unified' postfix will be added to the derived rasters file names: 'raster_1.tif' -> 'raster_1_unified.tif'
If CRSs of rasters won't match each other (and you will bypass built-in check) or an invalid pass will be provided a warning message will pop up and the execution will be cancelled:
Example of the warning message
When the execution is finished you will be notified and asked if rasters should be added to TOC:
Happy New Year!
P.S. I would like to thank Australian government for making the code they create an open source. Their kindness saved me a couple of hours.
I decided to make my script for counting unique values in raster more usable. Now you can use it via SEXTANTE in QGIS. Download script for SEXTANTE and extract archive to the folder that is is intended for your Python SEXTANTE scripts (for examlple ~./qgis2/processing/scripts). If you don't know where this folder is located go Processing -> Options and configuration -> Scripts -> Scripts folder, see the screenshot:
Now restart QGIS and in SEXTNTE Toolbox go to Scripts. You will notice new group named Raster processing. There the script named Unique values count will be located:
Launch it and you will see its dialogue window:
Main window
Note that Help is available:
Either single- or multi-band rasters are accepted for processing. After the raster is chosen (input field) one need to decide whether to round values for counting or not. If no rounding is needed - leave 'round values to ndigits' field blank. Otherwise enter needed value there. Note that negative values are accepted - this will round values to ndigits before decimal point.
When the set up is finished hit the Run button. You will get something like this:
Recently I demonstrated how to get histograms for the rasters in QGIS. But what if one need to count exact number of the given value in a raster (for example for assessment of classification results)? In this case the scripting is required.
UPD: the script below was transformed in more usable SEXTANTE script, see this post.
We'll use Python-GDAL for this task (thanks to such tutorials as this one) it is extremely easy to learn how to use it). The script I wrote seems to be not optimal (hence its working just fine and quickly so you may use it freely) due to it is not clear to me whether someone would like to use it for the floating data types (which seems to be not that feasible due to great number of unique values in this case) or such task is performed only for integers (and in this case code might be optimised). It works with single and multi band rasters.
How to use
The code provided below should be copied into a text file and saved with '.py' extension. Enable Python console in QGIS, hit "Show editor" button in console and open the script file. Then replace 'raster_path' with the actual path to the raster and hit 'Run script button'. In the console output you will see sorted list of unique values of raster per band:
Script body (to the right) and its output (to the left)
Script itself
#!/usr/bin/python -tt
#-*- coding: UTF-8 -*-
'''
#####################################################################
This script is designed to compute unique values of the given raster
#####################################################################
'''
from osgeo import gdal
import sys
import math
path = "raster_path"
gdalData = gdal.Open(path)
if gdalData is None:
sys.exit( "ERROR: can't open raster" )
# get width and heights of the raster
xsize = gdalData.RasterXSize
ysize = gdalData.RasterYSize
# get number of bands
bands = gdalData.RasterCount
# process the raster
for i in xrange(1, bands + 1):
band_i = gdalData.GetRasterBand(i)
raster = band_i.ReadAsArray()
# create dictionary for unique values count
count = {}
# count unique values for the given band
for col in range( xsize ):
for row in range( ysize ):
cell_value = raster[row, col]
# check if cell_value is NaN
if math.isnan(cell_value):
cell_value = 'Null'
# add cell_value to dictionary
try:
count[cell_value] += 1
except:
count[cell_value] = 1
# print results sorted by cell_value
for key in sorted(count.iterkeys()):
print "band #%s - %s: %s" %(i, key, count[key])
I fixed an issue of my utility for polygon width calculation in a given direction caused by redundant geometry created in certain circumstances. 'Fixed' - somewhat a too strong claim for a nasty about 100 lines long workaround-function, where I had to comment almost every line to trace what is going on (optimisation is left for the feature). Nevertheless full functionality of the utility is now available! Updated version is here.
Another note: check validity of your polygons geometries beforehand to obtain correct results.
I have created an "Azimuth-Width" utility recently. It is designed to compute (for now) maximum or minimum width of the given polygon(s) in the given direction. Corresponding QGIS plugin coming soon (or not too soon), so this utility may help these who are in need right now.
If you are non-Windows user ensure that command "QgsApplication.setPrefixPath( qgis_prefix, True)" in the utility file have a valid path to QGIS installation. If not - modify "qgis_prefix" to set correct value. It is "/usr" now - must work in most of cases.
General usage:
Copy width.py to the directory with a shh-file containing polygons.
Using console navigate to the directory. In console type:
:~> python ./width.py [file to analyse] [field to store values] [azimuth (decimal degrees)] [mode ('min' or 'max')] [mode-2 ('abs', or 'rel')] [algorithm ('byStep' or 'byVertex' or 'Mix')] [step (real numver, CRS units; for 'byStep' and 'Mix' only)]
Example:~> python ./width.py poly.shp width 285.9 max abs byStep 1.3
Where:
./width.py - name of this utility file.
[file to analyse] - a shp-file with polygons to analyse.
[field to store values] - if it does not exist it will be created.
[azimuth] - a direction for width calculation. Accepts decimal degrees from 0 to 360.
[mode] - type of width to calculate. Currently 'min' (returns minimum width value in given direction) and 'max' (returns maximum value in given direction) modes are available. Mode 'min' used with 'byVertex' algorithm will return minimum polygon width different then 0.0. This will be different (greatly most time) to 'min' mode and 'byStep' or 'Mix' algorithm.
[mode-2] - if polygon is not convex it may have several segments in given direction. If you want to take sum of the segments of the result - use 'abs', if you want only the longest (shortest) segment - use 'rel'. Fixed:Note that currently 'rel' sometimes will provide incorrect results for non-convex polygons and for convex polygons with redundant vertices on its edges. This issue will be solved when the behaviour of the intersection() command will be changed or when I will implement a workaround for it.
[algorithm] - algorithm that will be used: 'byStep', 'byVertex' or 'Mix'. "byStep" algorithm (improved version of the general idea that described here) will take provided value (shp-file CRS units) and swipe the polygons line with the given step along the bounding box. Precision depends on the step. Speed and precision depends on step - lower step mean more precise result but computation will take more time. "byVertex" (general idea described here) will intersect polygon with the line only in vertexes of the given polygon. For convex polygons this algorithm will be faster and more precise using 'max' mode then "byStep" algorithm. But this algorithm is not suitable for 'min' mode. "Mix" algorithm will use both "byVertex" and "byStep" algorithm so in some cases it will be most precise but will consume even more time than 'byStep'.
[step] - step for "byStep" and "Mix" algorithm. Takes decimal values in shp-file CRS units. Lower step means more precision but more computation time.
Some advises:
Use Equal Area projections.
Make sure your polygons have correct geometry. Otherwise the calculations will be incorrect.
'byStep' algorithm should be faster when polygon have enormous number of vertices.
If you have a large variety in polygons area e.g. a continent and an island to save computation time you may define big step (like 1000 m or so) to swipe through continent faster. It will save computation time and a width for the small island will be calculated even if step is grater then any side of island's' Bounding Box: step for the island in this case will be 1/100 of the shortest side of it's Bounding Box.
