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.
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()
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()
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()
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