Date created: June 10, 2022
Author: P. Alexander Burnham

Summary Last time we set up several cubes from geotif files and explored some of the basic querrying and cube building functionality of cubetime. However, all of the files we were working with, had the same scales and grid sizes. Spacetime has the functionality to rescale rasters of different extents, grid sizes, epsg codes etc. and create a cleaned cube where all of those factors are consistent accross the data layers.

Loading file names and creating a file object

Let’s load two rasters of data from India. We can then load the data up as a file object and compare some attributes between them.

# read files and load data
data = ["demoData/LULC_1995.tif", "demoData/India_cropped-area_1km_2016.tif"]
ds = read_data(data)
# check number of bands
ds.get_band_number()

# check the epsg codes
## [1, 1]
ds.get_epsg_code()

# check the raster dimensions
## ['EPSG:32644', 'EPSG:4326']
ds.get_dims()

# check pixel size
## [(29484, 33555), (2373, 3269)]
ds.get_pixel_size()

# check no data values
## [100.0, 0.008868160000000002]
ds.get_nodata_value()
## [127.0, -1.7e+308]

Let’s plot the two files and compare them.

import matplotlib.pyplot as plt

# Get the array data
arrays = ds.get_data_array()

# plot file 1
plt.imshow(arrays[0], vmin=0, vmax=17)
plt.show()

# plot file 2

plt.imshow(arrays[1], vmin=0, vmax=100)
plt.show()

Let’s rescale these files using some spacetime functions.

Apart from the fact that both datasets are from approximately the same region, we can see that almost all of the attributes between these two files are different. Using raster_scale() we can get these two files to have the same spatial reference system,, pixel size and nodata value with cells aligned to the same lat and lon values.

# lets align these files
alignedDS = raster_align(data=ds, resolution=.08, SRS=4326, noneVal=127)

Let’s look at these attributes now that we have run the align function.


# check epsg codes
alignedDS.get_epsg_code()
# get raster dimensions
## ['EPSG:4326', 'EPSG:4326']
alignedDS.get_dims()
# pixel size
## [(409, 383), (264, 363)]
alignedDS.get_pixel_size()
# check no data values
## [0.08, 0.08]
alignedDS.get_nodata_value()
## [127.0, 127.0]

Lets also look at the rescaled images.

# get rescaled data
scaledArray = alignedDS.get_data_array()

# plot file 1
plt.imshow(scaledArray[0], vmin=0, vmax=17)
plt.show()

# plot file 2

plt.imshow(scaledArray[1], vmin=0, vmax=100)
plt.show()

Cropping file objects

We can see now that the only difference is the extent of the files. The first file covers more surface area of the globe. Lets crop the files to have the same greatest common bounding box.

# trim the file object
trimmedDS = raster_trim(alignedDS)
# lets compare the dimensions of the new objects
trimmedDS.get_dims()
## [(264, 363), (264, 363)]

Comparing the plotted images reveals that the files are now propperly aligned and trimmed

# get data
trimmedArray = trimmedDS.get_data_array()

# plot file 1
plt.imshow(trimmedArray[0], vmin=0, vmax=17)
plt.show()

# plot file 2

plt.imshow(trimmedArray[1], vmin=0, vmax=100)
plt.show()

Make a cube

Now that we have cleaned data sets, we can assemble them into a cube as shown in the previous vignette. Lets start be creating a time object. cube_time() allows you to create a time vector that skips user-specified chunks of time using the skips arguemnt.

# create time vec
yearObj = cube_time(start="1995", length=2, scale = "year", skips = 21)

We will now create a cube by assigning files and bands to time.

cube = make_cube(data = trimmedDS, fileName = "indiaCube.nc4", organizeFiles = "filestotime", organizeBands="bandstotime", timeObj = yearObj)

Finally, lets look at the data

cube.get_data_array()
<xarray.DataArray (lat: 363, lon: 264, time: 2)>
array([[[127., 127.],
        [127., 127.],
        [127., 127.],
        ...,
        [127., 127.],
        [127., 127.],
        [127., 127.]],

       [[127., 127.],
        [127., 127.],
        [127., 127.],
        ...,
        [127., 127.],
        [127., 127.],
        [127., 127.]],

       [[127., 127.],
        [127., 127.],
        [127., 127.],
        ...,
...
        ...,
        [127., 127.],
        [127., 127.],
        [127., 127.]],

       [[127., 127.],
        [127., 127.],
        [127., 127.],
        ...,
        [127., 127.],
        [127., 127.],
        [127., 127.]],

       [[127., 127.],
        [127., 127.],
        [127., 127.],
        ...,
        [127., 127.],
        [127., 127.],
        [127., 127.]]], dtype=float32)
Coordinates:
  * lon      (lon) float32 68.08 68.16 68.24 68.32 ... 88.88 88.96 89.04 89.12
  * lat      (lat) float32 37.12 37.04 36.96 36.88 36.8 ... 8.4 8.32 8.24 8.16
  * time     (time) datetime64[ns] 1995-12-31 2016-12-31