Expressions
The power of Yirgacheffe comes from being able to operate on geospatial data as if it was an elemental type. You can combine and work with layers without worrying about individual pixels or how different layers are aligned spatially. Assuming your data is in the same projection and pixel scale, you can just get on working. For example, here's how a simple Area of Habitat calculation might be implemented:
import yirgaceffe as yg
with (
yg.read_raster("habitats.tif") as habitat_map,
yg.read_raster('elevation.tif') as elevation_map,
yg.read_shape('species123.geojson') as range_polygon
):
refined_habitat = habitat_map.isin([...species habitat codes...])
refined_elevation = (elevation_map >= species_min) & (elevation_map <= species_max)
aoh = refined_habitat * refined_elevation * range_polygon
print(f'area for species 123: {aoh.sum()}')
aoh.to_geotiff("result.tif")
Operators
Add, subtract, multiple, divide
Pixel-wise addition, subtraction, multiplication or division (both true and floor division), and remainder. Either between arrays, or with constants:
with (
yg.read_raster('test1.tif') as layer1,
yg.read_raster('test2.tif') as layer2
):
result = layer1 + layer2
result.to_geotiff("result.tif")
or
Boolean testing
Testing for equality, less than, less than or equal, greater than, and greater than or equal are supported on layers, along with logical or and logical and, as per this example, where elevation_upper and elevation_lower are scalar values:
filtered_elevation = (min_elevation_map <= elevation_upper) & (max_elevation_map >= elevation_lower)
Power
Pixel-wise raising to a constant power:
Log, Exp, Clip, etc.
The following math operators common to numpy and other libraries are currently supported:
- abs
- all
- any
- ceil
- clip
- exp
- exp2
- floor
- isin
- log
- log2
- log10
- maximum
- minimum
- nan_to_num
- round
- sum
Typically these can be invoked either on a layer as a method:
Or via the operators module, as it's sometimes nicer to do it this way when chaining together operations in a single expression:
2D matrix convolution
To facilitate image processing algorithms you can supply a weight matrix to generate a processed image. Currently this support only works for square weight matrices of an odd size.
For example, to apply a blur function to a raster:
blur_filter = np.array([
[0.0, 0.1, 0.0],
[0.1, 0.6, 0.1],
[0.0, 0.1, 0.0],
])
with yg.read_raster('original.tif') as layer1:
calc = layer1.conv2d(blur_filter)
calc.to_geotiff("blurred.tif")
Type conversion
Similar to numpy and other Python numerical libraries, Yirgacheffe will automatically deal with simple type conversion where possible, however sometimes explicit conversion is either necessary or desired. Similar to numpy, there is an as_type operator that lets you set the conversion:
from yirgacheffe.operations import DataType
with yg.read_raster('float_data.tif') as float_layer:
int_layer = float_layer.as_type(DataType.Int32)
Area conversion
In general when you are working on multiple data sources, Yirgacheffe will attempt to work out a minimal valid geospatial area of coverage for a given expression, based on the notion that for any region outside that defined in a GeoTIFF Yirgacheffe will fill in the gaps with zero. For example, if you multiply two areas together that have different geospatial areas, then it will take the intersection as it knows that pixels outside those defined in either layer will be zero, so the multiplication will result in zero, and so the intersection avoids unnecessary computation. However, for say addition, it will have to take the union of the areas as 0 plus a value is still a value.
However, sometimes you will want to override this behaviour, and specify a layer is cropped to a certain subregion of the overall calculation is done to match the area of some other data you have. For that you can use the as_area operator.
with yg.read_raster('sweden.tif') as swe:
stockholm = swe.as_area(yg.Area(...coordinates of stockholm...))
Projection conversion
Another problem frequently encountered is having two data sources in different map projections. In that case Yirgacheffe can be used to reproject the data. You can either do this at load time, using read_raster_like, or within an expression:
with yg.read_raster('sweden_epsg:6004.tif') as swe:
swe_wgs84 = swe.as_projection("eps:4326", yg.ResampleingMethod.Nearest)
swe_wgs84.to_geotiff('sweden_wgs84.tif')
You must specify a resampling method to be used.
Note: be careful using this as it is an expensive operation. If you use the same data in multiple calculations and reproject it each time you are probably better to reproject it and save it to disk, and then use that version multiple times.
Storing the results of expressions
There are three ways to store the result of a computation.
Saving to a GeoTIFF
In all the above examples we use the to_geotiff call, to which you pass a filename for a GeoTIFF, into which the results will be written. You can optionally pass a callback to save which will be called for each chunk of data processed and give you the amount of progress made so far as a number between 0.0 and 1.0:
def print_progress(p)
print(f"We have made {p * 100} percent progress")
...
calc.to_geotiff(result, callback=print_progress)
Aggregations
The alternative is to call aggregation functions such as sum, min, or max which will give you a single value by aggregating the data within the layer or expression:
with (
RasterLayer.layer_from_file(...) as area_layer,
VectorLayer(...) as mask_layer
):
intersection = yg.find_intersection([area_layer, mask_layer])
area_layer.set_intersection_window(intersection)
mask_layer.set_intersection_window(intersection)
calc = area_layer * mask_layer
total_area = calc.sum()
As numpy arrays
Finally, if you want to read the pixel values of either a layer or an expression, you can call read_array: