Misanthrope's Thoughts

Chaining AND operators in Django

noreply@blogger.com (Yury Ryabov) — Tue, 16 Oct 2018 15:10:00 +0000

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)
)

How to Configure Pulseaudio to Work With Subwoofer

noreply@blogger.com (Yury Ryabov) — Tue, 07 Feb 2017 18:41:00 +0000

What I hate about Pulseaudio is its crappy default configuration. After fresh installation of SUSE Leap 42.2 I encountered an issue that subwoofer didn't work with Amarok... again((( Lucky me, the first link that google provided me with was link to my own(!!!) thread that I started about 6 years ago (sick!!!). In 2011 I encountered the same issue, asked for help, but solved on my own and answered for my own request. And now future me was able to find that long forgotten thread of mine. That's so sweet to recieve a messege to miself )))

This time I will put the answer here to remember it better ;-)

The problem solved, tnx to this post.
In /etс/pulse/daemon.conf I've made following changes:
Code:
 
enable-lfe-remixing = yes
default-channel-map = front-left,front-right,rear-left,rear-right,front-center,lfe

Environmental Monitoring App

noreply@blogger.com (Yury Ryabov) — Fri, 23 Sep 2016 19:16:00 +0000

Python, R, Qt, peewee, bokeh, pandas, SQLite plus couple of sleepless nights and here you are a сute app for the environmental monitoring needs )))

Main window of the application

UUID Field for the Peewee ORM

noreply@blogger.com (Yury Ryabov) — Tue, 23 Aug 2016 22:41:00 +0000

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 str
    def python_value(self, value):
        return uuid.UUID(value)  # convert str to UUID

A satisfied customer

noreply@blogger.com (Yury Ryabov) — Sun, 08 Nov 2015 23:33:00 +0000

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 ;-)

My QGIS Processing Scripts at GitHub

noreply@blogger.com (Yury Ryabov) — Thu, 22 Oct 2015 21:38:00 +0000

This is probably my shortest post ever.

All my QGIS processing scripts (R and Python) and models that I already blogged about, plus some extra are now available at GitHub.

A Quick Map With QGIS and OSM

noreply@blogger.com (Yury Ryabov) — Thu, 18 Jun 2015 16:27:00 +0000

What I love about QGIS is that one is able to create a nice map quickly. The other day I was asked to make a situation map for the project we are working on to include it into presentation. Аll I had was a laptop with no relevant spatial data at all, but with QGIS installed (I even had no mouse to draw something). Though it was more than enough: I loaded OSM as a base layer and used annotation tool to add more sense to it. Voilà:

How to Import Layer into PostGIS Database Using PyQGIS

noreply@blogger.com (Yury Ryabov) — Wed, 10 Jun 2015 21:03:00 +0000

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:

uri = "dbname='test' host=localhost port=5432 user='user' password='password' key=gid type=POINT table=\"public\".\"test\" (geom) sql="
crs = None
# layer - QGIS vector layer
error = QgsVectorLayerImport.importLayer(layer, uri, "postgres", crs, False, False)
if error[0] != 0:
    iface.messageBar().pushMessage(u'Error', error[1], QgsMessageBar.CRITICAL, 5)

A Note About a "#IND"-Error in OfreoToolbox Classification Algorithms

noreply@blogger.com (Yury Ryabov) — Wed, 01 Apr 2015 13:07:00 +0000

One of the optional input for the creation of the classification models in OrfeoToolbox is XML image statistic file that is produced by Compute Image Second Order Statistics tool. If you opt to calculate this statistics - necessarily check the created XML file. If you will see "#INF" record in band statistics (i.e. "1.#INF") - replace "#INF" with something like "0e+32" ("1.#INF" -> "1.0e+32"), or (a better solution) - calculate the statistics for this problematic band independently (probably with another tool) and completely replace the value.

If you leave "#INF" record in XML file, models for classification that will be created based on it will have "#IND" records and when you will try to run classification based upon these models - the process will be terminated complaining about unrecognised "#IND"-values, because values must be numeric and not characters.

I don't know what causes the "#INF" records to appear in the first place, maybe not properly defined no-data values? In some cased it is impossible to provide correct no-data value due to limitations of the input field (I use OTB from within Processing module for QGIS).

Why are natural environments often degraded in the name of economic progress even when this invites disaster, and how can this be overcome?

noreply@blogger.com (Yury Ryabov) — Thu, 26 Mar 2015 10:57:00 +0000

This was one of discussion questions of the Disasters and Ecosystems MOOC.

Actually the answer is simple. The formula for successful environmental degradation consists of 2 variables - overpopulation and capitalism.

When there are a lot of people - most of them a poor, uneducated and hungry. When you are hungry you will do everything to become less hungry today even if it can potentially lead to negative consequences tomorrow, which you may not even foresee if you are uneducated.

