import rastereasy
Rastereasy
In rastereasy, one of the main class to deal with georeferenced images is Geoimage. It contains the main functions (for resampling, reprojection, stack, plotting, …). Let’s start by opening and plotting an image, some bands or making color compositions
1) first step: open and visualize an image
name_im='./data/demo/sentinel.tif'
image=rastereasy.Geoimage(name_im)
Print basic information of the image
image.info()
- Size of the image:
- Rows (height): 1000
- Cols (width): 1000
- Bands: 12
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.04099599, 38.39058840)
- Driver: GTiff
- Data type: int16
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, '9': 9, '10': 10, '11': 11, '12': 12}
As we can see, spectral bands car named '1', '2', ... in a dictionnary where the name is associated with its number in the image (from 1 to the number of bands).
Alternatively , we can give more explicit names to the spectral bands by specifying them :
# Alternatively, we can give names to bands
names = {"CO" : 1,"B": 2,"G":3,"R":4,"RE1":5,"RE2":6,"RE3":7,"NIR":8,"WA":9,"SWIR1":10,"SWIR2":11,"SWIR3":12}
image_names=rastereasy.Geoimage(name_im,names=names)
image_names.info()
# Note : this is equavalent to image.change_names(names)
- Size of the image:
- Rows (height): 1000
- Cols (width): 1000
- Bands: 12
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.04099599, 38.39058840)
- Driver: GTiff
- Data type: int16
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'CO': 1, 'B': 2, 'G': 3, 'R': 4, 'RE1': 5, 'RE2': 6, 'RE3': 7, 'NIR': 8, 'WA': 9, 'SWIR1': 10, 'SWIR2': 11, 'SWIR3': 12}
Create a color composite using bands 4, 3, and 2
image.colorcomp(['4', '3', '2'])
<Figure size 640x480 with 0 Axes>
2) Example of easy manipulations (resample, reproject. See dedicated examples for more details / options
Resample the image to a resolution of 2 meters
image_resampled = image.resampling(5)
image_resampled.info()
image.resampling(15,inplace=True) # inplace=True : modify directly the image
image.info()
- Size of the image:
- Rows (height): 2000
- Cols (width): 2000
- Bands: 12
- Spatial resolution: 5.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.04101857, 38.39056574)
- Driver: GTiff
- Data type: int16
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, '9': 9, '10': 10, '11': 11, '12': 12}
- Size of the image:
- Rows (height): 666
- Cols (width): 666
- Bands: 12
- Spatial resolution: 15.015015015015013 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.04097334, 38.39061114)
- Driver: GTiff
- Data type: int16
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, '9': 9, '10': 10, '11': 11, '12': 12}
Reproject the image to the EPSG:4326 coordinate system
image_reproject = image.reproject("EPSG:27700")
image_reproject.info()
image_reproject.colorcomp()
- Size of the image:
- Rows (height): 733
- Cols (width): 733
- Bands: 12
- Spatial resolution: 19.64992127401193 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.04099917, 38.39064373)
- Driver: GTiff
- Data type: int16
- Projection system: EPSG:27700
- Nodata: -32768.0
- Given names for spectral bands:
{'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, '9': 9, '10': 10, '11': 11, '12': 12}
<Figure size 640x480 with 0 Axes>
Save the image
image_reproject.save('modify_image.tif')
One can combine operations
names = {"CO" : 1,"B": 2,"G":3,"R":4,"RE1":5,"RE2":6,"RE3":7,"NIR":8,"WA":9,"SWIR1":10,"SWIR2":11,"SWIR3":12}
image_names=rastereasy.Geoimage(name_im,names=names)
image_names.info()
image_modified=image_names.crop(50,100,120,300).resampling(30).reproject("EPSG:27700")
image_modified.colorcomp()
image_modified.info()
- Size of the image:
- Rows (height): 1000
- Cols (width): 1000
- Bands: 12
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.04099599, 38.39058840)
- Driver: GTiff
- Data type: int16
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'CO': 1, 'B': 2, 'G': 3, 'R': 4, 'RE1': 5, 'RE2': 6, 'RE3': 7, 'NIR': 8, 'WA': 9, 'SWIR1': 10, 'SWIR2': 11, 'SWIR3': 12}
<Figure size 640x480 with 0 Axes>
- Size of the image:
- Rows (height): 23
- Cols (width): 61
- Bands: 12
- Spatial resolution: 39.353156172510154 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.07928894, 38.36432299)
- Driver: GTiff
- Data type: int16
- Projection system: EPSG:27700
- Nodata: -32768.0
- Given names for spectral bands:
{'CO': 1, 'B': 2, 'G': 3, 'R': 4, 'RE1': 5, 'RE2': 6, 'RE3': 7, 'NIR': 8, 'WA': 9, 'SWIR1': 10, 'SWIR2': 11, 'SWIR3': 12}
Some general functions
Sum, product, comparison, …
Weighted sum of images
# 1. Read and plot two images
name_im1='./data/demo/im1.tif'
name_im2='./data/demo/im2.tif'
image1=rastereasy.Geoimage(name_im1)
image2=rastereasy.Geoimage(name_im2)
image1.colorcomp([4,3,2])
image2.colorcomp([4,3,2])
<Figure size 640x480 with 0 Axes>
<Figure size 640x480 with 0 Axes>
# Compute a weighted sum
weighted1=3
weighted2=5
image_sum=(weighted1*image1+weighted2*image2)/(weighted1+weighted2)
image_sum.info()
image_sum.colorcomp([4,3,2])
- Size of the image:
- Rows (height): 300
- Cols (width): 300
- Bands: 12
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.07261482, 38.35885970)
- Driver: GTiff
- Data type: float64
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, '9': 9, '10': 10, '11': 11, '12': 12}
<Figure size 640x480 with 0 Axes>
Multiplication, division
# Multiplication of images
image_product = image1*image2
image_product.info()
image_product.colorcomp([4,3,2])
- Size of the image:
- Rows (height): 300
- Cols (width): 300
- Bands: 12
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.07261482, 38.35885970)
- Driver: GTiff
- Data type: float64
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, '9': 9, '10': 10, '11': 11, '12': 12}
<Figure size 640x480 with 0 Axes>
# Division of images
image_divide = image1/image2
image_divide.info()
image_divide.colorcomp([4,3,2])
- Size of the image:
- Rows (height): 300
- Cols (width): 300
- Bands: 12
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.07261482, 38.35885970)
- Driver: GTiff
- Data type: float64
- Projection system: EPSG:32637
- Nodata: -32768.0
- Given names for spectral bands:
{'1': 1, '2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, '9': 9, '10': 10, '11': 11, '12': 12}
<Figure size 640x480 with 0 Axes>
Compute ndvi
name_im='./data/demo/im2.tif'
image=rastereasy.Geoimage(name_im)
# Get nir band (number 8)
nir = image.select_bands(8)
# Get red band (number 4)
red = image.select_bands(4)
# Compute ndvi
ndvi=(nir-red)/(nir+red)
# check info
ndvi.info()
# plot
ndvi.visu(colorbar=True)
# Note : ndvi.colorcomp() works also
ndvi.hist(title='histogram of NDVI')
- Size of the image:
- Rows (height): 300
- Cols (width): 300
- Bands: 1
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.03647565, 38.39059829)
- Driver: GTiff
- Data type: float64
- Projection system: EPSG:32637
- Nodata: nan
- Given names for spectral bands:
{'8': 1}
<Figure size 640x480 with 0 Axes>
<Figure size 640x480 with 0 Axes>
Identifying vegetation with (use >, <)
threshold=0.3
veget=ndvi>=threshold
veget.info()
veget.visu()
- Size of the image:
- Rows (height): 300
- Cols (width): 300
- Bands: 1
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.03647565, 38.39059829)
- Driver: GTiff
- Data type: bool
- Projection system: EPSG:32637
- Nodata: 0
- Given names for spectral bands:
{'8': 1}
<Figure size 640x480 with 0 Axes>
Identifying vegetation with ndvi.where() (similar than with numpy)
ndvi.info()
veget=ndvi.where(ndvi>=threshold,1,0)
veget.info()
veget.visu()
- Size of the image:
- Rows (height): 300
- Cols (width): 300
- Bands: 1
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.03647565, 38.39059829)
- Driver: GTiff
- Data type: float64
- Projection system: EPSG:32637
- Nodata: nan
- Given names for spectral bands:
{'8': 1}
- Size of the image:
- Rows (height): 300
- Cols (width): 300
- Bands: 1
- Spatial resolution: 10.0 meters / degree (depending on projection system)
- Central point latitude - longitude coordinates: (7.03647565, 38.39059829)
- Driver: GTiff
- Data type: float64
- Projection system: EPSG:32637
- Nodata: nan
- Given names for spectral bands:
{'8': 1}
<Figure size 640x480 with 0 Axes>
Identifying vegetation with indexing
image2=image.copy()
image2[veget==0]=0
image2.colorcomp([4,3,2])
<Figure size 640x480 with 0 Axes>
Min, max of the image, per band, row, …
help(image.min)
print('Overall minimum of the image : ',image.min(),'\n')
help(image.max)
print('Maximum of the image for each band: ', image.max(axis = 'pixel'),'\n')
help(image.std)
print('Standard deviation of the image along the lines: ',image.std(axis = 'row'),'\n')
help(image.mean)
print('Overall mean of the image along the lines: ',image.mean(axis='row'),'\n')
print('Mean of the image along the bands: ',image.mean(axis='band'),'\n')
Help on method min in module rastereasy.rastereasy:
min(axis=None) method of rastereasy.rastereasy.Geoimage instance
Calculate the minimum value along a specified axis.
Parameters
----------
axis : {'band', 'row', 'col', None}, optional
The axis along which to compute the minimum:
- 'band': Minimum across spectral bands for each pixel
- 'row': Minimum across rows (lines) for each band and column
- 'col': Minimum across columns for each band and row
- None: Global minimum of the entire image
Default is None.
Returns
-------
float or numpy.ndarray
- If axis=None: A single value representing the global minimum
- If axis='band': Array with shape (nb_rows,nb_cols) containing mins along bands
- If axis='row': Array with shape (nb_bands,nb_cols) containing mins along rows
- If axis='col': Array with shape (nb_bands,nb_rows) containing mins along cols
- If axis='pixel': Array with shape (nb_bands) containing mins along all pixels for each band
Raises
------
ValueError
If an invalid axis is specified
Examples
--------
>>> min_value = image.min() # Global minimum value
>>> print(f"Minimum pixel value: {min_value}")
>>>
>>> band_mins = image.min(axis='pixel') # Minimum along all pixels for each band
Overall minimum of the image : 92
Help on method max in module rastereasy.rastereasy:
max(axis=None) method of rastereasy.rastereasy.Geoimage instance
Calculate the maximum value along a specified axis.
Parameters
----------
axis : {'band', 'row', 'col', None}, optional
The axis along which to compute the maximum:
- 'band': Maximum across spectral bands for each pixel
- 'row': Maximum across rows (lines) for each band and column
- 'col': Maximum across columns for each band and row
- None: Global maximum of the entire image
Default is None.
Returns
-------
float or numpy.ndarray
- If axis=None: A single value representing the global maximum
- If axis='band': Array with shape (nb_rows,nb_cols) containing max along bands
- If axis='row': Array with shape (nb_bands,nb_cols) containing max along rows
- If axis='col': Array with shape (nb_bands,nb_rows) containing max along cols
- If axis='pixel': Array with shape (nb_bands) containing maxs along all pixels for each band
Raises
------
ValueError
If an invalid axis is specified
Examples
--------
>>> max_value = image.max() # Global maximum value
>>> print(f"Maximum pixel value: {max_value}")
>>>
>>> band_maxs = image.max(axis='pixel') # Maximum along all pixels for each band
Maximum of the image for each band: [1438 2290 2478 2680 2958 4702 5118 5352 5281 4455 4268 3839]
Help on method std in module rastereasy.rastereasy:
std(axis=None) method of rastereasy.rastereasy.Geoimage instance
Calculate the standard deviation along a specified axis.
Parameters
----------
axis : {'band', 'row', 'col', None}, optional
The axis along which to compute the standard deviation:
- 'band': Std dev across spectral bands for each pixel
- 'row': Std dev across rows (lines) for each band and column
- 'col': Std dev across columns for each band and row
- None: Global standard deviation of the entire image
Default is None.
Returns
-------
float or numpy.ndarray
- If axis=None: A single value representing the global std
- If axis='band': Array with shape (nb_rows,nb_cols) containing std along bands
- If axis='row': Array with shape (nb_bands,nb_cols) containing std along rows
- If axis='col': Array with shape (nb_bands,nb_rows) containing std along cols
- If axis='pixel': Array with shape (nb_bands) containing std along all pixels for each band
Raises
------
ValueError
If an invalid axis is specified
Examples
--------
>>> std_value = image.std() # Global standard deviation
>>> print(f"Standard deviation of pixel values: {std_value}")
>>>
>>> band_stds = image.std(axis='pixel') # Standard deviation along all pixels for each band
Notes
-----
The standard deviation is a measure of the amount of variation or dispersion in the values.
A low standard deviation indicates that the values tend to be close to the mean, while a high
standard deviation indicates that the values are spread out over a wider range.
Standard deviation of the image along the lines: [[150.9863209 150.9863209 150.9863209 ... 11.8936491 11.8936491
11.8936491 ]
[234.94052004 239.04220185 240.4099581 ... 18.04782874 16.25493942
14.76004027]
[253.91344448 259.85286116 264.11683391 ... 16.11090176 15.35531034
15.88821051]
...
[131.55089302 131.55089302 131.55089302 ... 10.02921233 10.02921233
10.02921233]
[247.08075189 276.0303183 276.0303183 ... 8.34862863 8.34862863
7.42651705]
[396.38312851 423.73239274 423.73239274 ... 9.05840984 9.05840984
8.38442736]]
Help on method mean in module rastereasy.rastereasy:
mean(axis=None) method of rastereasy.rastereasy.Geoimage instance
Calculate the mean value along a specified axis.
Parameters
----------
axis : {'band', 'row', 'col', None}, optional
The axis along which to compute the mean:
- 'band': Mean across spectral bands for each pixel
- 'row': Mean across rows (lines) for each band and column
- 'col': Mean across columns for each band and row
- None: Global mean of the entire image
Default is None.
Returns
-------
float or numpy.ndarray
- If axis=None: A single value representing the global mean
- If axis='band': Array with shape (nb_rows,nb_cols) containing mean along bands
- If axis='row': Array with shape (nb_bands,nb_cols) containing mean along rows
- If axis='col': Array with shape (nb_bands,nb_rows) containing mean along cols
- If axis='pixel': Array with shape (nb_bands) containing mean along all pixels for each band
Raises
------
ValueError
If an invalid axis is specified
Examples
--------
>>> mean_value = image.mean() # Global mean value
>>> print(f"Mean pixel value: {mean_value}")
>>>
>>> band_means = image.mean(axis='pixel') # Mean along all pixels for each band
Notes
-----
This method uses np.nanmean, which ignores NaN values in the calculation.
If you have NaN values as nodata, they won't affect the mean calculation.
Overall mean of the image along the lines: [[ 844.97 844.97 844.97 ... 233.23333333
233.23333333 233.23333333]
[1026.59333333 1025.18 1029.07333333 ... 259.75666667
261.71666667 262.65666667]
[1345.12666667 1345.64333333 1346.77333333 ... 461.78666667
460.83333333 461.09 ]
...
[2702.07666667 2702.07666667 2702.07666667 ... 203.57
203.57 203.57 ]
[3376.90666667 3367.49333333 3367.49333333 ... 183.02
183.02 182.01333333]
[3127.89 3097.09666667 3097.09666667 ... 153.95666667
153.95666667 153.27333333]]
Mean of the image along the bands: [[2341.16666667 2325.75 2332.58333333 ... 249.16666667
249.41666667 251. ]
[2315.08333333 2321.16666667 2321.16666667 ... 251.16666667
251. 250.66666667]
[2315.25 2318.5 2319.08333333 ... 251.08333333
249.83333333 249.75 ]
...
[1803.16666667 1747.66666667 1760.66666667 ... 235.58333333
237.83333333 237.16666667]
[1801. 1750.5 1754.41666667 ... 233.25
236.83333333 235.75 ]
[1813.25 1823.91666667 1844.58333333 ... 229.58333333
233.33333333 236. ]]
# Extract the 30th line
image[30,:]
array([[ 933, 933, 933, ..., 221, 221, 221],
[1278, 1250, 1254, ..., 284, 271, 252],
[1678, 1624, 1664, ..., 456, 453, 448],
...,
[2828, 2828, 2828, ..., 209, 209, 209],
[3710, 3614, 3614, ..., 193, 193, 191],
[3570, 3470, 3470, ..., 169, 169, 161]],
shape=(12, 300), dtype=int16)
# all spectral values in a given pixel
row=20
col=10
print('spectral values in',row,',',col,' : ',image[row,col])
print('size of all spectral values in row',row,':', image[row,:].shape)
print('size of all spectral values in col',col,':', image[:,col].shape)
# check consistency
pixel_row = 15
pixel_col = 22
print('diff with table with values at row ',row,' in pixel col',pixel_col,':', image[row,:][:,pixel_col]-image[row,pixel_col])
print('diff with table with values at col ',col,' in pixel row',pixel_row,':', image[:,col][:,pixel_row]-image[pixel_row,col])
spectral values in 20 , 10 : [ 986 1254 1582 2088 2317 2445 2552 2646 2729 2727 3543 3365]
size of all spectral values in row 20 : (12, 300)
size of all spectral values in col 10 : (12, 300)
diff with table with values at row 20 in pixel col 22 : [0 0 0 0 0 0 0 0 0 0 0 0]
diff with table with values at col 10 in pixel row 15 : [0 0 0 0 0 0 0 0 0 0 0 0]
im1=rastereasy.Geoimage('./data/demo/sentinel.tif')
im2=rastereasy.Geoimage('./data/demo/sentinel.tif')+4