Monday 24 July 2017

Projecting coordinates in Python

Sometimes we have to manipulate huge spatial databases, and it uses to be a headache if we must to convert the spatial point to a different projection. Well, an excellent way to solve this sort of issues is taking advantage of the libraries of python. In this case, I am going to describe how to use the package utm for it.

First, we are going to simulate some spatial data in a decimal projection.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import utm

lat = np.linspace(20, 40, 100)
long = np.linspace(60, 70, 100)
index = np.arange(1,101)

d = {'x': pd.Series(lat, index=index),
     'y': pd.Series(long, index=index)}

df = pd.DataFrame(d)
print(df)

Done!, we have already created a Dataframe with one hundred of spatial points in decimal degrees. Now, we are going to project them to UTM. First, we make a temporal "temp" variable which will store our projected XY coordinates. Afterwards, we are going to use the function utm.from_latlon() in a loop. This instruction will put each project coordinate in your variable "temp" that you created recently.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import utm

lat = np.linspace(20, 40, 100)
long = np.linspace(60, 70, 100)
index = np.arange(1,101)

d = {'x': pd.Series(lat, index=index),
     'y': pd.Series(long, index=index)}

df = pd.DataFrame(d)
#print(df)

temp = [[None]*1 for ii in range(len(df.x))]

## Use the utm package to transform coordinates
for i in np.arange(0, len(df.x)):
    temp[i] = utm.from_latlon(df.y[i+1], df.x[i+1])

# Transform in a dataframe
data = pd.DataFrame(temp)
data.columns = ['X', 'Y', 'Zone', 'Code']
print(data)

Finally, we convert the projected data to a DataFrame and assign this to a new variable "data", and that's all. This procedure helps me to transform a database with two hundred thousand locations, so it is useful to project massive databases.





Thursday 20 July 2017

Diffusion equation

The diffusion equation is widely discussed in the literature, especially in physics and mathematics. However, the diffusion equation also has relevance in other fields, for instance, biology, where it can be useful to understand some biological phenomena at different scales, including population biology. Therefore, it is crucial to know how to solve it.

The diffusion equation is a partial differential equation which involves both a time and space. In 1D it would be:

$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}$

It can be discretised using numerical methods. In this case, we are going to use the finite difference method, so that the above equation can be discretised as follows:

$\frac{u^{t+1}_{i} - u^{t+1}_{i}}{\Delta t} = D \frac{u^t_{i+1} - 2u^t_{i} + u^t_{i-1}}{\Delta^2x}$

where, $t$ indicate the time set and $i$ the x position.

To obtain the steady-state, we are going to programming this in python. The below code solve the diffusion equation numerically
import numpy as np
import matplotlib.pyplot as plt

#Setting parameters
nx = 101
dx = 0.01
nt = 20000
D = .1
dt = 0.0001 #be aware setting this value

#Initial conditions
u = np.zeros(nx)
u[25:75] = 2
un = u.copy()

#Numerical scheme: Finite differences
for n in range(nt):
    for i in range(1, nx-1):
        un = u.copy()
        u[i] = un[i] + D * dt /dx**2* (un[i+1]-2*un[i]+un[i-1])

#Plot the results
plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()

The solution is



Where it is shown how the distribution diffuses across the domain; however, as time passes, the distribution density tends to zero, so that the mass conservative property is not satisfied. If we wish to maintain a constant mass, we have to apply boundary conditions, in this case, we need to set-up a zero flux condition in the extremes of the domain. Then:

$\frac{\partial u}{ \partial t} = 0 $


import numpy as np
import matplotlib.pyplot as plt

#Setting parameters
nx = 101
dx = 0.01
nt = 20000
D = .1
dt = 0.0001 #be aware setting this value

#Initial conditions
u = np.zeros(nx)
u[25:75] = 2
un = u.copy()

#Numerical scheme: Finite differences
for n in range(nt):
    for i in range(1, nx-1):
        un = u.copy()
        u[i] = un[i] + D * dt /dx**2* (un[i+1]-2*un[i]+un[i-1])                                                                                                       # Zero flux condition
        u[0] = u[1]
        u[-1]= u[-2]                                                            
#Plot the results
plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()

Here, zero flux condition is included straightforward by adding:  u[0] = u[1]   (left)  and
u[-1]= u[-2] (right). Then, the mass can not go outside of the domain.

But!, be careful, because this zero flux condition only works for this simple equation. If you are working with different models, for instance, the advection-diffusion equation, zero flux condition has to be implemented differently.

Sunday 16 July 2017

Updating spyder

If you are python users, maybe you use spyder as your IDLE. For update your spyder version, the common procedure is by the instruction:
conda update spyder

The above should work if you have conda package manager. Moreover, you can also upgrade spyder for specific environments, for instance, if we have a python environment called py27, we can upgrade spyder with:

conda update -n py27 spyder

However, whether the above does not work, another way is using the pip instruction:

pip install --upgrade spyder