Humans are good in adaptation. When the adaptation is strong enough it leads to abuse (for example, if you are well adopted at the stock market you start abusing it to increase your profit even if it will cost dearly to the other stakeholders - people value their own well-being much more than the other's and of course much more than the well-being of environment especially when they know that their own impact seems negligible compared to impact of the entire population). When you live in condition of free market of capitalistic world - you are your only hope for not being hungry (or being more wealthy) now. And as you know from the economic theory - the capitalist economy needs a constant grows of consumption and production - so you need more and more resources to just sustain the economy. In conditions of capitalist market people value today's profit much more than losses of tomorrow.

You see - the capitalist economy needs people to consume more and more; more people - more consumption; more people - more poverty and lack of education; more hungry uneducated people people - more people willing to do anything to survive now and don't even bother themselves about the future.

Overpopulation and a consumption society (created by capitalist economy) inter-stimulate each other and destroy the environment for the today's profits or food and doesn't care much of the consequences of tomorrow because most are either uneducated or doesn't care at all plus you have to live through today to face consequences of your actions tomorrow (a day-by-day living).

Obviously there are 3 steps to improve the situation:

Decrease the population.
Educate people.
Create new sustainable economy model that would equally value tomorrow's losses and today's profits, and would not rely on constantly increasing consumption.

Pan-sharpening Using R

noreply@blogger.com (Yury Ryabov) — Sun, 08 Feb 2015 14:39:00 +0000

In my previous post I described how to perform pan-sharpening using OrfeoToolbox and QGIS. This time I will show you how to do this in R. At the bottom you will find several functions I wrote on top of the 'raster' package that allow a convenient pan-sharpening in R.

Motivation

You may wonder why I even bothered myself with pan-sharpening in R when I already have a nice model for pan-sharpening in QGIS. See, one can't control the data-type of the imagery returned by pan-sharpening that involves OTB. This leads to some unpleasant consequences: during pan-sharpening one will get floating point pixel values even if in initial values were integers. So for example a 600 MiB multi-spectral imagery (with integer pixel values) after pan-sharpening will grow to 5.2 GB. But if we will change datatype of the resulting imagery to force it store only integers it size will be reduced from 5.2 to 2.8 GB which is a huge difference. 'raster' package in R allows to control output datatype. Plus using R you can play with different filtering options to play with.

The Theory

In OTB pan-sharpening is performed using following well-known formula:

Where i and j are pixels indices, PAN is the panchromatic image, XS is the multi-spectral image and PAN_smooth is the panchromatic image smoothed with a kernel to fit the multi-spectral image scale.

We will implement exact the same approach using 'raster' package for R.

Code Usage and Result

As pan-sharpening is the type of procedure that will reoccur over some time I decided to write generic functions for pan-sharpening itself and for saving the results to have easier time in future.

The usage is as simple as:

pan <- raster('pan.tif')
multi <- brick('multi.tif')
pansharp <- processingPansharp(pan, multi)
output_path <- 'path_and_filename_without_extention'
saveResult(pansharp, output_path)

Here you are the example results from the script and from the OTB model for one of the illegal landfills in Russia:

Initial multi-band raster

Initial panchromatic raster

Result of pan-sharpening using R script

Result of pan-sharpening using OTB

Which output do you like better: from OTB or R? Comparing both output results you can notice that the output from R bears heavier filtering markings than the one from OTB. On the other hand R output has more hues of the green colours which actually helps in distinguishing different types of vegetation. As you will see in the code - one can easily adjust or modify procedure of filtering panchromatic raster (extractLPF() function) in order to get desired output.

The code

library(raster)


# Create needed functions -------------------------------------------------

pansharpFun <- function(raster){
    '
    This function pansharpens a raster
    ' 
    # @param raster - Raster object with 3 bands (to-be-pansharpened, high-res and low-frequency component of the high-res image)
    # @param band - band numver, integer
    # @return pansharpened_raster - pansharpened Raster object
    # pansharp = Lowres * Highres / LPF[Highres]
    
    pansharpened_raster <- (raster[,1] * raster[,2]) / raster[,3]
}

extractLPF <- function(pan, multi, filter = 'auto', fun = mean) {
    '
    Returns a low-frequency component of the high-resolution raster by the
        filter adjusted to the low-resolution raster
    '
    # @param pan - a high-resolution panchromatic raster - Raster object
    # @param multi - low-resolution raster to be pansharpened - Raster object
    # @param filter - a smoothing wondow - matrix
    # @param fun - a function to process filter (part of the focal() function)
    # @return LPF - a low-frequency component of the high-resolution raster - Raster object
    
    # Adjust filter size 
    if (filter == 'auto') {
        pan_res <- res(pan) # (x, y) resolution of the panchromatic raster in CRS units (?)
        multi_res <- res(multi) # (x, y) resolution of the lowres raster in CRS units (?)
        x_res_ratio <- round(multi_res[1]/pan_res[1])
        y_res_ratio <- round(multi_res[2]/pan_res[2])
        total <- x_res_ratio + y_res_ratio
        filter <- matrix(1, nc = x_res_ratio, nr = y_res_ratio)
        
        # Enshure that the matrix has an uneven number of colums and rows (needed by focal())
        if (nrow(filter)%%2 == 0) {
            filter <- rbind(filter, 0)
        } 
        if (ncol(filter)%%2 == 0) {
            filter <- cbind(filter, 0)
        }
        
    LPF <- focal(pan, w = filter, fun = fun) # low-frequency component
    
}  


processingPansharp <- function(pan, multi, filter = 'auto', fun = mean){
    ' 
    Pansharpening routine
    '
    # @param pan - a high-resolution panchromatic raster - Raster object
    # @param multi - low-resolution raster to be pansharpened - Raster object
    # @param filter - a smoothing wondow - matrix
    # @param fun - a function to process filter (part of the focal() function)
    # @return pansharp - pansharpened 'multi' raster - Raster object

    # Check if input parameters are valid - we can loose a lot of time if some of the inputs is wrong
      
    LPF <- extractLPF(pan, multi, filter, fun)
        
    multi <- resample(multi, pan) # resample low-resolution image to match high-res one
        
    all <- stack(multi, pan, LPF)
    
    bands <- nbands(multi)
    pan_band <- bands + 1
    lpf_band <- bands + 2
    
    # Pansharpen layers from low-resolution raster one by one
    pansharp_bands <- list()
    for (band in 1:bands) {
        subset <- all[[c(band, pan_band, lpf_band)]]
        raster <- calc(subset, pansharpFun)
        pansharp_bands[[band]] <- raster
    }
    
    pansharp <- stack(pansharp_bands)
}

saveResult <- function(raster, path, format = 'GTiff', datatype = 'INT2S'){
    '
    Saves Raster object to location
    '
    # @param raster - raster to be saved - Raser object
    # @param path - path including filename without extention - string
    # @param format - format of the output raster accordingly to writeRaster() function - string
    # @param datatype - datatype of the raster accordingly to writeRaster() - string
    
    writeRaster(raster, 
                path, 
                format = format, 
                datatype = datatype, 
                overwrite = T)
}


# Do pansharpening --------------------------------------------------------


pan <- raster('pan.tif')
multi <- brick('multi.tif')
pansharp <- processingPansharp(pan, multi)
output_path <- 'r_pansharp-new' # includes path and filename but not the extention
saveResult(pansharp, output_path)

Pansharpening in QGIS Using QTB

noreply@blogger.com (Yury Ryabov) — Sat, 31 Jan 2015 16:25:00 +0000

UPD: also you may want to check out my post on Pansharpening using R.

For a long time I wanted to play with the OrfeoToolBox instruments (and its GUI named Monteverdi) - a set of powerfull tools to process remote-sensing imagery. Finally I got the opportunity - I needed to perform pansharpening of the World-View-2 scene and I decided to make it using OTB modules available via QGIS Processing toolbox.

To make OTB modules available in QGIS you need to install it on your system. Official site provides good instructions on how to do it. When you have OTB installed you need to enable it in QGIS. Go Processing -> Options and configuration -> Providers, activate OTB and provide folder information if needed:

Pansharpening is made in two steps: 1) resampling of the low resolution raster via Superimpose sensor tool; 2) pansharpening using Pansharpening tool. Alternatively you can do it using one command line (no QGIS or other GUI needed) in your OS's console:

otbcli_BundleToPerfectSensor -inp pan_image -inxs xs_image -out output_image.

