G2S: The GeoStatistical Server

A free and flexible multiple point (geo)statistics framework including state-of-the-art algorithms:
QuickSampling and Narrow Distribution Selection

View project on GitHub

Brief overview

This page provides some additional information on how to do common tasks. It also serves as personal note to myself 😁

Multiple realization

A question that is frequently asked is how to do multiple realizations. Currently there are three different ways to do it, each one of them has pros and cons.

The lazy solution

The simplest solution is to do a for-loop over the realizations. However, in each step the algorithm needs to wait for the data to load, and this creates an overhead


sims1=numpy.empty((250,250,numberOfSimulation));
for i in range(numberOfSimulation):
	sims1[:,:,i],_=g2s('-sa',computationServer,'-a','qs','-ti',ti,'-di',numpy.zeros((250,250))*numpy.nan,'-dt',numpy.ones((1,)),'-k',1.2,'-n',50,'-j',0.5,'-s',100+i);
				

For-loop without overhead

By using the G2S submission queue, it's possible to remove most of the overhead. This is the most versatile solution and it is recommended over the first solution, especially in case the computations are run on another machine.
Furthermore, in this solution the parameters can be changed for each realization.


ids=numpy.empty((numberOfSimulation,), dtype='long');
for i in range(numberOfSimulation):
	ids[i]=g2s('-sa',computationServer,'-a','qs','-ti',ti,'-di',numpy.zeros((250,250))*numpy.nan,'-dt',numpy.ones((1,)),'-k',1.2,'-n',50,'-j',0.5,'-s',100+i,'-submitOnly');

sims2=numpy.empty((250,250,numberOfSimulation));
for i in range(numberOfSimulation):
	sims2[:,:,i],_=g2s('-sa',computationServer,'-waitAndDownload',ids[i]);
				

Overdimensioning of the destination image

This is the third solution to get multiple simulations at once. Although this solution looks easier, it has more limitations (e.g., being limited same size, same parameters, etc.), and therefore it is ⚠️ not guaranteed to stay functional in the future. This last solution should only be considered in case of extremely parallelized simulations (number of threads >50) and/or extremely small simulations (less than 1e4 pixels per simulation)


sims3,_=g2s('-sa',computationServer,'-a','qs','-ti',ti,'-di',numpy.zeros((numberOfSimulation,250,250))*numpy.nan,'-dt',numpy.ones((1,)),'-k',1.2,'-n',50,'-j',0.5);
				

How to do a simulation by segment

Although QS tends to reduce issues related to non-stationarity, it won't remove them completely. Therefore in cases with high non-stationarity, the challenge is to respect each zone separately. Here we assume that we have different training images for each specific zone.

Using a secondary variable

The first solution for doing a simulation by segment is to add a secondary variable with the information of non-stationarity. In this case the setting of the algorithm has to be identical for both training images.


ti1=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/bangladesh.tiff');
ti2=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff');

dt=[0];
N=300;

ti1_bis=cat(3,ti1,zeros(size(ti1)));
ti2_bis=cat(3,ti2,ones(size(ti2)));

localKernel=cat(3,kernel,padarray(15,25*[1,1]));
% localKernel=cat(3,kernel,kernel);

sim=g2s('-a','qs','-ti',ti1_bis,ti2_bis,'-di',cat(3,nan(size(tiSelection)),tiSelection),'-dt',[0,1],'-n',50,'-k',1.5,'-j',0.5,'-ki',localKernel);
imshow(medfilt2(sim(:,:,1)))
drawnow;
pause(5)

Using training image indices

The second solution for doing a simulation by segment is to use a training image index map. This requires that each training image represents a stationary subdomain. Also in this case the setting of the algorithm needs to be identical between training images.


ti1=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/bangladesh.tiff');
ti2=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff');

dt=[0];
N=300;

sim=g2s('-a','qs','-ti',ti1,ti2,'-di',nan(size(tiSelection)),'-dt',dt,'-n',50,'-k',1.5,'-j',0.5,'-ki',kernel,'-ii',tiSelection);
imshow(medfilt2(sim))
pause(5)

Simulating each subdomain sequentially

The last solution for doing a simulation by segment is to do simulations per subdomain. This allows for changing the settings between each subdomain. However each subdomain is simulated sequentialy and therefore can be at the origin of some artifacts, regarding the order of simulations.


ti1=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/bangladesh.tiff');
ti2=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff');

dt=[0];
N=300;

path=reshape(randperm(numel(tiSelection)),size(tiSelection));
p1=path;
p1(tiSelection==0)=-inf;
p2=path;
p2(tiSelection==1)=-inf;
sim=g2s('-a','qs','-ti',ti1,'-di',nan(size(tiSelection)),'-dt',dt,'-n',50,'-k',1.5,'-j',0.5,'-ki',kernel,'-sp',p1);
sim=g2s('-a','qs','-ti',ti2,'-di',sim,'-dt',dt,'-n',50,'-k',1.5,'-j',0.5,'-ki',kernel,'-sp',p2);
imshow(medfilt2(sim))
pause(5)