How to use a raster dataset to store and retrieve "vectors" using Python?

Question

For a rectangular area of about 10 by 10 km, I have to generate "displacement vectors" with a 1 meter resolution (i.e. a delta-x and delta-y value for every cell of 1 by 1 meter). I will end up with 100 million tuples (x, y, delta_x, delta_y). In a plain text file, this would look like:

x y delta_x delta_y
------ ------ ------- -------
228682 184554 0.429 1.221 
228682 184555 0.428 1.222 
228682 184556 0.428 1.223 
228682 184557 0.427 1.224 
228682 184558 0.427 1.225 
...
229526 184822 -0.084 0.593 
229526 184823 0.185 0.002 
229526 184824 0.185 0.001 
229526 184825 0.186 -0.001 
229526 184826 0.186 -0.002 
...

I was wondering if there is an efficient way to store these values as a dataset of type raster/coverage...? The x, y would be the pixel or cell position, and the delta_x and delta_y would be the values. Is HDF5 the way to go? I have no experience using HDF5.

I need to store the values in Python and retrieve them in Python in an efficient way (more efficient then first reading 100 million tuples from a text file, then storing them in a R-Tree index and finally using the index doing intersects).

I am looking for some guidance, tips, best practices, small examples or use cases where this has been done before.

Answer 1

Your first chance at reducing the size of the data is to remove the need to explicitly store x and y. Since your data is regularly gridded you can implicitly store these values as the grid indices.

This way you can generate a 3 dimensional numpy array with the shape (2, n_x, n_y).

delta_x and delta_y for the location x = 20 and y = 10 would then be stored at array(:, 10, 20).

There are a number of options to store such an array on disk - numpy.savez and h5py probably being the most popular ones.

import numpy as np
import h5py
test_data = np.arange(200000000).reshape(2, 10000, 10000).astype(np.float32)
np.savez("/data/test_tmp/numpy.npz", test_data)
np.savez_compressed("/data/test_tmp/numpy_compressed.npz", test_data)
with h5py.File("/data/test_tmp/h5py.h5", "w") as hdf5:
 hdf5.create_dataset("innit", data=test_data, compression="gzip")

>ls -al /data/test_tmp/
-rw-rw-r-- 1 dl dl 77886268 Jan 13 17:31 h5py.h5
-rw-rw-r-- 1 dl dl 72289983 Jan 13 17:31 numpy_compressed.npz
-rw-rw-r-- 1 dl dl 800000212 Jan 13 17:30 numpy.npz

On this test dataset you can see the effects of compression. The uncompressed array numpy.npz has a size of ~800 Mb while both compressed arrays are around a 10th of the size at ~75 Mb.

I prefer numpy.savez for its ease of use but h5py definitely gives you a lot more options in regards of compression and multiple datasets in a container.

• user52724: your solution seems perfect, and thank you for the elaborate answer. In the mean time, I was working on a similar solution using Numpy and storing the values in a GeoTiff having 2 bands, one for delta_x and one for delta_y. The resulting uncompressed GeoTiff is about 10 % the size of my text version (storing delta's in centimeters, hence integers, not floats). I still have to do the reading part. What's your opinion regarding my solution compared to yours? In particular with regard to reading the values.
• Kersten: Your solution basically achieves the same results, although I would say with unnecessary complication (you don't really need GDAL since you don't really need the spatial information for storage - or do you?) Can you provide a subset? Then we would be able to do a performance comparison.
• user52724: Your solution is marked as the answer because it was very to the point and more generic than what I am doing. However, I will go with GeoTiff solution. It's working now and fast. Moreover, this was not mentioned but a colleague of mine needs to use the dataset in Smallworld coding in Magic. I have no idea what his capabilities are, but my guess is that consulting the value of a GeoTiff is a safer option than your output that he might not be able to read. But thanks a bunch for the help!

Answer 2

The alternative answer using a GeoTiff for storing the delta_x and delta_y values uses Numpy and GDAL. Here is the outline for the solution.

The implementation for creating the GeoTiff uses the code in this recipe:

https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#create-raster-from-array

The delta_x and delta_y values are stored in separate bands. I chose to store the displacements in centimeters to reduce the size of the GeoTiff (signed 16-bit integers).

I read the values from the GeoTiff in Python using the code in this recipe:

https://gist.github.com/stefanocudini/5201689

The GeoTiff solution is less generic having a dependency on GDAL, but this comes in handy if the output need to be read in other GIS.

• Kersten: +1, this is indeed a nice solution if someone needs to load the data into a GIS or related geospatial software.