Performing two steps instead of just one is super-boring! Lets create a model that will allow us to perform pansharpening in QGIS in one step (virtually). Go to Processing Toolbox -> Models -> Tools and start creating new model:

In the model creation window fill the fields of the model name and group name. Then add two rasters as inputs and name them: High_resolution_panchrom_raster and Low_resolution_raster. Now add the module named Superimpose sensor (Orfeo Toolbox -> Geomentry -> Superimpose sensor) and configure it to resample the Low_resolution raster: use High_resolution_panchrom_raster as input in Reference input field and Low_resolution_raster as input in The image to reproject field.

Add Pansharpening (rsc) (Orfeo Toolbox -> Geomentry -> Superimpose sensor) module to model. Use High_resolution_panchrom_raster as input in the Input PAN Image and output of the Superimpose sensor as the input in the Input XS Image field. Give a name (pansharpened_OTB) for the Output image to let model know that this is the final stage of the processing.

Here how our model looks like in Model builder:

And this is how it looks like when you launch it:

Also you can skip model creation process and just download model that I created. You need to paste files from archive's models folder into your /.qgis2/processing/models directory.

If you will have issues running OTB Pansharpening it is likely that they will be covered in following gis.stackexchange topics: OTB Pansharpening Error: Adapter for adaptPansharpening-bayes not found and What causes OTB pansharpening ERROR: Inputs do not occupy the same physical space?

A Dialogue With Bot

noreply@blogger.com (Yury Ryabov) — Wed, 28 Jan 2015 22:13:00 +0000

Back in the days when ICQ was popular spam bots added you to their contact lists dozens times a day. But they were brute and ugly. At first they spew their spam links at your face when you added them, then they simply started to insert spam links in the body of the invitation. ICQ is long dead and spam-bots never bothered me in modern protocols... until some time ago.

A contact with the nick-name dazzling.chick9 (according to the profile - a 22 (or so) years old girl) added my on skype. Well, I felt it was bot, but added it anyway. For reasons unknown it was unable to chat with me several times. It sent me messages like 'hi' or 'hello', but didn't respond to me further. These greetings were sent the second I connected to the skype - it was obvious that these messages could be sent only by bot. But I waited patiently and today it finally started to talk to me:

[13:51:52] Amanda: hi

[13:52:26] Me: hi

[13:52:55] Amanda: how are you? i'm Amanda, I saw your Skype in the skype directory. wanna chat ?

[13:53:26] Me: chat about what?

[13:53:56] Amanda: I cannot wait until the summer is actually here, i can tone up my body more

[13:55:46] Me: wow, you really think you can pass turing test?

[13:55:53] Me: lol

[13:56:14] Amanda: i love too show what i have, i love the attention from guys

[13:56:15] Amanda: what should i wear? yoga pants or spandex?

[13:58:48] Me: a spanch bob outfit maybe? this bot is lame cause it doesn't attempt to understand what I'm talking about

[13:59:18] Amanda: well, i sort of have my just a top on atm so you wanna c how it looks?

[14:00:06] Me: yeah baby, send me that spam/scam/fishing link you prepared for me

[14:00:38] Amanda: my brothers computer does have a cam all you have 2 do is go two http://goo.gl/scam_or_virus_link

[14:00:55] Me: bingo!

[14:01:27] Amanda: I love this way because its always works + its free. this way no one else sees. If my mom saw what i want too show you i would kill myself.

[14:01:50] Me: I would love to see that

[14:02:03] Me: (how you kill yourself)

[14:02:22] Amanda: i sorta want you two chat to me kinda dirty. i love being talked down to.

[14:03:17] Me: Am I not dirty enough, stupid bot?

[14:03:49] Amanda: you make me really wet, listen i got it all ready, i will b there 4 you

[14:05:05] Me: great! You know, I knew you are bot before you asked me "how are you?"

[14:05:21] Amanda: k

[14:05:31] Me: do you want to know how I knew it?

[14:06:33] Me: if you want - just ask me ;-)

[14:07:26] Me: Any way it was fun, dude)))

I wonder how stupid people could be to sink that some random chick would approach them with a free striptease offer? Though I was amused with participation in the real-life Turing test. This bot imitated typing: it inserted message and waited about 30 second to actually send it so I would see 'typing...' indication in chat window, though this disguise was ruined in its 4-th and 5-th messages which had only a second between them. Also the story told in conversation does not correspond to age indicated in profile - only a teenagerwould be concerned about her mom and use a brother's computer, not a 22+ woman. Amanda stopped responding and I blocked her and reported abuse. Also I never dared to follow the link despite I'm using Linux, what if I really missed a free striptease? ;-)

How to Predict Where Will Next Disaster Strike?

noreply@blogger.com (Yury Ryabov) — Sun, 18 Jan 2015 20:23:00 +0000

It is amusing coincidence that another MOOC that I took this week (Geospatial Intelligence & the Geospatial revolution) mentioned [natural] disasters. About the other course see my recent Disasters: Myth or the Reality post.

In Geospatial Intelligence they gave a weird assignment: one need to mark the location on the world map where the next international natural disaster will occur O_o. This is not and easy task by any means and the lecturer suggested to use one's 'gut feeling' if one's knowledge is insufficient (I suppose it is close to impossible to find someone who can make such a prediction taking into account all the types of the disasters). Though the link to the International Disasters Database was given, so I accepted the challenge (to make a data-driven prediction). To predict the exact location of the next disaster one would need a lot of data - far more that you can get out of that database so my goal was to make prediction at the country level. (BTW the graphs from my post about disasters seems to be based on the data from this database - I saw one of them at that site)

I passed a query to the database and saved the output to process it with R. The dataframe looks like this:

Example of disasters dataset

So how to predict the country with the next disaster? I came up with the idea to calculate cumulative average occurrence of disasters per country per year and plot it on the graph to see the trends. If I would just calculate average occurrence of disasters per country for the whole time of the observations I would have significant issues choosing from countries that would have close numbers. Plus the total average disasters per year can be misleading by itself due to it can be high because of high amount of disasters in the beginning of XX century but relatively low number in XXI.

The formula for the calculation of the cumulative average for the given year that I used was:

Cumulative_Average = Total_Occurences / ( Given_Year - (Starting_Year - 1) ) ,

where: Total_Occurrences is the sum of occurrences of disasters for given country in time interval between the starting year and the given year (inclusive).

Here is the plot I got for the short-list countries (plotting the results for all the 180 countries from the dataset makes plot unreadable):

Cumulative average number of disasters

