Create a Raster from Array

Create a Raster from Array#

Introduction#

It is sometimes useful to create a test image of known pixel values for testing geospatial algorithms. Other times, you may have an array of pixels from a non-spatial packag that needs to be turned into a georeferenced raster. In this tutorial, we will see how XArray and rioxarray can be used to accomplish this.

Overview of the Task#

We use xarray to create a DataArray from an array of pixels and use rioxarray to assign a CRS and create a GeoTIFF image.

Inputs:

  • An array of pixel values, coordinates of the upper left pixel, pixel resolution and the CRS.

Output:

  • image.tif : A georeferenced image with the supplied pixel values.

Setup and Data Download#

The following blocks of code will install the required packages and download the datasets to your Colab environment.

%%capture
if 'google.colab' in str(get_ipython()):
    !pip install rioxarray

import os
import numpy as np
import xarray as xr
import rioxarray as rxr

output_folder = 'output'

if not os.path.exists(output_folder):
    os.mkdir(output_folder)

We want to create an image with a resolution of 100 meters with the upper-left pixel at coordinates (780850,1432187) in the CRS EPSG:32643 (WGS84 / UTM Zone 43N).

upper_left_x, upper_left_y = (780850,1432187)
resolution = 1000 # 1 km
crs = 'EPSG:32643'

If we wanted the resulting image with known pixel values, we can define a 2-dimentional array. Since we are storing small integers values set the data type to Byte (uint8). If you had larger integers or floating point numbers, you can use appropriate data type.

array = np.array([
    [0, 0, 1, 1],
    [0, 0, 1, 1],
    [0, 2, 2, 2],
    [2, 2, 3, 3]
], dtype=np.uint8)

Another option is to create an array of random values. The following block create an array of 100 x 100 pixels with random values between 0 to 4.

array = np.random.randint( 0, 4, size=(100, 100)).astype(np.uint8)

Next, we need to assign X and Y coordinates for each pixels. We use np.linspace function to create a sequence of x and y coordinates for each pixel of the image.

num_pixels = array.shape[0]

x_coords = np.linspace(
    start=upper_left_x,
    stop=upper_left_x + (resolution*(num_pixels-1)),
    num=num_pixels, dtype=np.uint)
y_coords = np.linspace(
    start=upper_left_y,
    stop=upper_left_y - (resolution*(num_pixels-1)),
    num=num_pixels, dtype=np.uint)
x_coords, y_coords

(array([780850, 781850, 782850, 783850], dtype=uint64),
 array([1432187, 1431187, 1430187, 1429187], dtype=uint64))

Now we create a DataArray and assign x and y coordinates.

da = xr.DataArray(
    data=array,
    coords={
        'y': y_coords,
        'x': x_coords
    }
)
da

<xarray.DataArray (y: 4, x: 4)> Size: 16B
array([[0, 0, 1, 1],
       [0, 0, 1, 1],
       [0, 2, 2, 2],
       [2, 2, 3, 3]], dtype=uint8)
Coordinates:
  * y        (y) uint64 32B 1432187 1431187 1430187 1429187
  * x        (x) uint64 32B 780850 781850 782850 783850

Next, we assign a CRS. The rioxarray extension provides a rio accessor that allows us to set a CRS.

Even though we are not using rioxarray directly, we still need to import it which activates the rio accessor in xarray.

da = da.rio.write_crs(crs)
da

<xarray.DataArray (y: 4, x: 4)> Size: 16B
array([[0, 0, 1, 1],
       [0, 0, 1, 1],
       [0, 2, 2, 2],
       [2, 2, 3, 3]], dtype=uint8)
Coordinates:
  * y            (y) uint64 32B 1432187 1431187 1430187 1429187
  * x            (x) uint64 32B 780850 781850 782850 783850
    spatial_ref  int64 8B 0

Now we can save the DataArray in any of the supported raster formats.

output_file = 'image.tif'
output_path = os.path.join(output_folder, output_file)
da.rio.to_raster(output_path)

If you want to give feedback or share your experience with this tutorial, please comment below. (requires GitHub account)