Georeferencing a numpy tensor (spatio-temporal data)

Question

I have a CSV file in which each row represents a record of an event with the following columns:

latitude, longitude, time

I want to build a numpy tensor with n x m x 8760 (hours in a year) dimensions where each cell will represent the sum of the events in the corresponding space and time.

So, in order to correlate the CSV events with the numpy tensor, I need to georeference the numpy tensor in a spatial reference system (ideally WGS 84).

Until now I have defined the numpy tensor as follows:

import numpy as np
from numpy import array
# define a spatial grid with time dimension (8760 hours)
c = np.zeros((8760, 100, 100)).astype(int) + 1

My first five lines from my CSV file are presented below:

Longitude, Latitude, Date, Time
-1.60202910794,53.7518017092,3/17/2017,0815 
-1.53349564433,53.7944009244,1/14/2017,1330
-1.5624701361,53.7642562523,1/1/2017,0805
-1.5624701361,53.7642562523,1/1/2017,0805

So how can I georeference a numpy tensor with n x m x 8760 ?

Answer 1

I hope I understood your problem correctly (although still don't know why you added timestamps at the end if your csv clearly shows the time being two string columns).

First of all, you have to define an origin in space for your tensor. This origin consists of a lon, lat pair of coordinates that represent the upper left corner of your array. Let's assume this origin is -1.712604, 53.891014, just northwest of Leeds, UK. Now, you have to define the size of each cell of your array. Let's assume it is 0.005° in x and 0.005° in y (you'll have square pixels where each side is roughly 500m). The following image shows a grid representig each one of your cells (without the time dimension) and your four points (one is a duplicate so there are only 3 points visible) in red.

What you want to accomplish can be outlined in these two steps:

• 1) For each point, get the index of its location in your array. Because you have three dimensions in your array (time, y, x) you'll need to get a three-value index for each point.
• 2) For each index, get the cell of the array and add 1.

The first step is the complex one. To get the index of each point, you need to find what hour from the 8760 in the year your date belongs to. Then, you have to find what row and column your latitude and longitude values belong to. Here is a function that accomplish this, leveraging numpy vectorization. This means this functions is called only once and get the index of every point in your csv file instead of requiring you to do a for loop for each point.

def get_indices(t, time, x, y, ox, oy, pw, ph):
    """
    Gets the band (k), row (i) and column (j) indices in an array for a
    given set of timestamps and coordinates. Partly based on
    https://gis.stackexchange.com/a/92015/86131

    :param t:   array of datetime values
    :param t:   array with a range of datetime values
    :param x:   array of x coordinates (longitude)
    :param y:   array of y coordinates (latitude)
    :param ox:  raster x origin (left boundary)
    :param oy:  raster y origin (upper boundary)
    :param pw:  raster pixel width
    :param ph:  raster pixel height
    :return:    band (k), row (i) and column (j) indices
    """
    k = np.searchsorted(time, t)
    i = np.floor((oy-y) / ph).astype('int')
    j = np.floor((x-ox) / pw).astype('int')

    return k, i, j

So far, so good. Now you have to get the values you are going to pass to this function in order to get all the indices. First, let's declare the values we mentioned above:

# specify array origin and pixel resolution
ox = -1.712604
oy = 53.891014
pw = 0.005
ph = 0.005

Now, lets create an array with a datetime value for each hour in the year (using pandas):

# create an array with 8760 datetime values (one for each hour)
time = pd.date_range('2017-01-01', '2018-01-01', freq='H', closed='right')
time = pd.Series(time).values

You still need to get your actual values: the ones stored in the csv. I created a test.csv file with the table you provided (and removed the white spaces it had). Here is a snippet to read the csv as a pandas DataFrame, create a new Datetime column with actual datetime values and then getting all the needeed values (time, longitude and latitude) into numpy arrays. Note that you have to use the dt.floor('H') method on your datetime values so they are closed hours (e.g 08:00 instead of 08:15) and you can match the dates created above with pd.date_range()

# read data and create a datetime column
df = pd.read_csv('test.csv', dtype={'Time': str})
df['Datetime'] = df['Date'] + ' ' + df['Time']
df['Datetime'] = pd.to_datetime(df['Datetime'], format='%m/%d/%Y %H%M')
df['Datetime'] = df['Datetime'].dt.floor('H')

# get values as numpy arrays
t = df['Datetime'].values
x = df['Longitude'].values
y = df['Latitude'].values

Now it's time to call the function and get the indices with the following line:

idx = get_indices(t, time, x, y, ox, oy, pw, ph)

If you inspect the contents of idx, you'll get the following tuple of 1D arrays:

(array([1807,  324,    7,    7]),
 array([27, 19, 25, 25]),
 array([22, 35, 30, 30]))

For example, your first point corresponds to the 1807^th hour (band), the 27^th row and the 22^nd column in your tensor.

Finally, you have to create an array full of zeros, and add 1 to each index in your tensor (leveraging again, numpy vectorization). This operation will take into account repeated indices, as you can see from the result of summing all the cells in your tensor:

arr = np.zeros((8760, 100, 100), dtype=int)
np.add.at(arr, idx, 1)
arr.sum()  # 4

Note: the get_indices() function will only work for points that lie within your defined grid. If they are outside you'll get either negative indices or indices out of bounds.