It is clear that China and Indonesia are the two most likely candidates for the next disaster to strike, with a China having a lead. I'm not ready to provide insight on the reasons of the increasing number of natural disasters in the countries at the plot now (especially for Turkey and Iran). Maybe it is just that the events become documented more often?... It should be investigated further.

The code

Here is the code to create the plot above. 'sqldf' package was really helpful for divide data for the short list countries from the rest of 180 countries.

library(ggplot2)
library(sqldf)
library(grid)
#library(gridExtra)


# Load natural disasters data ---------------------------------------------

dis <- read.csv("~/R/Disasters/Natural_disasters.csv")

# Create data frame with average number of disasters per year -------------

average_events <- data.frame(country = character(),
                             year = numeric(),
                             disasters_per_year = numeric(),
                             stringsAsFactors = F)

countries <- unique(dis$country)

starting_year <- min(dis$year) - 1 # we subtract 1 year to have numbers greater than 0 further on

for (country in countries) {
    data <- dis[dis$country == country,] # we need data for one country at a time
    disasters_count <- 0
    years <- unique(data$year)
    
    for (year in years) {
        total_years <- year - starting_year 
        y_data <- data[data$year == year,]
        n_disasters <- sum(y_data$occurrence)
        disasters_count <- disasters_count + n_disasters
        average_disasters <- disasters_count / total_years
        row <- data.frame(country = country, year = year, disasters_per_year = average_disasters)
        average_events <- rbind(average_events, row)
    }
    
}


# Plot data about average number of disasters per country per year --------


# Data for 180 countries is hard to plot, lets filter mots affected.
# Let's use SQL to query data: subset data for countries that had more than 0.6 disasters per year
# in any year after 2000
danger <- sqldf('SELECT * FROM average_events WHERE country IN 
      (SELECT DISTINCT country FROM average_events WHERE disasters_per_year >= 0.6 AND year > 2000)')

p <- ggplot(danger, aes (x = year, y = disasters_per_year)) + 
            geom_line(size = 1.2, aes(colour = country,  linetype = country)) +
            labs(title = 'Cumulative average number of disasters per year',
                 x = 'Year',
                 y = 'Average number of disasters cumulative') +
            guides(guide_legend(keywidth = 3, keyheight = 1)) +
            theme(axis.text.x = element_text(angle=0, hjust = NULL),
                  axis.title = element_text(face = 'bold', size = 14),
                  title = element_text(face = 'bold', size = 16),
                  legend.position = 'right',
                  legend.title = element_blank(),
                  legend.text = element_text(size = 12),
                  legend.key.width = unit(1.5, 'cm'),
                  legend.key.height = unit(1, 'cm'))
           
plot(p)

Disasters: Myth or the Reality?

noreply@blogger.com (Yury Ryabov) — Fri, 16 Jan 2015 20:56:00 +0000

I enrolled a MOOC titled "Disasters and Ecosystems: Resilience in a Changing Climate" which is organised by the UNEP (and other organisations... which names I'm going to learn by heart cause they have like 2 minutes of credits after each lecture O_o ). Not that I know nothing about disasters, risks or climate change (I'm a geographer and ecologist after all), but I was curious about the product that was made by organisation of this class.

The third video (and first video that is not an introduction) they teach us about the disasters; differences between hazard and disaster; and risks. Well... the thing they told, the graphs they showed - that what inspired the title of this post.

Terminology

Here see some definitions they use.

Disaster. When they say "disaster" they mean "natural disaster" that was enhanced by human [mismanagement].

Risk - a potential losses due to disasters.

Hazard - A dangerous phenomenon, substance, human activity or condition that may cause loss of life, injury or other health impacts, property damage, loss of livelihoods and services, social and economic disruption, or environmental damage.

Exposure - People, property, systems, or other elements present in hazard zones that are thereby subject to potential losses.

Vulnerability - the characteristics and circumstances of a community, system or asset that make it susceptible to the damaging effects of a hazard

Fails

The risk

They presented a "great" formula for (a disaster) risk evaluation that they use in the UN:

Risk = Hazard * Exposure * Vulnerability

where: Exposure = People * ExposureTime

Vulnarability - succeptability to hazard.

Well these characteristics do correspond to the risk, but the formula is stupid! I already wrote about that: Risk = Probability * Damage. And this formula actually corresponds to the definition they give (see Terminology section). We can't get a monetary outcome from their formula. We can't get numeric numeric output out of that formula at all: can you multiply flood by people? Can you???!!!

A Disaster with Disasters

The fail with the risk evaluation is a common mistake, but the fail with disaster - that is what really cool!

Take a look at this plot (which is from reading materials from the course):

What can you conclude from this plot? That the world is doing to hell and we all will fall to disaster? Let's look closer. The exposure is growing faster for poorer countries (and it is the only conclusion they make in lecture)... but the total number of people exposed (and for each type of countries) seems to be the almost unchanged! Interesting... This means (see the definition for the exposure) that there are just a 150% increase of property value in the dangerous area of the poorer countries (and 25% for the richest) on a span of 30 years. Does this graph shows us only the economic grows? I think it does... (reminds me of my previous post).

Now to the most delicious part. Take a look at this two graphs from the lecture readings:

Deaths dynamics

Damage dynamics

This is interesting. Despite the population growth and all that questionable "climate change" staff people die less (in total numbers), see fig. 1, but the damage increases, see fig. 2. Did they take inflation into account for the damage graph? Do not know... I think they didn't, otherwise they would use "discounted damage" term instead of just "damage" and would indicate the base year. So the second graph seems to demonstrate inflation and may be the economic grows.

Clearly disasters are not that disastrous. Despite the new on the TV on the subject the nature's wrath even enhanced by human is less and less dangerous for human lives. The pockets are to suffer: the storm in port wrecking the humble fisherman's boat or a trawler - that's the difference.

Conclusion

From these graphs I can conclude one thing - it is safer to live now than in the past, a disaster should not be feared as a deadly havoc. To my mind the disaster nowadays is entirely economic issue. See, if we loose less people and (maybe) more money - we should just develop more advanced insurance techniques to cover economic damage and relax. The disasters should just be studied as phenomena to develop cheap early warning systems, let the property be destroyed (just cover the losses with insurance) and additional employment to be created (rebuilding).

This is my conclusion form the graphs I showed here: disasters are an ancient myth! Just buy insurance! LOL

Do You Know What You Show at Your Map?

noreply@blogger.com (Yury Ryabov) — Sat, 10 Jan 2015 11:28:00 +0000

As access to the GIS and mapping is becoming easier every year the more people and companies create maps. Unfortunately often they just do not know what they are actually showing at their maps. This issue is being mentioned over and over again.

Here is the example that I discovered recently: Cyberthreat Real-Time Map by Kaspersky antivirus company. Here how it looks like:

Amongst the other info they show the Infection rank for each country... based on total threats detected.... You may have already guessed what is the fail, but I let me explain it anyway.

See, the №1 infected country is Russia, which is the home country for Kaspersky and where this antivirus is quite popular. So we can conclude that the rankings that supposed to demonstrate the severity of virus activities merely demonstrates the number of Kaspersky software installations across the globe.

Lets test this hypothesis. I don't have the data about the number of installation of Kaspersky software per country, but it is safe to assume that this number is proportional to the population of the given country. Also it is easier to get infection rankings for countries from the map than the number of the threats detected. If I had total threats data per country I would compare it to the population. Having infection rankings it is more rational to compare it to the population rankings instead. So I picked 27 random countries and compared their infection and population rankings. The result is demonstrated at the plot below:

Infection rank vs. Population rank

The linear model is fairly close to Inrection rank = Population rank. It is clear that the phenomena that is presented as an Infection rank just reflects a total software installations per country and not the severity of the 'cyberthreat'. In order to get the actual Infection rank the number of detected threats have to be normalised by the number of software installations.

How to Get Accurate Measurements in QGIS

noreply@blogger.com (Yury Ryabov) — Sun, 14 Dec 2014 10:35:00 +0000

This post is mainly my answer to the corresponding question.

The main rule for accurate measurements is to define correct CRS for your layers and project. This is somewhat wide theme to talk about so I won't ;-)

Lets just discuss what you should do when the project CRS can't be easily used for measurement (like Web Mercator) and you want to use Measure tool to calculate distance or area. Go Settings -> Project Properties, in CRS tab enable on the fly transformation; in General tab in Measure tool menu choose ellipsoid for distance calculations. See picture:

When you start measurements with the Measure tool move your cursor over the Measure window - an information about measurement settings will pop up just like that:

Now you can make your measurements and be sure that they are accurate. NOTE that you should always check Measure tool settings like it shown at the picture above because Ellipsoid settings for the Measure tool are not always stored or displayed in Project settings correctly (QGIS 2.6.1).

A Flawed Research on a TSP Algorithm

noreply@blogger.com (Yury Ryabov) — Sun, 30 Nov 2014 18:02:00 +0000

A Travelling Salesman Problem (TSP) is a well known computational challenge. A lot of algorithms were developed to solve it or its special cases.

I came around an article authored by Fang Liu 'A dual population parallel ant colony optimization algorithm for solving the travelling salesman problem'. In this article he proposed a modification of an Ant Colony System algorithm for solving TSP and presented results obtained by his algorithm. In the table with results all looked fine - his algorithm was able to provide very good solutions for the TSP instances from TSPLIB (which is the common testing ground for TSP algorithms).

So the researcher presented good results... it seems. But then he decided to show the best routes his algorithm was able to find and annotated them with the corresponding route costs. Lets take a look at one of them. Here you are his best route for the 'att48' instance from the TSPLIB:

Route that claims to be optimal (but the cost is very wrong)

The optimal route for 'att48' and its cost is well-known (it applies to all TSPLIB instanses). Its cost is approximately 33523 (there are different approaches to rounding distances between points). So what we see at the picture above should be the optimal route (or extremely close to it). But dear reader, do you think that you see optimal route? Humans are able to provide very good solutions to TSP instances that consists of not too many points. I bet you can draw far better route yourself. The route from this picture is 1, 8, 46, 33, 20, 17, 43, 27, 19, 37, 6, 30, 36, 28, 7, 18, 44, 31, 38, 9, 40, 15, 12, 11, 47, 21, 13, 25, 14, 23, 3, 22, 16, 41, 34, 2, 29, 5, 48, 39, 32, 24, 42, 10, 45, 35, 4, 26, 1 and its cost is 41052 which is whooping 22% far from optimal! The same story for another illustration in the article.

Here take a look at the optimal route which cost is really 33523:

Actually optimal route for 'att48' with Cost = 33523

So what we can conclude? I do believe that the routes demonstrated in the article are the best routes found by given algorithm, but costs are put-up for the routes and for the table of results in the article as well. I think that author developed an algorithm that wasn't able to find good solutions and provided fraud table with the put-up testing results. And clearly this article wasn't reviewed by scientist that have knowledge in TSP area because these plots are so obviously flawed that one can't overlook it!

No one usually shows plots of the routes they find with their algorithms. I wonder if there are more modern algorithms with the put-up results?

How to Create Delauney Triangulation Graph from a .shp-file Using NetworkX

noreply@blogger.com (Yury Ryabov) — Sun, 02 Nov 2014 17:27:00 +0000

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

Interactive Visualisation of the Profitable Amount of Waste to Dispose Illegally

noreply@blogger.com (Yury Ryabov) — Tue, 30 Sep 2014 20:57:00 +0000

"Wow!" - I said to myself after reading R Helps With Employee Churn post - "I can create interactive plots in R?!!! I have to try it out!"

I quickly came up with an idea of creating interactive plot for my simple model for assessment of the profitable ratio between the volume waste that could be illegally disposed and costs of illegal disposal [Ryabov Y. (2013) Rationale of mechanisms for the land protection from illegal dumping (an example from the St.-Petersburg and Leningrad region). Regional Researches. №1 (39), p. 49-56]. The conditions for profitable illegal dumping can be describes as follows:

Here: k - the probability of being fined for illegal disposal of waste;
P - maximum fine for illegal disposal of waste (illegal dumping);
V - volume of waste to be [illegally] disposed by the waste owner;
E - costs of illegal disposal of waste per unit;
T - official tax for waste disposal per unit.

The conditions for the profitable landfilling can be described as follows:

Here: V1 - total volume of waste that is supposed to be disposed at illegal landfill;
Tc - tax for disposal of waste at illegal landfill per unit;
P1 - maximum fine for illegal landfilling;
E1 - expenditures of the illegal landfill owner for disposal of waste per unit.

Lets plot the graphs (with some random numbers (except for fines) for a nice looking representation) to have a clue how it looks like.

Note that there is a footnote (this post provides nice examples on how to do it) with the values used for plotting - it is important to have to have this kind of indication if we want to create a series of plots.

Now I will show you the result and then will provide the code and some tips.

Playing with the plot

Tips and Tricks

Before I will show you code I want to share my hardly earned knowledge about nuances of the manipulate library. There are several ways to get static plot like that using ggplot, but some of them will fail to be interactive with manipulate.

All the data for the plot must be stored in one dataframe.
All data for plots must be derived from the dataframe (avoid passing single variables to ggplot).
Do not use geom_hline() for the horizontal line - generate values for this line and store them inside dataframe and draw as a regular graph.
To create a footnote (to know exactly which parameters were used for the current graph) use arrangeGrob() function from the gridExtra library.
Always use $ inside aes() settings to address columns of your dataframe if you want plots to be interactive

The Code

library(ggplot2)
library(grid)
library(gridExtra)
library(manipulate)
library(scales)
library(reshape2)

## Ta --- official tax for waste utilisation per tonne or cubic metre.
## k --- probability of getting fined for illegal dumping the waste owner (0 V, y <=> E

max_waste_volume <- 2000 
Illegal_dumping_fine_P <- 300000
Illigal_landfilling_fine_P1 <- 500000
Fine_probability_k <- 0.5
Official_tax_Ta <- 600

# mwv = max_waste_volume
# P = Illegal_dumping_fine_P
# P1 = Illigal_landfilling_fine_P1
# k = Fine_probability_k
# Ta = Official_tax_Ta

updateData <- function(mwv, k, P1, P, Ta){
    
    # creates and(or) updates global data frame to provide data for the plot
    
    new_data <<- NULL
    new_data <<- as.data.frame(0:mwv)
    names(new_data) <<- 'V'
    new_data$IlD <<- k*P1/new_data$V
    new_data$IlD_fill <<- new_data$IlD
    new_data$IlD_fill[new_data$IlD_fill > Ta] <<- NA # we don't want ribbon to fill area above Official tax 
    new_data$IlL <<- Ta-k*P/new_data$V
    
    new_data$Ta <<- Ta
    new_data$zero <<- 0
    dta <<- melt(new_data, id.vars="V", measure.vars=c("IlD", "IlL", "Ta"))
    dta.lower <<- melt(new_data, id.vars="V", measure.vars=c("IlD_fill", "zero", "Ta"))
    dta.upper <<- melt(new_data, id.vars="V", measure.vars=c("Ta", "IlL", "Ta"))
    dta <<- cbind(dta, lower=dta.lower$value, upper=dta.upper$value)
    dta$name <<- factor(NA, levels=c("Illegal landfill owner's\nprofitable ratio",
                                    "Waste owner's\nprofitable ratio", 
                                    "Official tax"))
    dta$name[dta$variable=="IlD"] <<- "Illegal landfill owner's\nprofitable ratio"
    dta$name[dta$variable=="IlL"] <<- "Waste owner's\nprofitable ratio"
    dta$name[dta$variable=="Ta"] <<- "Official tax"
}

updateLabels <- function(k, P1, P, Ta){
    
    ### creates footnote caption for the plot
    
    prob <- paste('Fining probability = ', k, sep = '')
    landfilling_fine <- paste('Illegal landfilling fine = ', P1, sep = '')
    dumping_fine <- paste('Illegal dumping fine = ', P, sep = '')
    tax <- paste('Official tax = ', Ta, sep = '')
    note <<- paste(prob, landfilling_fine, sep = '; ')
    note <<- paste(note, dumping_fine, sep = '; ')
    note <<- paste(note, tax, sep = '; ')
    note
}


plotDumping <- function(mwv, P, P1, k, Ta){
    
    ### this function draws the plot
    
   # initialise plot data
    updateData(mwv, k, P1, P, Ta)
    updateLabels(k, P1, P, Ta)
    
    # draw the plot
    profit <- ggplot(dta, aes(x=dta$V, y=dta$value, ymin=dta$lower, ymax=dta$upper, 
                              color=dta$name, fill=dta$name, linetype=dta$name)) +
        geom_line(size=1.2) + 
        geom_ribbon(alpha=.25, linetype=0) +
        theme(axis.text.x = element_text(angle=0, hjust = 0),
              axis.title = element_text(face = 'bold', size = 14),
              title = element_text(face = 'bold', size = 16),
              legend.position = 'right',
              legend.title = element_blank(),
              legend.text = element_text(size = 12),
              legend.key.width = unit(2, 'cm'),
              legend.key.height = unit(1.2, 'cm'))+
        scale_linetype_manual(values=c(4, 5, 1)) +
        scale_fill_manual(values = c("#F8766D","#00BFC4",NA)) +
        scale_color_manual(values = c("#F8766D","#00BFC4", '#66CD00')) + 
        labs(title="Profitable ratio between the volume \nof illegally disposed waste \nand costs of illegal disposal of waste",
             x="Waste volume, cubic meters",
             y="Cost per cubic meter, RUB") +
        xlim(c(0, max(new_data$V)))+
        ylim(c(0, Ta*1.5))
        
    
    # add a footnote about paramaters used for the current plot
    profit <- arrangeGrob(profit, 
                          sub = textGrob(note, 
                                         x = 0, 
                                         hjust = -0.1, 
                                         vjust=0.1, 
                                         gp = gpar(fontface = "italic", fontsize = 12)))
    
    # show plot
    print(profit)

}


simDumping <- function(max_waste_volume = 2000, 
                       Illegal_dumping_fine_P = 300000,
                       Illigal_landfilling_fine_P1 = 500000,
                       Fine_probability_k = 0.5,
                       Official_tax_Ta = 600) {
    
    ### this function creates interactive plot
    
    max_waste_volume <<- max_waste_volume 
    Illegal_dumping_fine_P <<- Illegal_dumping_fine_P
    Illigal_landfilling_fine_P1 <<- Illigal_landfilling_fine_P1
    Fine_probability_k <<- Fine_probability_k
    Official_tax_Ta <<- Official_tax_Ta
    
    manipulate(suppressWarnings(plotDumping(max_waste_volume, 
                           Illegal_dumping_fine_P,
                           Illigal_landfilling_fine_P1,
                           Fining_probability_k,
                           Official_tax_Ta)
                           ),
               
               # set up sliders
               max_waste_volume = slider(0, 50000, 
                                         initial = max_waste_volume, 
                                         step = 100,
                                         label = 'X axis range'),
               Illegal_dumping_fine_P = slider(0, 5000000, 
                                               initial = Illegal_dumping_fine_P,
                                               step = 10000, 
                                               label = 'Illegal dumping fine (P)'),
               Illigal_landfilling_fine_P1 = slider(0, 5000000, 
                                                    initial = Illigal_landfilling_fine_P1, 
                                                    step = 10000,
                                                    label = 'Illegal landfilling fine (P1)'),
               Fining_probability_k = slider(0, 1, 
                                             initial = 0.5,
                                             step = 0.01,
                                             label = 'Probability of being fined (k)'),
               Official_tax_Ta = slider(0, 3000, 
                                        initial = Official_tax_Ta, 
                                        step = 50,
                                        label = 'Official waste disposal tax (T)')
    )
}

simDumping() # for reasons unknown I have to run this function twice to get proper interactive plot

The Most Hilarious Rendering "Bug" I've Ever Seen

noreply@blogger.com (Yury Ryabov) — Fri, 12 Sep 2014 22:22:00 +0000

Can you guess the projection that I used to display a world map at the picture below (EPSG code is left there because I doubt someone would recall it anyway)?

It is the World Bonne projection in QGIS (original shp's CRS is EPSG:4326). I already had some fun with Bonne in this post. The world got torn apart at scale 1:75 000 000. At bigger scales map renders normally:

Fortunately this issue can be easily resolved by disabling "feature simplification by default for newly added layers" under Settings -> Options -> Rendering.

Schema of the Conservation Areas in Leningradskaya Region: Some Notes for Beginner Mappers

noreply@blogger.com (Yury Ryabov) — Sun, 07 Sep 2014 22:27:00 +0000

Schema of the Conservation Areas in Leningradskaya Region

About a year ago I was asked to create a small (a B5 size) and simple schema of the conservation areas in Leningradskaya region. I did it using QGIS. Here you are the author version of the schema and several notes that might be helpful for a beginner map-maker.

There was a huge disproportion between areas of different objects and both polygon and point markers were needed to show them. I decided to use Centroid Fill option in polygon style to be able to use only one layer (polygon) instead of two (polygon and point). Using Centroid Fill makes points in centres of the small polygons overlap and mask these tiny polygons.

All the administrative borers were stored in one layer (and there are far more borders than one see here). They are drawn using rule-based symbology so I didn't even need to subset this layer to get rid of the rest of the polygons in this layer.

All the names of the surrounding countries, regions, city and the water bodies are not labels derived from layers, but labels created inside the map composer. It was quicker and their position was more precise which is crucial for such a small schema.

There was a lack of space for the legend so I had to utilise every bit of canvas to place it. I had to use 3 legend items here. One of them is actually overlapping the other and setting a transparent background for the legends was helping with that.

Finally labels for the conservation areas (numbers) were outlined with white colour to be perfectly readable. Some of them were moved manually (with storing coordinates in special columns of the layer) to prevent overlapping with other labels and data.

P.S. Don't be afraid to argue with the client about the workflow. Initially I was asked it manually digitise a 20 000 x 15 000 pixels raster map of the Leningrad region to extract the most precise borders of the conservation areas (and districts of the region). Of course I could do it and charge more money for the job, but what's the point if some of that borders are not even to be seen at this small schema? So I convinced client to use data from OSM and saved myself a lot of time.

Donetsk Label at OpenStreetMap

noreply@blogger.com (Yury Ryabov) — Sun, 31 Aug 2014 13:32:00 +0000

Now, when Donetsk became the hot spot of the civil war at the Ukraine I decided to take a glance at the map of the war zone. Of course I used OSM for this. But despite I knew the approximate location of the city I couldn't find it for several minuets. I was confused... The reason was that "Donetsk" label (label of the city which is the administrative center of the Donetskaya region BTW) at the scales smaller than 10 km is suppressed by the label of adjusted town "Makeevka". WTF?!

"Donetsk" ("Донецьк") label at OSM is suppressed by the label of adjusted town

Random VS Pseudo Random PROC or Why Does Axe's Win Rate Dropped in 6.81 Patch (Dota 2)

noreply@blogger.com (Yury Ryabov) — Wed, 28 May 2014 20:25:00 +0000

Introduction

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.

average difference: 0.34
#######################################################################################################
|   (1) cumulative hits |   (2) random mean PROCs |   (3) pseudo random mean PROCs |   (4) difference |
|-----------------------+-------------------------+--------------------------------+------------------|
|                     1 |                    0.17 |                           0    |             0.17 |
|                     2 |                    0.45 |                           0.09 |             0.36 |
|                     3 |                    0.56 |                           0.12 |             0.44 |
|                     4 |                    0.71 |                           0.26 |             0.45 |
|                     5 |                    0.81 |                           0.36 |             0.45 |
|                     6 |                    1.06 |                           0.75 |             0.31 |
|                     7 |                    1.37 |                           0.88 |             0.49 |
|                     8 |                    1.35 |                           0.97 |             0.38 |
|                     9 |                    1.47 |                           1.22 |             0.25 |
|                    10 |                    1.64 |                           1.33 |             0.31 |
|                    11 |                    1.95 |                           1.65 |             0.3  |
|                    12 |                    2.1  |                           1.68 |             0.42 |
|                    13 |                    2.36 |                           1.79 |             0.57 |
|                    14 |                    2.56 |                           2.16 |             0.4  |
|                    15 |                    2.61 |                           2.26 |             0.35 |
|                    16 |                    2.61 |                           2.37 |             0.24 |
|                    17 |                    2.87 |                           2.43 |             0.44 |
|                    18 |                    2.77 |                           2.83 |            -0.06 |
|                    19 |                    3.22 |                           2.84 |             0.38 |
|                    20 |                    3.03 |                           2.84 |             0.19 |
|                    21 |                    3.6  |                           3.22 |             0.38 |
|                    22 |                    3.71 |                           3.32 |             0.39 |
|                    23 |                    3.68 |                           3.68 |             0    |
|                    24 |                    3.93 |                           3.72 |             0.21 |
|                    25 |                    4.54 |                           3.92 |             0.62 |
|                    26 |                    4.65 |                           4.04 |             0.61 |
|                    27 |                    4.47 |                           4.27 |             0.2  |
|                    28 |                    4.83 |                           4.39 |             0.44 |
|                    29 |                    4.78 |                           4.48 |             0.3  |
|                    30 |                    4.93 |                           4.8  |             0.13 |
|                    31 |                    5.3  |                           4.92 |             0.38 |
|                    32 |                    5.1  |                           5.04 |             0.06 |
|                    33 |                    5.79 |                           5.34 |             0.45 |
|                    34 |                    5.82 |                           5.54 |             0.28 |
|                    35 |                    6.04 |                           5.52 |             0.52 |
|                    36 |                    5.67 |                           5.7  |            -0.03 |
|                    37 |                    6.64 |                           5.98 |             0.66 |
|                    38 |                    6.4  |                           6.03 |             0.37 |
|                    39 |                    6.72 |                           6.41 |             0.31 |
|                    40 |                    7.07 |                           6.46 |             0.61 |
|                    41 |                    7.04 |                           6.59 |             0.45 |
|                    42 |                    7.18 |                           6.81 |             0.37 |
|                    43 |                    7.08 |                           6.9  |             0.18 |
|                    44 |                    7.61 |                           7.09 |             0.52 |
|                    45 |                    7.72 |                           7.21 |             0.51 |
|                    46 |                    7.73 |                           7.66 |             0.07 |
|                    47 |                    8    |                           7.68 |             0.32 |
|                    48 |                    8.41 |                           7.7  |             0.71 |
|                    49 |                    8.57 |                           7.93 |             0.64 |
|                    50 |                    8.54 |                           8.11 |             0.43 |
|                    51 |                    8.16 |                           8.31 |            -0.15 |
|                    52 |                    9    |                           8.4  |             0.6  |
|                    53 |                    9.01 |                           8.79 |             0.22 |
|                    54 |                    8.98 |                           8.88 |             0.1  |
|                    55 |                    9.54 |                           8.99 |             0.55 |
|                    56 |                    9.4  |                           9.13 |             0.27 |
|                    57 |                    9.75 |                           9.08 |             0.67 |
|                    58 |                   10.1  |                           9.42 |             0.68 |
|                    59 |                    9.71 |                           9.64 |             0.07 |
|                    60 |                   10.31 |                           9.87 |             0.44 |
|                    61 |                   10.19 |                          10.31 |            -0.12 |
|                    62 |                   10.29 |                          10.21 |             0.08 |
|                    63 |                   10.76 |                          10.55 |             0.21 |
|                    64 |                   10.82 |                          10.48 |             0.34 |
|                    65 |                   10.7  |                          10.77 |            -0.07 |
|                    66 |                   11.27 |                          10.94 |             0.33 |
|                    67 |                   11.81 |                          11.06 |             0.75 |
|                    68 |                   11.54 |                          11.34 |             0.2  |
|                    69 |                   11.98 |                          11.26 |             0.72 |
|                    70 |                   12.26 |                          11.66 |             0.6  |
|                    71 |                   11.35 |                          12.01 |            -0.66 |
|                    72 |                   12.03 |                          11.64 |             0.39 |
|                    73 |                   12.37 |                          11.94 |             0.43 |
|                    74 |                   12.74 |                          12.16 |             0.58 |
|                    75 |                   13.16 |                          12.66 |             0.5  |
|                    76 |                   12.77 |                          12.59 |             0.18 |
|                    77 |                   13.1  |                          12.65 |             0.45 |
|                    78 |                   13.2  |                          12.95 |             0.25 |
|                    79 |                   13.59 |                          13.06 |             0.53 |
|                    80 |                   13.32 |                          13.27 |             0.05 |
|                    81 |                   13.74 |                          13.25 |             0.49 |
|                    82 |                   13.98 |                          13.47 |             0.51 |
|                    83 |                   14.86 |                          13.74 |             1.12 |
|                    84 |                   14.53 |                          13.8  |             0.73 |
|                    85 |                   14.54 |                          14.29 |             0.25 |
|                    86 |                   14.49 |                          14.22 |             0.27 |
|                    87 |                   15.03 |                          14.46 |             0.57 |
|                    88 |                   15.49 |                          14.55 |             0.94 |
|                    89 |                   15.51 |                          14.8  |             0.71 |
|                    90 |                   15.58 |                          15    |             0.58 |
|                    91 |                   16.17 |                          15.01 |             1.16 |
|                    92 |                   14.89 |                          15.43 |            -0.54 |
|                    93 |                   15.73 |                          15.4  |             0.33 |
|                    94 |                   16.16 |                          15.67 |             0.49 |
|                    95 |                   16.11 |                          16.17 |            -0.06 |
|                    96 |                   16.03 |                          16.17 |            -0.14 |
|                    97 |                   16.41 |                          16.07 |             0.34 |
|                    98 |                   17.26 |                          16.27 |             0.99 |
|                    99 |                   15.65 |                          16.7  |            -1.05 |
|                   100 |                   16.15 |                          16.96 |            -0.81 |

Classification of the Hyper-Spectral and LiDAR Imagery using R (mostly). Part 3: Shadow Removal

noreply@blogger.com (Yury Ryabov) — Thu, 08 May 2014 17:02:00 +0000

This is my third post in this series, here are the links to the Part 1 and Part 2. Sorry for the delay, but at that time I got completely consumed with my PhD thesis and this story just went out of my mind.

To the point. There were 2 types of shadows I had to deal with: cloud shadows and cast shadows (from building and trees). In scientific literature you can find several relatively complicated algorithms to deal with them (including DEM usage for cast shadow detection) and I definitely should have used some of them, but... I was too lazy to implement them when there were couple of quick'n'dirty workarounds available.

Cloud shadow removal

To remove cloud shadows I decided to calculate ratio between mean pixel value of the lighted up and shadowed parts of the raster for each band and multiply shadowed pixels by this ratio. I manually created a mask for shadowed areas. There were 3 areas and their borders were quite distinct in NIR channels:

Cloud shadows in one of the NIR bands

I drew several polygons over shadows in QGIS and converted them to raster (Raster -> Conversion -> Rasterise). Now everything was ready for patching:

library(raster)
library(rgdal)

# load rasters
image <- stack('2013_IEEE_GRSS_DF_Contest_CASI.tif') # hyper spectral imagery
s_mask <- raster('shadow_mask.tif') # mask for cloud shadow

# extract shadowed area
m_cloud <- s_mask
m_cloud[m_cloud < 1] <- NA
shadowed <- mask(image, m_cloud)

# calculate zonal statistics
zd <- zonal(image, s_mask, fun = 'mean', na.rm = T)
zd <- as.data.frame(zd)

# calculate ratio between lighted up and shadowed zones
zd[3,] <- zd[1,] / zd[2,] 

# recalculae data in shadows
multiplyer <- zd[3,]
multiplyer <- multiplyer[,2:ncol(zd)]
multiplyer <- as.numeric(multiplyer)
enlight <- shadowed*multiplyer # patch layer
restored <- cover(enlight, image, progress = 'text')

# save result
wr <- writeRaster(restored, 
                  filename = 'restored.tif', 
                  datatype = 'GTiff', overwrite = TRUE)

Now it's time to check out results:

The same NIR band after cloud shadow removal

Yes, we can easily find patched areas, but overall result is good for such simple approach.

Cast shadow removal

Shadows from buildings can be processed the same way as cloud shadows as shown above except mask creation approach is different. In case of cast shadows they have to be detected automatically (for there are a bit more than 3 as in case with clouds). Using one of the NIR bands I created a learning sample of points and passed it to Random Forest algorithm to classify imagery into shadowed and non-shadowed areas. I will not describe here this classification process for Random Forest was used for overall classification of the imagery and it is the subject of my next post of this series.