modules¶
generic¶
input-output of mapdata¶
- dhdt.generic.mapping_io.make_geo_im(Z, R, crs, fName, meta_descr='project Eratosthenes', no_dat=nan, sun_angles='az:360-zn:90', date_created='-0276-00-00')¶
Create georeferenced tiff file (a GeoTIFF)
- Parameters
Z (numpy.array, size=(m,n)) – band image
R (list, size=(1,6)) – GDAL georeference transform of an image
crs (string) – coordinate reference string
fname (string) – filename for the image with extension
no_dat (datatype, integer) – no data value
sun_angles (string) – string giving meta data about the illumination angles
date_created (string) – string given the acquistion date in +YYYY-MM-DD
Examples
>>> import os >>> fpath = os.path.join(os.getcwd(), "data.jp2")
>>> (I, spatialRef, geoTransform, targetPrj) = read_geo_image(fpath)
>>> Z_ones = np.zeros(Z.shape, dtype=bool) >>> make_geo_im(Z_ones, geoTransformM, spatialRefM, 'ones.tif')
- dhdt.generic.mapping_io.make_multispectral_vrt(df, fpath=None, fname='multispec.vrt')¶
virtual raster tile (VRT) is a description of datasets written in an XML format, it eases the display of multi-spectral data or other means.
- Parameters
df (pandas.DataFrame) – organization of the different spectral bands
fpath (string) – path of the directory of interest
fname (string) – file name of the virtual raster tile
Examples
>>> s2_df = list_central_wavelength_msi() >>> s2_df = s2_df[s2_df['gsd']==10] >>> s2_df, datastrip_id = get_s2_image_locations(IM_PATH, s2_df)
create file in directory where imagery is situated >>> make_multispectral_vrt(s2_df, fname=’multi_spec.vrt’)
- dhdt.generic.mapping_io.read_geo_image(fname, boi=array([], dtype=float64), no_dat=None)¶
This function takes as input the geotiff name and the path of the folder that the images are stored, reads the image and returns the data as an array
- Parameters
fname (string) – geotiff file name and path.
boi (numpy.array, size=(k,1)) – bands of interest, if a multispectral image is read, a selection can be specified
no_dat ({integer,float}) – no data value
- Returns
data ({numpy.array, numpy.masked.array}, size=(m,n), ndim=2) – data array of the band
spatialRef (string) – osr.SpatialReference in well known text
geoTransform (tuple, size=(8,)) – affine transformation coefficients.
targetprj (osgeo.osr.SpatialReference() object) – coordinate reference system (CRS)
See also
make_geo_imbasic function to write out geographic data
read_geo_infobasic function to get meta data of geographic imagery
read_nc_imbasic function to read geographic raster out of nc-file
Examples
>>> import os >>> from dhdt.generic.mapping_io import read_geo_image >>> fpath = os.path.join(os.getcwd(), "data.jp2" )
>>> I, spatialRef, geoTransform, targetPrj = read_geo_image(fpath)
>>> I_ones = np.zeros(I.shape, dtype=bool) >>> make_geo_im(I_ones, geoTransformM, spatialRefM, "ones.tif")
- dhdt.generic.mapping_io.read_geo_info(fname)¶
This function takes as input the geotiff name and the path of the folder that the images are stored, reads the geographic information of the image
- Parameters
fname (string) – path and file name of a geotiff image
- Returns
spatialRef (string) – osr.SpatialReference in well known text
geoTransform (tuple, size=(8,1)) – affine transformation coefficients, but also giving the image dimensions
targetprj (osgeo.osr.SpatialReference() object) – coordinate reference system (CRS)
rows (integer, {x ∈ ℕ | x ≥ 0}) – number of rows in the image, that is its height
cols (integer, {x ∈ ℕ | x ≥ 0}) – number of collumns in the image, that is its width
bands (integer, {x ∈ ℕ | x ≥ 1}) – number of bands in the image, that is its depth
See also
read_geo_imagebasic function to import geographic imagery data
handle Sentinel-2 imagery¶
- dhdt.generic.handler_sentinel2.get_crs_from_mgrs_tile(tile_code)¶
- Parameters
tile_code (string) – US Military Grid Reference System (MGRS) tile code
- Returns
crs – target projection system
- Return type
osgeo.osr.SpatialReference
See also
get_utmzone_from_mgrs_tile,get_utmzone_from_mgrs_tileNotes
- The tile structure is a follows “AABCC”
“AA” utm zone number, starting from the East, with steps of 8 degrees
“B” latitude zone, starting from the South, with steps of 6 degrees
- dhdt.generic.handler_sentinel2.get_epsg_from_mgrs_tile(tile_code)¶
- Parameters
tile_code (string, e.g.: '05VMG') – MGRS tile coding
- Returns
epsg – code used to denote a certain database entry
- Return type
integer
See also
dhdt.generic.get_utmzone_from_mgrs_tile,dhdt.generic.get_crs_from_mgrs_tileNotes
- The tile structure is a follows “AABCC”
“AA” utm zone number, starting from the East, with steps of 8 degrees
“B” latitude zone, starting from the South, with steps of 6 degrees
The following acronyms are used:
CRS : coordinate reference system
EPSG : european petroleum survey group (a coordinate refenence database)
MGRS : US military grid reference system
UTM : universal transverse mercator
WGS : world geodetic system
- dhdt.generic.handler_sentinel2.get_generic_s2_raster(tile_code, spac=10, tile_path=None)¶
Create spatial metadata of a Sentinel-2, so no downloading is needed.
- Parameters
tile_code (string) – mgrs tile coding, which is also given in the filename (“TXXXXX”)
spac (integer, {10,20,60}, unit=m) – pixel spacing of the raster
tile_path (string) – path to the MGRS tiling file
- Returns
tuple – georeference transform of an image, including shape
pyproj.crs.crs.CRS – coordinate reference system (CRS)
Notes
- The tile structure is a follows “AABCC”
“AA” utm zone number, starting from the East, with steps of 8 degrees
“B” latitude zone, starting from the South, with steps of 6 degrees
The following acronyms are used:
CRS : coordinate reference system
MGRS : US military grid reference system
s2 : Sentinel-2
- dhdt.generic.handler_sentinel2.get_s2_dict(s2_df)¶
given a dataframe with a filename within, create a dictionary with scene and satellite specific metadata.
- Parameters
s2_df (pd.dataframe) – metadata and general multi spectral information about the MSI instrument that is onboard Sentinel-2
- Returns
s2_dict – dictionary with scene specific meta data, giving the following information: * ‘COSPAR’ : string
id of the satellite
- ’NORAD’string
id of the satellite
- ’instruments’{‘MSI’}
specific instrument
- ’launch_date’: string, e.g.:’+2015-06-23’
date when the satellite was launched into orbit
- ’orbit’: string
specifies the type of orbit
- ’full_path’: string
the location where the root folder is situated
- ’MTD_DS_path’: string
the location where metadata about the datastrip is present
- ’MTD_TL_path’: string
the location where metadata about the tile is situated
- ’QI_DATA_path’: string
the location where metadata about detector readings are present
- ’IMG_DATA_path’: string
the location where the imagery is present
- ’date’: string, e.g.: ‘+2019-10-25’
date of acquisition
- ’tile_code’: string, e.g.: ‘05VMG’
MGRS tile coding
- ’relative_orbit’: integer, {x ∈ ℕ}
orbit from which the imagery were taken
- Return type
dictionary
Notes
The metadata is scattered over a folder structure for Sentinel-2, it has typically the following set-up:
* S2X_MSIL1C_20XX... <- path is given in s2_dict["full_path"] ├ AUX_DATA ├ DATASTRIP │ └ DS_XXX_XXXX... <- path is given in s2_dict["MTD_DS_path"] │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... <- path is given in s2_dict["MTD_TL_path"] │ ├ AUX_DATA │ ├ IMG_DATA │ │ ├ TXXX_X..XXX.jp2 │ │ └ TXXX_X..XXX.jp2 │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
The following acronyms are used:
AUX : auxiliary
COSPAR : committee on space research international designator
DS : datastrip
IMG : imagery
L1C : product specification,i.e.: level 1, processing step C
MGRS : military grid reference system
MTD : metadata
MSI : multi spectral instrument
NORAD : north american aerospace defense satellite catalog number
TL : tile
QI : quality information
s2 : Sentinel-2
See also
dhdt.generic.get_s2_image_locations
- dhdt.generic.handler_sentinel2.get_s2_image_locations(fname, s2_df)¶
The Sentinel-2 imagery are placed within a folder structure, where one folder has an ever changing name. Fortunately this function finds the path from the metadata
- Parameters
fname (string) – path string to the Sentinel-2 folder
s2_df (pandas.dataframe) – index of the bands of interest
- Returns
s2_df_new (pandas.dataframe) – dataframe series with relative folder and file locations of the bands
datastrip_id (string) – folder name of the metadata
Examples
>>> import os >>> fpath = '/Users/Data/' >>> sname = 'S2A_MSIL1C_20200923T163311_N0209_R140_T15MXV_20200923T200821.SAFE' >>> fname = 'MTD_MSIL1C.xml' >>> full_path = os.path.join(fpath, sname, fname) >>> im_paths,datastrip_id = get_s2_image_locations(full_path) >>> im_paths ['GRANULE/L1C_T15MXV_A027450_20200923T163313/IMG_DATA/T15MXV_20200923T163311_B01', 'GRANULE/L1C_T15MXV_A027450_20200923T163313/IMG_DATA/T15MXV_20200923T163311_B02'] >>> datastrip_id 'S2A_OPER_MSI_L1C_DS_VGS1_20200923T200821_S20200923T163313_N02.09'
- dhdt.generic.handler_sentinel2.get_utmzone_from_tile_code(tile_code)¶
- Parameters
tile_code (string, e.g.: '05VMG') – MGRS tile coding
- Returns
utmzone – code used to denote the number of the projection column of UTM
- Return type
integer
Notes
- The tile structure is a follows “AABCC”
“AA” utm zone number, starting from the East, with steps of 8 degrees
“B” latitude zone, starting from the South, with steps of 6 degrees
The following acronyms are used:
CRS : coordinate reference system
MGRS : US military grid reference system
UTM : universal transverse mercator
WGS : world geodetic system
- dhdt.generic.handler_sentinel2.meta_s2string(s2_str)¶
get meta information of the Sentinel-2 file name
- Parameters
s2_str (string) – filename of the L1C data
- Returns
s2_time (string) – date “+YYYY-MM-DD”
s2_orbit (string) – relative orbit “RXXX”
s2_tile (string) – tile code “TXXXXX”
Examples
>>> s2_str = 'S2A_MSIL1C_20200923T163311_N0209_R140_T15MXV_20200923T200821.SAFE' >>> s2_time, s2_orbit, s2_tile = meta_s2string(s2_str) >>> s2_time '+2020-09-23' >>> s2_orbit 'R140' >>> s2_tile 'T15MXV'
handle Landsat imagery¶
- dhdt.generic.handler_landsat.get_wrs_url(version=2)¶
get the location where the geometric data of Landsats path-row footprints are situated, for info see [wwwUSGS]
References
- dhdt.generic.handler_landsat.meta_LSstring(LSstr)¶
get meta data of the Landsat file name input: LSstr string filename of the L1TP data output: LStime string date “+YYYY-MM-DD”
LSpath string satellite path “PXXX” LSrow string tile row “RXXX”
handle RapidEye imagery¶
- dhdt.generic.handler_rapideye.meta_REstring(REstr)¶
get meta information of the RapidEye file name
- Parameters
REstr (string) – filename of the rapideye data
- Returns
REtime (string) – date “+YYYY-MM-DD”
REtile (string) – tile code “TXXXXX”
REsat (string) – which RapedEye satellite
Examples
>>> REstr = '568117_2012-09-10_RE2_3A_Analytic.tif' >>> REtime, REtile, REsat = meta_REstring(REstr) >>> REtime '+2012-09-10' >>> REtile '568117' >>> REsat 'RE2'
input¶
read Sentinel-2 (meta)data¶
Sentinel-2 is a high-resolution (10-60 m) multi-spectral satellite constellation from European Space Agency and the European Union. The following functions make is easier to read such data.
- dhdt.input.read_sentinel2.dn2toa_s2(dn)¶
convert the digital numbers of Sentinel-2 to top of atmosphere (TOA)
Notes
https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-2-msi/level-1c/algorithm
- dhdt.input.read_sentinel2.get_flight_bearing_from_detector_mask_s2(Det)¶
- Parameters
Det (numpy.array, size=(m,n), ndim={2,3}, dtype=integer) – array with numbers of the different detectors in Sentinel-2
- Returns
ψ – general bearing of the satellite
- Return type
float, unit=degrees
- dhdt.input.read_sentinel2.get_flight_bearing_from_gnss_s2(path, spatialRef, rec_tim, fname='MTD_DS.xml')¶
get the direction/argument/heading of the Sentinel-2 acquisition
- Parameters
path (string) – directory where the meta-data is located
spatialRef (osgeo.osr.SpatialReference) – projection system used
rec_tim ({integer,float}) – time stamp of interest
fname (string) – filename of the meta-data file
- Returns
az – array with argument values, based upon the map projection given
- Return type
numpy.array, size=(m,1), float
See also
get_s2_image_locationsto get the datastrip id
Notes
The metadata is scattered over the file structure of Sentinel-2, L1C
* S2X_MSIL1C_20XX... ├ AUX_DATA ├ DATASTRIP │ └ DS_XXX_XXXX... │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... │ ├ AUX_DATA │ ├ IMG_DATA │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
The following acronyms are used:
DS : datastrip
TL : tile
QI : quality information
AUX : auxiliary
MTD : metadata
MSI : multi spectral instrument
L1C : product specification,i.e.: level 1, processing step C
- dhdt.input.read_sentinel2.get_flight_orientation_s2(ds_path, fname='MTD_DS.xml', s2_dict=None)¶
get the flight path and orientations of the Sentinel-2 satellite during acquisition.
It is also possible to only give the dictionary, as this has such metadata within.
- Parameters
ds_path (string) – directory where the meta-data is located
fname (string, default='MTD_DS.xml') – filename of the meta-data file
s2_dict (dictonary, default=None) – metadata of the Sentinel-2 platform
- Returns
sat_time (numpy.array, size=(m,1), unit=nanosec) – satellite timestamps
sat_angles (numpy.array, size=(m,3)) – satellite orientation angles, given in ECEF
s2_dict (dictonary) – updated with keys: “time”, “quat”
See also
get_flight_path_s2get positions of the flight path
- dhdt.input.read_sentinel2.get_flight_path_s2(ds_path, fname='MTD_DS.xml', s2_dict=None)¶
It is also possible to only give the dictionary, as this has such metadata within.
- Parameters
ds_path (string) – location of the metadata file
fname (string, default='MTD_DS.xml') – name of the xml-file that has the metadata
s2_dict (dictonary) – metadata of the Sentinel-2 platform
- Returns
sat_time (numpy.array, size=(m,1), dtype=np.datetime64, unit=ns) – time stamp of the satellite positions
sat_xyz (numpy.array, size=(m,3), dtype=float, unit=meter) – 3D coordinates of the satellite within an Earth centered Earth fixed (ECEF) frame.
sat_err (numpy.array, size=(m,3), dtype=float, unit=meter) – error estimate of the 3D coordinates given by “sat_xyz”
sat_uvw (numpy.array, size=(m,3), dtype=float, unit=meter sec-1) – 3D velocity vectors of the satellite within an Earth centered Earth fixed (ECEF) frame.
s2_dict (dictonary) – updated with keys: “gps_xyz”, “gps_uvw”, “gps_tim”, “gps_err”, “altitude”, “speed”
See also
get_flight_orientation_s2get the quaternions of the flight path
Examples
Following the file and metadata structure of scihub:
>>> import os >>> import numpy as np
>>> S2_dir = '/data-dump/examples/' >>> S2_name = 'S2A_MSIL1C_20200923T163311_N0209_R140_T15MXV_20200923T200821.SAFE' >>> fname = os.path.join(S2_dir, S2_name, 'MTD_MSIL1C.xml') >>> s2_df = list_central_wavelength_s2()
>>> s2_df, datastrip_id = get_s2_image_locations(fname, s2_df) >>> path_det = os.path.join(S2_dir, S2_name, 'DATASTRIP', datastrip_id[17:-7])
>>> sat_tim, sat_xyz, sat_err, sat_uvw = get_flight_path_s2(path_det)
Notes
The metadata structure of MTD_DS looks like:
* MTD_DS.xml └ n1:Level-1C_DataStrip_ID ├ n1:General_Info ├ n1:Image_Data_Info ├ n1:Satellite_Ancillary_Data_Info │ ├ Time_Correlation_Data_List │ ├ Ephemeris │ │ ├ GPS_Number_List │ │ ├ GPS_Points_List │ │ │ └ GPS_Point │ │ │ ├ POSITION_VALUES │ │ │ ├ POSITION_ERRORS │ │ │ ├ VELOCITY_VALUES │ │ │ ├ VELOCITY_ERRORS │ │ │ ├ GPS_TIME │ │ │ ├ NSM │ │ │ ├ QUALITY_INDEX │ │ │ ├ GDOP │ │ │ ├ PDOP │ │ │ ├ TDOP │ │ │ ├ NOF_SV │ │ │ └ TIME_ERROR │ │ └ AOCS_Ephemeris_List │ ├ Attitudes │ ├ Thermal_Data │ └ ANC_DATA_REF │ ├ n1:Quality_Indicators_Info └ n1:Auxiliary_Data_Info
The following acronyms are used:
AOCS : attitude and orbit control system
CAMS : Copernicus atmosphere monitoring service
DEM : digital elevation model
GDOP : geometric dilution of precision
GPS : global positioning system
GRI : global reference image
IERS : international earth rotation and reference systems service
IMT : instrument measurement time
NSM : navigation solution method
NOF SV : number of space vehicles
PDOP : position dilution of precision
TDOP : time dilution of precision
- dhdt.input.read_sentinel2.get_integration_and_sampling_time_s2(ds_path, fname='MTD_DS.xml', s2_dict=None)¶
- Parameters
ds_path (string) – location where metadata of datastrip is situated
fname (string, default='MTD_DS.xml') – metadata filename
s2_dict (dictionary, default=None) – metadata dictionary, can be used instead of the string based method
Notes
The metadata is scattered over the file structure of Sentinel-2, L1C
* S2X_MSIL1C_20XX... ├ AUX_DATA │ └ DS_XXX_XXXX... │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... │ ├ AUX_DATA │ ├ IMG_DATA │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
The following acronyms are used:
DS : datastrip
TL : tile
QI : quality information
AUX : auxiliary
MTD : metadata
MSI : multi spectral instrument
L1C : product specification,i.e.: level 1, processing step C
- dhdt.input.read_sentinel2.get_intrinsic_camera_mat_s2(s2_df, det, boi)¶
- Parameters
s2_df (pandas.DataFrame) – metadata and general multispectral information about the MSI instrument that is onboard Sentinel-2
det (list) – list with detector numbers that are within the tile
boi (list) – list with bands of interest
- Returns
F – fundamental matrix
- Return type
numpy.array, size=(3,3)
Notes
The pushbroom detector onboard Sentinel-2 has the following configuration:
nadir det1 | ====# #===# #===# #===# #===# #===# #===# #===# #===# #===# #===# #==== | det12 ===== : 2592 pixels for 10 meter, 1296 pixels for 20-60 meter # : 98 pixels overlap for 20 meter dataReferences
- HG94
Hartley & Gupta, “Linear pushbroom cameras” Proceedings of the European conference on computer vision, 1994.
- Mo10
Moons et al. “3D reconstruction from multiple images part 1: Principles.” Foundations and trends® in computer graphics and vision, vol.4(4) pp.287-404, 2010.
- dhdt.input.read_sentinel2.list_central_wavelength_msi()¶
create dataframe with metadata about Sentinel-2
- Returns
df – metadata and general multispectral information about the MSI instrument that is onboard Sentinel-2, having the following collumns:
center_wavelength, unit=µm : central wavelength of the band
full_width_half_max, unit=µm : extent of the spectral sensativity
bandid : number for identification in the meta data
resolution, unit=m : spatial resolution of a pixel
along_pixel_size, unit=µm : physical size of the sensor
across_pixel_size, unit=µm : size of the photosensative sensor
focal_length, unit=m : lens characteristic of the telescope
field_of_view, unit=degrees : angle of swath of instrument
crossdetector_parallax, unit=degress : in along-track direction
name : general name of the band, if applicable
solar_illumination, unit=W m-2 µm-1 :
mostly following the naming convension of the STAC EO extension [stac]
- Return type
pandas.dataframe
See also
dhdt.input.read_landsat.list_central_wavelength_olifor the instrument data of Landsat 8 or 9
dhdt.input.read_aster.list_central_wavelength_asfor the instrument data of ASTER on the Terra satellite
dhdt.input.read_venus.list_central_wavelength_vsscfor the instrument data of VSSC on the VENµS satellite
Notes
The Multi Spectral instrument (MSI) has a Tri Mirror Anastigmat (TMA) telescope, which is opened at F/4 and has a focal length of 60 centimeter.
The crossdetector parallex is given here as well, see [La15]. While the crossband parallax is 0.018 for VNIR and 0.010 for the SWIR detector, respectively. The dimensions of the Silicon CMOS detector are given in [MG10]. While the SWIR detector is based on a MCT sensor.
The following acronyms are used:
MSI : multi-spectral instrument
MCT : mercury cadmium telluride, HgCdTe
TMA : trim mirror anastigmat
VNIR : very near infra-red
SWIR : short wave infra-red
CMOS : silicon complementary metal oxide semiconductor, Si
Examples
make a selection by name:
>>> boi = ['red', 'green', 'blue', 'near infrared'] >>> s2_df = list_central_wavelength_msi() >>> s2_df = s2_df[s2_df['common_name'].isin(boi)] >>> s2_df wavelength bandwidth resolution name bandid B02 492 66 10 blue 1 B03 560 36 10 green 2 B04 665 31 10 red 3 B08 833 106 10 near infrared 7
similarly you can also select by pixel resolution:
>>> s2_df = list_central_wavelength_msi() >>> tw_df = s2_df[s2_df['gsd']==20] >>> tw_df.index Index(['B05', 'B06', 'B07', 'B8A', 'B11', 'B12'], dtype='object')
References
- La15
Languille et al. “Sentinel-2 geometric image quality commissioning: first results” Proceedings of the SPIE, 2015.
- MG10
Martin-Gonthier et al. “CMOS detectors for space applications: From R&D to operational program with large volume foundry”, Proceedings of the SPIE conference on sensors, systems, and next generation satellites XIV, 2010.
- stac
- dhdt.input.read_sentinel2.read_band_s2(path, band=None)¶
This function takes as input the Sentinel-2 band name and the path of the folder that the images are stored, reads the image and returns the data as an array
- Parameters
path (string) – path of the folder, or full path with filename as well
band (string, optional) – Sentinel-2 band name, for example ‘04’, ‘8A’.
- Returns
data (numpy.array, size=(_,_)) – array of the band image
spatialRef (string) – projection
geoTransform (tuple) – affine transformation coefficients
targetprj (osr.SpatialReference object) – spatial reference
See also
list_central_wavelength_msicreates a dataframe for the MSI instrument
read_stack_s2reading several Sentinel-2 bands at once into a stack
dhdt.generic.read_geo_imagebasic function to import geographic imagery
dataExamples
>>> path = '/GRANULE/L1C_T15MXV_A027450_20200923T163313/IMG_DATA/' >>> band = '02' >>> _,spatialRef,geoTransform,targetprj = read_band_s2(path, band)
>>> spatialRef 'PROJCS["WGS 84 / UTM zone 15S",GEOGCS["WGS 84",DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93], PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000], PARAMETER["false_northing",10000000],UNIT["metre",1, AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH], AUTHORITY["EPSG","32715"]]' >>> geoTransform (600000.0, 10.0, 0.0, 10000000.0, 0.0, -10.0) >>> targetprj <osgeo.osr.SpatialReference; proxy of <Swig Object of type 'OSRSpatialReferenceShadow *' at 0x7f9a63ffe450> >
- dhdt.input.read_sentinel2.read_cloud_mask(path_meta, geoTransform)¶
- Parameters
path_meta (string) – directory where meta data is situated, ‘MSK_CLOUDS_B00.gml’ is typically the file of interest
geoTransform (tuple) – affine transformation coefficients
- Returns
msk_clouds – array which highlights where clouds could be
- Return type
numpy.array, size=(m,n), dtype=int
- dhdt.input.read_sentinel2.read_detector_mask(path_meta, boi, geoTransform)¶
create array of with detector identification
Sentinel-2 records in a pushbroom fasion, collecting reflectance with a ground resolution of more than 270 kilometers. This data is stacked in the flight direction, and then cut into granules. However, the sensorstrip inside the Multi Spectral Imager (MSI) is composed of several CCD arrays.
This function collects the geometry of these sensor arrays from the meta- data. Since this is stored in a gml-file.
- Parameters
path_meta (string) – path where the meta-data is situated.
boi (pandas.DataFrame) – list with bands of interest
geoTransform (tuple, size={(6,),(8,)}) – affine transformation coefficients
- Returns
det_stack – array where each pixel has the ID of the detector, of a specific band
- Return type
numpy.ma.array, size=(msk_dim[0:1],len(boi)), dtype=int8
See also
list_central_wavelength_msicreates a dataframe for the MSI instrument
read_view_angles_s2read grid of Sentinel-2 observation angles
Notes
The metadata is scattered over the file structure of Sentinel-2, L1C
* S2X_MSIL1C_20XX... ├ AUX_DATA ├ DATASTRIP │ └ DS_XXX_XXXX... │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... │ ├ AUX_DATA │ ├ IMG_DATA │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
The following acronyms are used:
DS : datastrip
TL : tile
QI : quality information
AUX : auxiliary
MTD : metadata
MSI : multi spectral instrument
L1C : product specification,i.e.: level 1, processing step C
Examples
>>> path_meta = '/GRANULE/L1C_T15MXV_A027450_20200923T163313/QI_DATA' >>> boi = ['red', 'green', 'blue', 'near infrared'] >>> s2_df = list_central_wavelength_msi() >>> boi_df = s2_df[s2_df['common_name'].isin(boi)] >>> geoTransform = (600000.0, 10.0, 0.0, 10000000.0, 0.0, -10.0) >>> >>> det_stack = read_detector_mask(path_meta, boi_df, geoTransform)
- dhdt.input.read_sentinel2.read_detector_time_s2(path, fname='MTD_DS.xml', s2_df=None)¶
get the detector metadata of the relative detector timing
- Parameters
path (string) – path where the meta-data is situated
fname (string) – file name of the metadata.
- Returns
det_time (numpy.array, numpy.datetime64) – detector time
det_name (list, size=(13,)) – list of the different detector codes
det_meta (numpy.array, size=(13,4)) –
- metadata of the individual detectors, each detector has the following:
spatial resolution [m]
minimal spectral range [nm]
maximum spectral range [nm]
mean spectral range [nm]
line_period (float, numpy.timedelta64) – temporal sampling distance
Notes
The sensor blocks are arranged as follows, with ~98 pixels overlap:
┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ |DET02| |DET04| |DET06| |SCA08| |SCA10| |SCA12| └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ |DET01| |DET03| |SCA05| |SCA07| |SCA09| |SCA11| └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ └-----┘a forward looking array, while band order is reverse for aft looking <-2592 10m pixels-> ┌-----------------┐ #*# satellite | B02 | | ├-----------------┤ | flight | B08 | | direction ├-----------------┤ | | B03 | v ├-----------------┤ etc. the detector order is B02, B08, B03, B10, B04, B05, B11, B06, B07, B8A, B12, B01 and B09
- dhdt.input.read_sentinel2.read_geotransform_s2(path, fname='MTD_TL.xml', resolution=10)¶
- Parameters
path (string) – location where the meta data is situated
fname (string) – file name of the meta-data file
resolution ({float,integer}, unit=meters, default=10) – resolution of the grid
- Returns
geoTransform – affine transformation coefficients
- Return type
tuple, size=(1,6)
Notes
The metadata is scattered over the file structure of Sentinel-2, L1C
* S2X_MSIL1C_20XX... ├ AUX_DATA ├ DATASTRIP │ └ DS_XXX_XXXX... │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... │ ├ AUX_DATA │ ├ IMG_DATA │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
- The metadata structure is as follows:
MTD_TL.xml
- └ n1:Level-1C_Tile_ID
├ n1:General_Info ├ n1:Geometric_Info │ ├ Tile_Geocoding │ │ ├ HORIZONTAL_CS_NAME │ │ ├ HORIZONTAL_CS_CODE │ │ ├ Size : resolution={“10”,”20”,”60”} │ │ │ ├ NROWS │ │ │ └ NCOLS │ │ └ Geoposition │ │ ├ ULX │ │ ├ ULY │ │ ├ XDIM │ │ └ YDIM │ └ Tile_Angles └ n1:Quality_Indicators_Info
The following acronyms are used:
DS : datastrip
TL : tile
QI : quality information
AUX : auxiliary
MTD : metadata
MSI : multi spectral instrument
L1C : product specification,i.e.: level 1, processing step C
- dhdt.input.read_sentinel2.read_mean_sun_angles_s2(path, fname='MTD_TL.xml')¶
Read the xml-file of the Sentinel-2 scene and extract the mean sun angles.
- Parameters
path (string) – path where xml-file of Sentinel-2 is situated
- Returns
Zn (float) – Mean solar zentih angle of the scene, in degrees.
Az (float) – Mean solar azimuth angle of the scene, in degrees
Notes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--|--> + | └----> East & xThe angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / |-- zenith angle | / | / | /| |/ |/ | elevation angle └---- └------
The metadata is scattered over the file structure of Sentinel-2, L1C
* S2X_MSIL1C_20XX... ├ AUX_DATA ├ DATASTRIP │ └ DS_XXX_XXXX... │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... │ ├ AUX_DATA │ ├ IMG_DATA │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
The following acronyms are used:
s2 : Sentinel-2
DS : datastrip
TL : tile
QI : quality information
AUX : auxiliary
MTD : metadata
MSI : multi spectral instrument
L1C : product specification,i.e.: level 1, processing step C
- dhdt.input.read_sentinel2.read_orbit_number_s2(path, fname='MTD_MSIL1C.xml')¶
read the orbit number of the image
- dhdt.input.read_sentinel2.read_sensing_time_s2(path, fname='MTD_TL.xml')¶
- Parameters
path (string) – path where the meta-data is situated
fname (string) – file name of the metadata.
- Returns
rec_time – time of image acquisition
- Return type
numpy.datetime64
See also
get_s2_granual_idNotes
The recording time is the average sensing within the tile, which is a combination of forward and backward looking datastrips. Schematically, this might look like:
x recording time of datastrip, i.e.: the start of the strip × recording time of tile, i.e.: the weighted average _________ datastrip ____x____/ /_________ / / / / / / / / / / / / / / / / / / / / / / / / / ┌---/------┐ / / / | / |/ / / | / × / / / |/ /| / / /--------/-┘ / / / tile / / / / / / / / / / / / / / / / / / / _________/ / _________ /________/Examples
demonstrate the code when in a Scihub data structure
>>> import os >>> fpath = '/Users/Data/' >>> sname = 'S2A_MSIL1C_20200923T163311_N0209_R140_T15MXV_20200923T200821.SAFE' >>> fname = 'MTD_MSIL1C.xml' >>> s2_path = os.path.join(fpath, sname, fname) >>> s2_df,_ = get_s2_image_locations(fname, s2_df) >>> s2_df['filepath'] B02 'GRANULE/L1C_T15MXV_A027450_20200923T163313/IMG_DATA/T115MXV...' B03 'GRANULE/L1C_T15MXV_A027450_20200923T163313/IMG_DATA/T15MXV...' >>> full_path = '/'.join(s2_df['filepath'][0].split('/')[:-2]) 'GRANULE/L1C_T15MXV_A027450_20200923T163313' >>> rec_time = read_sensing_time_s2(path, fname='MTD_TL.xml') >>> rec_time
- dhdt.input.read_sentinel2.read_stack_s2(s2_df)¶
read imagery data of interest into an three dimensional np.array
- Parameters
s2_df (pandas.Dataframe) – metadata and general multispectral information about the MSI instrument that is onboard Sentinel-2
- Returns
im_stack (numpy.array, size=(_,_,_)) – array of the band image
spatialRef (string) – projection
geoTransform (tuple, size=(8,1)) – affine transformation coefficients
targetprj (osr.SpatialReference object) – spatial reference
See also
list_central_wavelength_msicreates a dataframe for the MSI instrument
get_s2_image_locationsprovides dataframe with specific file locations
read_band_s2reading a single Sentinel-2 band
Examples
>>> s2_df = list_central_wavelength_msi() >>> s2_df = s2_df[s2_df['gsd'] == 10] >>> s2_df,_ = get_s2_image_locations(IM_PATH, s2_df)
>>> im_stack, spatialRef, geoTransform, targetprj = read_stack_s2(s2_df)
- dhdt.input.read_sentinel2.read_sun_angles_s2(path, fname='MTD_TL.xml')¶
This function reads the xml-file of the Sentinel-2 scene and extracts an array with sun angles, as these vary along the scene.
- Parameters
path (string) – path where xml-file of Sentinel-2 is situated
- Returns
Zn (numpy.array, size=(m,n), dtype=float) – array of the solar zenith angles, in degrees.
Az (numpy.array, size=(m,n), dtype=float) – array of the solar azimuth angles, in degrees.
Notes
The computation of the solar angles is done in two steps: 1) Computation of the solar angles in J2000; 2) Transformation of the vector to the mapping frame.
The outputs of the first step is the solar direction normalized vector with the Earth-Sun distance, considering that the direction of the sun is the same at the centre of the Earth and at the centre of the Sentinel-2 satellite.
Attitude of the satellite platform are used to rotate the solar vector to the mapping frame. Also Ground Image Calibration Parameters (GICP Diffuser Model) are used to transform from the satellite to the diffuser, as Sentinel-2 has a forward and backward looking sensor configuration. The diffuser model also has stray-light correction and a Bi-Directional Reflection Function model.
The angle(s) are declared in the following coordinate frame:
^ North & y | - <--|--> + | +----> East & xThe angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / |-- zenith angle | / | / | /| |/ |/ | elevation angle +---- +------
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | i | x --------+--------> --------+--------> | | | | image | j map | based v based |The metadata is scattered over the file structure of Sentinel-2, L1C
* S2X_MSIL1C_20XX... ├ AUX_DATA ├ DATASTRIP │ └ DS_XXX_XXXX... │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... │ ├ AUX_DATA │ ├ IMG_DATA │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
The metadata structure looks like:
* MTD_TL.xml └ n1:Level-1C_Tile_ID ├ n1:General_Info ├ n1:Geometric_Info │ ├ Tile_Geocoding │ └ Tile_Angles │ ├ Sun_Angles_Grid │ │ ├ Zenith │ │ │ ├ COL_STEP │ │ │ ├ ROW_STEP │ │ │ └ Values_List │ │ └ Azimuth │ │ ├ COL_STEP │ │ ├ ROW_STEP │ │ └ Values_List │ ├ Mean_Sun_Angle │ └ Viewing_Incidence_Angles_Grids └ n1:Quality_Indicators_Info
The following acronyms are used:
s2 : Sentinel-2
DS : datastrip
TL : tile
QI : quality information
AUX : auxiliary
MTD : metadata
MSI : multi spectral instrument
L1C : product specification,i.e.: level 1, processing step C
- dhdt.input.read_sentinel2.read_view_angles_s2(path, fname='MTD_TL.xml', det_stack=array([], dtype=float64), boi_df=None)¶
This function reads the xml-file of the Sentinel-2 scene and extracts an array with viewing angles of the MSI instrument.
- Parameters
path (string) – path where xml-file of Sentinel-2 is situated
fname (string) – the name of the metadata file, sometimes this is changed
boi_df (pandas.dataframe, default=4) – each band has a somewhat minute but different view angle
- Returns
Zn (numpy.ma.array, size=(m,n), dtype=float) – masked array of the solar zenith angles, in degrees.
Az (numpy.ma.array, size=(m,n), dtype=float) – masked array of the solar azimuth angles, in degrees.
See also
list_central_wavelength_s2Notes
The azimuth angle is declared in the following coordinate frame:
^ North & y | - <--|--> + | └----> East & xThe angles related to the satellite are as follows:
#*# #*# satellite ^ / ^ /| | / | / | nadir |-- zenith angle | / v | / | /| |/ |/ | elevation angle └----- surface └------
The metadata is scattered over the file structure of Sentinel-2, L1C
* S2X_MSIL1C_20XX... ├ AUX_DATA ├ DATASTRIP │ └ DS_XXX_XXXX... │ └ QI_DATA │ └ MTD_DS.xml <- metadata about the data-strip ├ GRANULE │ └ L1C_TXXXX_XXXX... │ ├ AUX_DATA │ ├ IMG_DATA │ ├ QI_DATA │ └ MTD_TL.xml <- metadata about the tile ├ HTML ├ rep_info ├ manifest.safe ├ INSPIRE.xml └ MTD_MSIL1C.xml <- metadata about the product
The metadata structure looks like:
* MTD_TL.xml └ n1:Level-1C_Tile_ID ├ n1:General_Info ├ n1:Geometric_Info │ ├ Tile_Geocoding │ └ Tile_Angles │ ├ Sun_Angles_Grid │ ├ Mean_Sun_Angle │ └ Viewing_Incidence_Angles_Grids │ ├ Zenith │ │ ├ COL_STEP │ │ ├ ROW_STEP │ │ └ Values_List │ └ Azimuth │ ├ COL_STEP │ ├ ROW_STEP │ └ Values_List └ n1:Quality_Indicators_Info
The following acronyms are used:
s2 : Sentinel-2
DS : datastrip
TL : tile
QI : quality information
AUX : auxiliary
MTD : metadata
MSI : multi spectral instrument
L1C : product specification,i.e.: level 1, processing step C
read Landsat (meta)data¶
Landsat 8 & 9 are the latest fleet of satellites that are part of a legacy. Both satellites are build by NASA and data provision is done through USGS. The following function makes reading such data easier.
- dhdt.input.read_landsat.get_sca_numbering_ls(ang, boi='BAND04')¶
- Parameters
ang (dictonary) – metadata of the Landsat scene
boi (string) – band of intertest, typically the data of band 4 is given by VAA.tif
- Returns
Sca – array with detector numbers
- Return type
np.array, size=(m,n), dtype=int, range=0…14
Notes
The sensor blocks are arranged as follows, with ~28 pixels overlap:
┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ |SCA02| |SCA04| |SCA06| |SCA08| |SCA10| |SCA12| |SCA14| └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ ┌-----┐ |SCA01| |SCA03| |SCA05| |SCA07| |SCA09| |SCA11| |SCA13| └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ └-----┘ └-----┘See also
- dhdt.input.read_landsat.list_central_wavelength_oli()¶
create dataframe with metadata about OLI on Landsat8 or Landsat9
- Returns
df – metadata and general multispectral information about the OLI instrument that is onboard Landsat8, having the following collumns:
wavelength, unit=µm : central wavelength of the band
bandwidth, unit=µm : extent of the spectral sensativity
bandid : number for identification in the meta data
resolution, unit=m : spatial resolution of a pixel
field_of_view, unit=degrees : angle of swath of instrument
name : general name of the band, if applicable
- Return type
pandas.DataFrame
See also
dhdt.input.read_sentinel2.list_central_wavelength_msifor the instrument data of MSI on Sentinel-2
dhdt.input.read_aster.list_central_wavelength_asfor the instrument data of ASTER on the Terra satellite
dhdt.input.read_venus.list_central_wavelength_vsscfor the instrument data of VSSC on the VENµS satellite
Notes
The time lag between panchromatic and red is 0.52 seconds, while the time difference between panchromatic and the green band is 0.65 sec., see [dM16]. The detector samping time is 4.2 msec, while the total staring time is 1.2 seconds. The detector arrays (or senso chip assemblies: SCA) of landsat are configured as follows:
a forward looking SCA, while bandorder is reverse for aft looking <--494 pixels--> ┌--------------┐ #*# satellite | green | | ├--------------┤ | flight | red | | direction ├--------------┤ | |near infrared | v ├--------------┤ | not in use | ├--------------┤ | blue | ├--------------┤ | panchromatic | └--------------┘ <--988 pixels-->
References
- dM16
de Michele, et al. “Volcanic plume elevation model and its velocity derived from Landsat 8” Remote sensing of environment, vol.176 pp.219-224, 2016.
Examples
make a selection by name:
>>> boi = ['red', 'green', 'blue', 'near infrared'] >>> ls_df = list_central_wavelength_oli() >>> ls_df = ls_df[ls_df['name'].isin(boi)] >>> ls_df wavelength bandwidth resolution common_name bandid B02 482 60 30 blue 2 B03 561 57 30 green 3 B04 655 37 30 red 4 B05 865 28 30 near infrared 5
similarly you can also select by pixel resolution:
>>> ls_df = list_central_wavelength_oli() >>> pan_df = ls_df[ls_df['resolution']==15] >>> pan_df.index Index(['B08'], dtype='object')
- dhdt.input.read_landsat.read_band_ls(path, band='B8')¶
read landsat image from path
- pathstring
path of the folder, or full path with filename as well
- bandstring, optional
Landsat8 band name, for example ‘4’, ‘B8’.
- datanp.array, size=(_,_)
array of the band image
- spatialRefstring
projection
- geoTransformtuple
affine transformation coefficients
- targetprjosr.SpatialReference object
spatial reference
list_central_wavelength_ls
>>> path = '/LC08_L1TP_018060_20200904_20200918_02_T1/' >>> band = '2' >>> _,spatialRef,geoTransform,targetprj = read_band_ls(path, band)
>>> spatialRef
- ‘PROJCS[“WGS 84 / UTM zone 15N”,GEOGCS[“WGS 84”,DATUM[“WGS_1984”, …
>>> geoTransform (626385.0, 30.0, 0.0, 115815.0, 0.0, -30.0) >>> targetprj <osgeo.osr.SpatialReference; proxy of <Swig Object of type 'OSRSpatial ...
- dhdt.input.read_landsat.read_metadata_ls(dir_path, suffix='MTL.txt')¶
Convert Landsat MTL file to dictionary of metadata values
read Terra ASTER data¶
The ASTER instrument onboard of the Terra satellite has been collecting data since 1999, but will be decommisioned soon. Nonetheless, its archive is very valuable and the following functions make reading such data easier.
- dhdt.input.read_aster.get_flight_path_as(as_path, fname='AST_L1A_*.Ancillary_Data.txt', as_dict=None)¶
- Parameters
as_path (string) –
fname (string) – file name of the accillary data
as_dict (dictionary) –
- dhdt.input.read_aster.list_central_wavelength_as()¶
create dataframe with metadata about Terra’s ASTER
- Returns
df – metadata and general multispectral information about the MSI instrument that is onboard Sentinel-2, having the following collumns:
wavelength : central wavelength of the band
bandwidth : extent of the spectral sensativity
bandid : number for identification in the meta data
resolution : spatial resolution of a pixel
name : general name of the band, if applicable
irradiance :
- Return type
pandas.DataFrame
See also
dhdt.input.read_sentinel2.list_central_wavelength_msifor the instrument data of MSI on Sentinel-2
dhdt.input.read_landsat.list_central_wavelength_olifor the instrument data of Landsat 8 or 9
dhdt.input.read_venus.list_central_wavelength_vsscfor the instrument data of VSSC on the VENµS satellite
Examples
make a selection by name:
>>> boi = ['red', 'green', 'blue'] >>> as_df = list_central_wavelength_as() >>> as_df = as_df[as_df['name'].isin(boi)] >>> as_df wavelength bandwidth resolution name bandid irradiance B01 560 80 15 blue 0 1848.0 B02 660 60 15 green 1 1549.0 B3B 810 100 15 red 2 1114.0 B3N 810 100 15 red 3 1114.0
similarly you can also select by pixel resolution:
>>> as_df = list_central_wavelength_as() >>> sw_df = as_df[as_df['gsd']==30] >>> sw_df.index Index(['B04', 'B05', 'B06', 'B07', 'B08', 'B09'], dtype='object')
- dhdt.input.read_aster.read_band_as_hdf(path, band='3N')¶
This function takes as input the ASTER L1T hdf-filename and the path of the folder that the images are stored, reads the image and returns the data as an array, plus associated geographical/projection metadata.
- Parameters
path (string) – path of the folder, or full path with filename as well
band (string, optional) – Terra ASTER band index, for example ‘02’, ‘3N’.
- Returns
data (numpy.ndarray, size=(_,_)) – array of the band image
spatialRef (string) – projection
geoTransform (tuple) – affine transformation coefficients
targetprj (osr.SpatialReference object) – spatial reference
References
- dhdt.input.read_aster.read_stack_as_l1(path, as_df)¶
read the band stack of Level-1 data
- Parameters
path (string) – location where the imagery data is situated
as_df (pandas.DataFrame) – metadata and general multispectral information about the instrument that is onboard Terra
- Returns
im_stack (numpy.ndarray, size=(_,_,_)) – array of the band image
spatialRef (string) – projection
geoTransform (tuple, size=(8,1)) – affine transformation coefficients
targetprj (osr.SpatialReference object) – spatial reference
See also
read RapidEye data¶
The RapidEye constellation composed of five commercial mini-satellites, image collection started in 2008 and halted in 2020. The following functions make reading such data easier.
- dhdt.input.read_rapideye.list_central_wavelength_re()¶
create dataframe with metadata about RapidEye
- Returns
df – metadata and general multispectral information about the MSI instrument that is onboard Sentinel-2, having the following collumns:
center_wavelength, unit=µm : central wavelength of the band
full_width_half_max, unit=µm : extent of the spectral sensativity
bandid : number for identification in the meta data
resolution, unit=m : spatial resolution of a pixel
common_name : general name of the band, if applicable
irradiance, unit=W m-2 μm-1 : exo-atmospheric radiance
relative_timing, unit=ms : relative sampling time difference
- Return type
pandas.dataframe
Notes
The detector arrays of the satellite are configured as follows[Kr13]_:
78 mm <--------------> ┌--------------┐ #*# satellite | red | | ├--------------┤ | flight | red edge | | direction ├--------------┤ | |near infrared | v ├--------------┤ | | | | | | | | ├--------------┤ ^ | not in use | | 6.5 μm ├--------------┤ v | green | ├--------------┤ | blue | └--------------┘
References
- Kr13
Krauß, et al. “Traffic flow estimation from single satellite images” Archives of the ISPRS, vol.XL-1/WL3, 2013.
Examples
make a selection by name:
>>> boi = ['red', 'green', 'blue'] >>> re_df = list_central_wavelength_re() >>> re_df = re_df[re_df['common_name'].isin(boi)] >>> re_df wavelength bandwidth resolution name bandid irradiance B01 475 70 5 blue 0 1997.8 B02 555 70 5 green 1 1863.5 B03 657 55 5 red 2 1560.4
similarly you can also select by bandwidth:
>>> re_df = list_central_wavelength_re() >>> sb_df = re_df[re_df['bandwidth']<=60] >>> sb_df.index Index(['B03', 'B04'], dtype='object')
- dhdt.input.read_rapideye.read_band_re(band, path)¶
read specific band of RapidEye image
This function takes as input the RapidEye band number and the path of the folder that the images are stored, reads the image and returns the data as an array
- Parameters
band (string) – RapidEye band number.
path (string) – path to folder with imagery.
- Returns
data (np.array, size=(m,n), dtype=integer) – data of one RadidEye band.
spatialRef (string) – osr.SpatialReference in well known text
geoTransform (tuple, size=(6,1)) – affine transformation coefficients.
targetprj (osgeo.osr.SpatialReference() object) – coordinate reference system (CRS)
- dhdt.input.read_rapideye.read_stack_re(path)¶
read specific band of RapidEye image
This function takes as input the RapidEye band number and the path of the folder that the images are stored, reads the image and returns the data as an array
- Parameters
path (string) – path to folder with imagery.
- Returns
data (np.array, size=(m,n,_), dtype=integer) – data of one RadidEye band.
spatialRef (string) – osr.SpatialReference in well known text
geoTransform (tuple, size=(6,1)) – affine transformation coefficients.
targetprj (osgeo.osr.SpatialReference() object) – coordinate reference system (CRS)
- dhdt.input.read_rapideye.read_sun_angles_re(path)¶
read the sun angles, corresponding to the RapidEye acquisition
- Parameters
path (string) – path where xml-file of RapidEye is situated.
- Returns
Zn (np.array, size=(m,n), dtype=float) – array of the Zenith angles
Az (np.array, size=(m,n), dtype=float) – array of the Azimtuh angles
Notes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--┼--> + | +----> East & xThe angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / ├-- zenith angle | / | / | /| |/ |/ | elevation angle +---- surface +--┴---
read PlanetScope data¶
The PlanetScope constellation is a commercial fleet of micro-satellites. The following function makes reading such data easier.
- dhdt.input.read_planetscope.read_band_ps(dir_path, fname, no_dat=None)¶
This function takes as input the PlanetScope file name and the path of the folder that the images are stored, reads the image and returns the data as an array
- Parameters
dir_path (string) – path of the folder
fname (string) – PlanetScope image name
- Returns
data (np.array, size=(_,_)) – array of the band image
spatialRef (string) – projection
geoTransform (tuple) – affine transformation coefficients
targetprj (osr.SpatialReference object) – spatial reference
See also
list_central_wavelength_ps
- dhdt.input.read_planetscope.read_detector_type_ps(path, fname='*_metadata.xml')¶
Notes
Dove Pilot has the following sensor array:
<---- 4032 pixels ----> ┌------------------------┐ ^ | | | | | | 2672 pixels | RGB | | | (Bayer pattern) | | | | | └------------------------┘ v
PlanetScope (PS Instrument) has the following sensor array:
<---- 6600 pixels ----> ┌------------------------┐ ^ | | | | RGB | | 4400 pixels | (Bayer pattern) | | ├------------------------┤ | | | | | NIR | | | | | └------------------------┘ v
Dove-R (PS2.SD Instrument) has the following sensor array:
<---- 6600 pixels ----> ┌------------------------┐ ^ | Green II (B04) | | ├------------------------┤ | 4400 pixels | Red (B05) | | ├------------------------┤ | | NIR (B13) | | ├------------------------┤ | ^ | Blue (B02) | | | 1100 pixels └------------------------┘ v v
Superdove(PSB.SD Instrument) has the following bands: Coastal Blue, Blue, Green, Green I, Green II, Yellow, Red, Rededge, NIR
<---- 8800 pixels ----> ┌------------------------+ ^ | Blue (B02) | | ├------------------------┤ | 5280 pixels | Red (B05) | | ├------------------------┤ | | Green (B03) | | ├------------------------┤ | | GreenII(B04) | | ├------------------------┤ | | Yellow (B06) | | ├------------------------┤ | | Rededge(B10) | | ├------------------------┤ | | NIR (B13) | | ├------------------------┤ | ^ | Coastal Blue | | | 660 pixels └------------------------┘ v v
read VENµS data¶
The VENµS satellite is a demonstration satellite, that acquires over specific pre-defined regions. It has a high repeat rate in the order of one or two days, if cloud cover permits. The following functions make reading of such data easier.
- dhdt.input.read_venus.list_central_wavelength_vssc()¶
create dataframe with metadata about VENµS
- Returns
df – metadata and general multispectral information about the SuperSpectral Camera that is onboard VENµS, having the following collumns:
center_wavelength, unit=µm : central wavelength of the band
full_width_half_max, unit=µm : extent of the spectral sensativity
bandid : number for identification in the meta data
gsd, unit=m : spatial resolution of a pixel
common_name : general name of the band, if applicable
- Return type
pandas.dataframe
See also
dhdt.input.read_sentinel2.list_central_wavelength_msifor the instrument data of MSI on Sentinel-2
dhdt.input.read_landsat.list_central_wavelength_olifor the instrument data of Landsat 8 or 9
dhdt.input.read_aster.list_central_wavelength_asfor the instrument data of ASTER on the Terra satellite
Notes
The sensor is composed of four detectors, which each have three arrays. These multi-spectral bands (BXX) are arranged as follows:
┌-----B06------┐ #*# satellite ├-----B04------┤ | flight └-----B03------┘ detector id: IV | direction ┌-----B09------┐ v ├-----B08------┤ └-----B07------┘ detector id: III ┌-----B11------┐ ├-----B12------┤ └-----B10------┘ detector id: II ┌-----B01------┐ ├-----B02------┤ └-----B05------┘ detector id: I
References
- Di22
Dick, et al. “VENμS: mission characteristics, final evaluation of the first phase and data production” Remote sensing vol.14(14) pp.3281, 2022.
- dhdt.input.read_venus.read_stack_vn(vn_df)¶
read imagery data of interest into an three dimensional np.array
- Parameters
vn_df (pandas.Dataframe) – metadata and general multispectral information about the VSSC instrument that is onboard VENµS
- Returns
im_stack (numpy.array, size=(_,_,_)) – array of the band image
spatialRef (string) – projection
geoTransform (tuple, size=(8,1)) – affine transformation coefficients
targetprj (osr.SpatialReference object) – spatial reference
preprocessing¶
image transforms¶
functions adjusting the gray scale of an image
- dhdt.preprocessing.image_transforms.cum_hist(img)¶
compute the cumulative histogram of an image
- Parameters
img (numpy.array, size=(m,n), ndim=2) – grid with intensity values
- Returns
cdf – emperical cumulative distribution function, i.e. a cumulative histogram.
- Return type
numpy.array, size=(2**n,), dtype=integer
- dhdt.preprocessing.image_transforms.gamma_adjustment(Z, gamma=1.0)¶
transform intensity in non-linear way
enhances high intensities when gamma>1
- Parameters
Z (numpy.ndarray, size=(m,n), dtype=integer) – array with intensity values
gamma (float) – power law parameter
- Returns
Z_new – array with transform
- Return type
numpy.ndarray, size=(m,n)
- dhdt.preprocessing.image_transforms.general_midway_equalization(Z)¶
equalization of multiple imagery, based on [De04].
- Parameters
Z (numpy.array, dim=3, type={uint8,uint16}) – stack of imagery
- Returns
Z_new – harmonized stack of imagery
- Return type
numpy.array, dim=3, type={uint8,uint16}
References
- De04
Delon, “Midway image equalization”, Journal of mathematical imaging and Vision”, vol.21(2) pp.119-134, 2004
Notes
The following acronyms are used:
CDF - cumulative distribution function
- dhdt.preprocessing.image_transforms.get_histogram(img)¶
calculate the normalized histogram of an image
- Parameters
img (numpy.ndarray, size=(m,n)) – gray scale image array
- Returns
hist – histogram array
- Return type
numpy.ndarray, size={(2**8,2), (2**16,2)}
- dhdt.preprocessing.image_transforms.gompertz_adjustment(Z, a=None, b=None)¶
transform intensity in non-linear way, so high reflected surfaces like snow have more emphasized.
- Parameters
Z (numpy.ndarray, size=(m,n), dtype=integer) – array with intensity values
a (float) – parameters for the Gompertz function
b (float) – parameters for the Gompertz function
- Returns
Z_new – grid with transform
- Return type
numpy.ndarray, size=(m,n)
- dhdt.preprocessing.image_transforms.high_pass_im(Im, radius=10)¶
- Parameters
Im (numpy.array, ndim={2,3}, size=(m,n,_)) – data array with intensities
radius (float, unit=pixels) – hard threshold of the sphere in Fourier space
- Returns
Inew – data array with high-frequency signal
- Return type
numpy.array, ndim={2,3}, size=(m,n,_)
- dhdt.preprocessing.image_transforms.histogram_equalization(img, img_ref)¶
transform the intensity values of an array, so they have a similar histogram as the reference.
- Parameters
img (numpy.array, size=(m,n), ndim=2) – asd
img_ref (numpy.array, size=(m,n), ndim=2) –
- Returns
img – transformed image
- Return type
numpy.array, size=(m,n), ndim=2
- dhdt.preprocessing.image_transforms.hyperbolic_adjustment(Z, intercept)¶
transform intensity through an hyperbolic sine function
- Parameters
Z (numpy.ndarray, size=(m,n), dtype=float) – grid with intensity values
- Returns
Z_new – grid with transform
- Return type
numpy.ndarray, size=(m,n), dtype=float, range=-1…+1
- dhdt.preprocessing.image_transforms.inverse_tangent_transformation(x)¶
enhances high and low intensities
- Parameters
x (numpy.ndarray, size=(m,n)) – grid with intensity values
- Returns
fx – grid with transform
- Return type
numpy.ndarray, size=(_,_)
Examples
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> x = np.arange(-1,1,.001) >>> y = inverse_tangent_transformation(x) >>> plt.figure(), plt.plot(x,y), plt.show() shows the mapping function
- dhdt.preprocessing.image_transforms.log_adjustment(Z)¶
transform intensity in non-linear way
enhances low intensities
- Parameters
Z (numpy.ndarray, size=(m,n)) – grid with intensity values
- Returns
Z_new – grid with transform
- Return type
numpy.ndarray, size=(m,n)
- dhdt.preprocessing.image_transforms.mat_to_gray(Z, notZ=None, vmin=None, vmax=None)¶
transform matix array to float, omitting nodata values
- Parameters
Z (numpy.ndarray, size=(m,n), dtype={integer,float}) – matrix of integers with data
notZ (numpy.ndarray, size=(m,n), dtype=bool) – matrix assigning which is a no data. The default is None.
- Returns
Znew – array with normalized intensity values
- Return type
numpy.ndarray, size=(m,n), dtype=float, range=0…1
See also
gamma_adjustment,log_adjustment,s_curve,inverse_tangent_transformationExamples
>>> import numpy as np >>> from skimage import data >>> img = data.astronaut() >>> np.max(img) 255 >>> gray = mat_to_gray(img[:,:,0], img[:,:,0]==0) >>> np.max(gray) 1.0
- dhdt.preprocessing.image_transforms.multi_spectral_grad(Ims)¶
get the dominant multi-spectral gradient of the stack, stripped down version of the fusion technique, described in [SW02].
- Parameters
Ims (numpy.array, dim=3, dtype=float) – stack of imagery
- Returns
grad_I – fused gradient of the image stack
- Return type
numpy.array, dim=2, dtype=complex
References
- SW02
Socolinsky & Wolff, “Multispectral image visualization through first-order fusion”, IEEE transactions on image fusion, vol.11(8) pp.923-931, 2002.
- dhdt.preprocessing.image_transforms.normalize_histogram(img)¶
transform image to have a uniform distribution
- Parameters
img (numpy.ndarray, size=(m,n), dtype={uint8,unit16}) – image array
- Returns
new_img – transfrormed image
- Return type
numpy.ndarray, size=(m,n), dtype={uint8,unit16}
See also
- dhdt.preprocessing.image_transforms.psuedo_inv(cumhist, val)¶
- Parameters
cumhist (numpy.array, size=(2**n,), dtype={integer, float}) – emperical cumulative distribution function, i.e. a cumulative histogram.
val ({integer,float}) –
- Returns
res – the pseudo-inverse of cumhist
- Return type
{integer,float}
Notes
The following acronyms are used:
CDF - cumulative distribution function
color transforms¶
functions transforming to and from different color spaces
- dhdt.preprocessing.color_transforms.erdas2hsi(Blue, Green, Red)¶
transform red, green, blue arrays to hue, saturation, intensity arrays
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Red (numpy.array, size=(m,n)) – red band of satellite image
- Returns
Hue (numpy.array, size=(m,n), float) – hue
Sat (numpy.array, size=(m,n), float) – saturation
Int (numpy.array, size=(m,n), float) – intensity
See also
References
- ER13
ERDAS, “User handbook”, 2013.
- dhdt.preprocessing.color_transforms.hcv2rgb(H, C, V)¶
transform a color bubble model to red green blue arrays, following [Sh95].
- Parameters
H (numpy.array, size=(m,n)) – array with amount of color
C (numpy.array, size=(m,n)) – luminance
V (numpy.array, size=(m,n)) – array with dominant frequency
- Returns
Red, Green, Blue – array with colours
- Return type
numpy.array, size=(m,n)
References
- Sh95
Shih, “The reversibility of six geometric color spaces”, Photogrammetric engineering & remote sensing, vol.61(10) pp.1223-1232, 1995.
- dhdt.preprocessing.color_transforms.lab2lch(L, a, b)¶
transform XYZ tristimulus arrays to Lab values
- Parameters
L (numpy.array, size=(m,n)) –
a (numpy.array, size=(m,n)) –
b (numpy.array, size=(m,n)) –
- Returns
C,h
- Return type
numpy.array, size=(m,n)
References
- FR98
Ford & Roberts. “Color space conversion”, pp. 1–31, 1998.
- Si18
Silva et al. “Near real-time shadow detection and removal in aerial motion imagery application” ISPRS journal of photogrammetry and remote sensing, vol.140 pp.104–121, 2018.
- dhdt.preprocessing.color_transforms.lms2lab(L, M, S)¶
transform L, M, S arrays to lab color space
- Parameters
L (numpy.array, size=(m,n)) –
M (numpy.array, size=(m,n)) –
S (numpy.array, size=(m,n)) –
- Returns
l,a,b
- Return type
numpy.array, size=(m,n)
References
- Re01
Reinhard et al. “Color transfer between images” IEEE Computer graphics and applications vol.21(5) pp.34-41, 2001.
- dhdt.preprocessing.color_transforms.rgb2hcv(Blue, Green, Red)¶
transform red green blue arrays to a color space
- Parameters
Blue (numpy.array, size=(m,n)) – Blue band of satellite image
Green (numpy.array, size=(m,n)) – Green band of satellite image
Red (numpy.array, size=(m,n)) – Red band of satellite image
- Returns
H (numpy.array, size=(m,n)) – array with amount of color
C (numpy.array, size=(m,n)) – luminance
V (numpy.array, size=(m,n)) – array with dominant frequency
References
- Sm93
Smith, “Putting colors in order”, Dr. Dobb’s Journal, pp 40, 1993.
- Ts06
Tsai, “A comparative study on shadow compensation of color aerial images in invariant color models”, IEEE transactions in geoscience and remote sensing, vol. 44(6) pp. 1661–1671, 2006.
- dhdt.preprocessing.color_transforms.rgb2hsi(Red, Green, Blue)¶
transform red, green, blue arrays to hue, saturation, intensity arrays
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Blue (numpy.array, size=(m,n)) – blue band of satellite image
- Returns
Hue (numpy.array, size=(m,n), range=0…1) – Hue
Sat (numpy.array, size=(m,n), range=0…1) – Saturation
Int (numpy.array, size=(m,n), range=0…1) – Intensity
References
- Ts06
Tsai, “A comparative study on shadow compensation of color aerial images in invariant color models”, IEEE transactions in geoscience and remote sensing, vol. 44(6) pp. 1661–1671, 2006.
- Pr91
Pratt, “Digital image processing” Wiley, 1991.
- dhdt.preprocessing.color_transforms.rgb2lms(Red, Green, Blue)¶
transform red, green, blue arrays to XYZ tristimulus values
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Blue (numpy.array, size=(m,n)) – blue band of satellite image
- Returns
L,M,S
- Return type
numpy.array, size=(m,n)
References
- Re01
Reinhard et al. “Color transfer between images” IEEE Computer graphics and applications vol.21(5) pp.34-41, 2001.
- dhdt.preprocessing.color_transforms.rgb2xyz(Red, Green, Blue, method='reinhardt')¶
transform red, green, blue arrays to XYZ tristimulus values
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Blue (numpy.array, size=(m,n)) – blue band of satellite image
method –
- ‘reinhardt’
XYZitu601-1 axis
- ’ford’
D65 illuminant
- Returns
X,Y,Z
- Return type
numpy.array, size=(m,n)
References
- Re01
Reinhard et al. “Color transfer between images” IEEE Computer graphics and applications vol.21(5) pp.34-41, 2001.
- FR98
Ford & Roberts. “Color space conversion”, pp. 1–31, 1998.
- dhdt.preprocessing.color_transforms.rgb2ycbcr(Red, Green, Blue)¶
transform red, green, blue arrays to luna and chroma values
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Blue (numpy.array, size=(m,n)) – blue band of satellite image
- Returns
Y (numpy.array, size=(m,n)) – luma
Cb (numpy.array, size=(m,n)) – chroma
Cr (numpy.array, size=(m,n)) – chroma
References
- Ts06
Tsai, “A comparative study on shadow compensation of color aerial images in invariant color models”, IEEE transactions in geoscience and remote sensing, vol. 44(6) pp. 1661–1671, 2006.
- dhdt.preprocessing.color_transforms.rgb2yiq(Red, Green, Blue)¶
transform red, green, blue to luminance, inphase, quadrature values
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Blue (numpy.array, size=(m,n)) – blue band of satellite image
- Returns
Y (numpy.array, size=(m,n)) – luminance
I (numpy.array, size=(m,n)) – inphase
Q (numpy.array, size=(m,n)) – quadrature
References
- GW92
Gonzalez & Woods “Digital image processing”, 1992.
- dhdt.preprocessing.color_transforms.xyz2lab(X, Y, Z, th=0.008856)¶
transform XYZ tristimulus arrays to Lab values
- Parameters
X (numpy.array, size=(m,n)) –
Y (numpy.array, size=(m,n)) –
Z (numpy.array, size=(m,n)) –
- Returns
L,a,b
- Return type
numpy.array, size=(m,n)
References
- FR98
Ford & Roberts. “Color space conversion”, pp. 1–31, 1998.
- Si18
Silva et al. “Near real-time shadow detection and removal in aerial motion imagery application” ISPRS journal of photogrammetry and remote sensing, vol.140 pp.104–121, 2018.
- dhdt.preprocessing.color_transforms.xyz2lms(X, Y, Z)¶
transform XYZ tristimulus arrays to LMS values
- Parameters
X (numpy.array, size=(m,n)) – modified XYZitu601-1 axis
Y (numpy.array, size=(m,n)) – modified XYZitu601-1 axis
Z (numpy.array, size=(m,n)) – modified XYZitu601-1 axis
- Returns
L,M,S
- Return type
numpy.array, size=(m,n)
References
- Re01
Reinhard et al. “Color transfer between images” IEEE Computer graphics and applications vol.21(5) pp.34-41, 2001.
- dhdt.preprocessing.color_transforms.yiq2rgb(Y, Inph, Quadr)¶
transform luminance, inphase, quadrature values to red, green, blue
- Parameters
Y (numpy.array, size=(m,n)) – luminance
Inph (numpy.array, size=(m,n)) – inphase
Quadr (numpy.array, size=(m,n)) – quadrature
- Returns
Red (numpy.array, size=(m,n)) – red band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Blue (numpy.array, size=(m,n)) – blue band of satellite image
See also
References
- GW92
Gonzalez & Woods “Digital image processing”, 1992.
acquisition geometry¶
functions of use for geometric information extraction and retrival
- dhdt.preprocessing.acquisition_geometry.compensate_ortho_offset(Img, Z, dx, dy, obs_az, obs_zn, geoTransform)¶
- Parameters
Img (numpy.array, size=(m,n), ndim=2, dtype={float,integer}) – grid with intensities to be compensated
Z (numpy.array, size=(m,n), ndim=2, dtype={float,integer}, unit=meters) – grid with elevations
dx (float, unit=meter) – mis-registration offset
dy (float, unit=meter) – mis-registration offset
obs_az (numpy.array, size=(m,n), unit=degrees) – azimuth angle of observation
obs_zn (numpy.ndarray, size=(m,n), unit=degrees) – zenith angle of observation
geoTransform1 (tuple) – affine transformation coefficients of array Z
- Returns
Img_cor – corrected intensities
- Return type
numpy.array, size=(m,n), ndim=2, dtype=float
See also
- dhdt.preprocessing.acquisition_geometry.get_ortho_offset(Z, dx, dy, obs_az, obs_zn, geoTransform)¶
get terrain displacement due to miss-registration
- Parameters
Z (numpy.array, size=(m,n), unit=meter) – elevation model
dx (float, unit=meter) – mis-registration
dy (float, unit=meter) – mis-registration
obs_az (numpy.array, size=(m,n), unit=degrees) – azimuth angle of observation
obs_zn (numpy.ndarray, size=(m,n), unit=degrees) – zenith angle of observation
geoTransform (tuple, size={(6,), (8,)}) – affine transformation coefficients of array Z
- Returns
dI, dJ – terrain displacements
- Return type
numpy.ndarray, size=(m,n), unit=pixels
Notes
It is important to know what type of coordinate systems exist, hence:
coordinate | coordinate ^ y system 'ij'| system 'xy' | | | | j | x --------+--------> --------+--------> | | | | | i | v |The azimuth angle is declared in the following coordinate frame:
^ North & y | - <--|--> + | +----> East & xThe angles related to the satellite are as follows:
#*# #*# satellite ^ / ^ /| | / | / | nadir |-- zenith angle | / v | / | /| |/ |/ | elevation angle +----- surface +------
- dhdt.preprocessing.acquisition_geometry.get_template_acquisition_angles(Az, Zn, Det, i_samp, j_samp, t_size)¶
- Parameters
Az (numpy.ndarray, size=(m,n), float) – array with sensor azimuth angles
Zn (numpy.ndarray, size=(m,n), float) – array with sensor zenith angles
Det (numpy.ndarray, size=(m,n), bool) – array with sensor detector ids
i_samp (numpy.ndarray, size=(k,l), integer) – array with collumn coordinates of the template centers
j_samp (numpy.ndarray, size=(k,l), integer) – array with row coordinates of the template centers
t_size (integer, {x ∈ ℕ | x ≥ 1}, unit=pixels) – size of the template
- Returns
Azimuth, Zenith – mean observation angles, that is, azimuth and zenith
- Return type
np.array, size=(k,l), float
Notes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--|--> + | +----> East & x
- dhdt.preprocessing.acquisition_geometry.get_template_aspect_slope(Z, i_samp, j_samp, t_size, spac=10.0)¶
- Parameters
Z (numpy.ndarray, size=(m,n), float, unit=meters) – array with elevation values
i_samp (numpy.ndarray, size=(k,l), integer, unit=pixels) – array with collumn coordinates of the template centers
j_samp (numpy.ndarray, size=(k,l), integer, unit=pixels) – array with row coordinates of the template centers
t_size (integer, {x ∈ ℕ | x ≥ 1}, unit=pixels) – size of the template
spac (float, unit=meters) – spacing of the elevation grid
- Returns
Slope, Aspect – mean slope and aspect angle in the template
- Return type
numpy.ndarray, size=(k,l), float
See also
get_aspect_slopeNotes
It is important to know what type of coordinate systems exist, hence:
coordinate | coordinate ^ y system 'ij'| system 'xy' | | | | j | x --------+--------> --------+--------> | | | | | i | v |
- dhdt.preprocessing.acquisition_geometry.slope_along_perp(Z, Az, spac=10)¶
get the slope of the terrain along a specific direction
- Parameters
Z (numpy.ndarray, size=(m,n), units=meters) – elevation model
Az ({numpy.ndarray, float}, units=degrees) – argument or azimuth angle of interest
spac (float, units=meters) – spacing used for the elevation model
- Returns
Slp_para (numpy.ndarray, size=(m,n), units=degrees) – slope in the along direction of the azimuth angle
Slp_perp (numpy.ndarray, size=(m,n), units=degrees) – slope in the perpendicular direction of the azimuth angle
aquatic transforms¶
functions of use for enhancement of water or suspended material
- dhdt.preprocessing.aquatic_transforms.alternative_floating_algae_index(Red, Rededge, Near)¶
transform near and short wave infrared arrays to get a algae-index. Based upon [Wa18]
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Rededge (numpy.ndarray, size=(m,n)) – rededge band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
AFAI – array with floating algae transform
- Return type
numpy.ndarray, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[AFAI=(B_{06} - B_{04})/(B_{8A} - B_{04})\]References
- Wa18
Wang, et al. “Remote sensing of Sargassum biomass, nutrients, and pigments” Geophysical Research Letters, vol.45(22) pp.12-359, 2018.
- dhdt.preprocessing.aquatic_transforms.atmospherically_resistant_vegetation_index(Blue, Red, Near, gamma=1.0, L=0.5)¶
Based upon [KT92]
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
ARVI – array with atmospherically resistant vegetation index
- Return type
numpy.ndarray, size=(m,n)
References
- KT92
Kaufman & Tanre, “Atmospherically resistant vegetation index (ARVI)” IEEE transactions in geoscience and remote sensing, vol.30 pp.261–270, 1992.
- dhdt.preprocessing.aquatic_transforms.floating_algae_index(Red, Near, Swir, df)¶
Following [Hu09]
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
Swir (numpy.ndarray, size=(m,n)) – shortwave infrared band of satellite image
df (datafram) –
- Returns
FAI – array with floating algae index
- Return type
numpy.ndarray, size=(m,n)
References
- Hu09
Hu, “A novel ocean color index to detect floating algae in the global oceans” Remote sensing of the enviroment, vol.113(10) pp.2188-2129, 2009
- dhdt.preprocessing.aquatic_transforms.infrared_percentage_vegetation_index(Red, Near)¶
adjustment over the NDVI, see [Cr90].
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
NPVI – array with infrared percentage vegetation index
- Return type
numpy.ndarray, size=(m,n), range=0…1
References
- Cr90
Crippen “Calculating the vegetation index faster” Remote sensing of environment, vol.34 pp.71-73, 1990.
- dhdt.preprocessing.aquatic_transforms.modified_chlorophyll_absorption_index(Green, Red, Rededge)¶
Based upon [Da00]
- Parameters
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Rededgde (numpy.ndarray, size=(m,n)) – rededge band of satellite image
- Returns
MCARI – array with atmospherically resistant vegetation index
- Return type
numpy.ndarray, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[MCARI= ((B_{05} - B_{04}) - .2(B_{05} - B_{03})) * (B_{04})/(B_{03})\]References
- Da00
Daughtry et al. “Estimating corn leaf chlorophyll concentration from leaf and canopy reflectance” Remote sensing of environment vol.74 pp.229–239, 2000.
- dhdt.preprocessing.aquatic_transforms.modified_simple_ratio(Red, Near)¶
Based upon [Ch96]
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
SR – array with simple ratio
- Return type
numpy.ndarray, size=(m,n)
References
- Ch96
Chen, “Evaluation of vegetation indices and modified simple ratio for boreal applications” Canadian journal of remote sensing, vol.22 pp.229–242, 1996.
- dhdt.preprocessing.aquatic_transforms.modified_soil_adjusted_vegetation_index(Red, Near)¶
Based upon [Qi94]
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
MSAVI – array with modified soil adjusted vegetation index
- Return type
numpy.ndarray, size=(m,n)
References
- Qi94
Qi et al., “A modified soil vegetation adjusted index” Remote sensing of environment, vol.48, pp.119–126, 1994.
- dhdt.preprocessing.aquatic_transforms.normalized_difference_vegetation_index(Red, Near)¶
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
NDVI – array with normalized difference vegetation index
- Return type
numpy.ndarray, size=(m,n)
References
- Ro74
Rouse et al. “Monitoring the vernal advancements and retrogradation of natural vegetation” in NASA/GSFC final report, pp.1–137, 1974.
- dhdt.preprocessing.aquatic_transforms.renormalized_difference_vegetation_index(Red, Near)¶
Based upon [RB95]
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
NDVI – array with renormalized difference vegetation index
- Return type
numpy.ndarray, size=(m,n)
References
- RB95
Rougean & Breon, “Estimating PAR absorbed by vegetation from bidirectional reflectance measurements” Remote sensing of environment, vol.51 pp.375–384, 1995.
- dhdt.preprocessing.aquatic_transforms.simple_ratio(Red, Near)¶
Based upon [Jo69]
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
SR – array with simple ratio
- Return type
numpy.ndarray, size=(m,n)
References
- Jo69
Jordan, “Derivation of leaf area index from quality of light on the forest floor” Ecology vol.50 pp.663–666, 1969.
- dhdt.preprocessing.aquatic_transforms.soil_adjusted_vegetation_index(Red, Near, L=0.5)¶
Based upon [Hu88]
- Parameters
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
L (float, default=.5) – multiplication factor
- Returns
SAVI – array with soil adjusted vegetation index
- Return type
numpy.ndarray, size=(m,n)
References
- Hu88
Huete, “A soil vegetation adjusted index (SAVI)” Remote sensing of environment, vol.25 pp.295–309, 1988.
- dhdt.preprocessing.aquatic_transforms.triangular_vegetation_index(Green, Red, Rededge)¶
Based upon [BL00]
- Parameters
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Rededgde (numpy.ndarray, size=(m,n)) – rededge band of satellite image
- Returns
TVI – array with triangular vegetation index
- Return type
numpy.ndarray, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[TVI= .5 (120(B_{06} - B_{03}) - 200(B_{04} - B_{03}))\]References
- BL00
Broge & Leblanc, “Comparing prediction power and stability of broadband and hyperspectral vegetation indices for estimation of green leaf area index and canopy chlorophyll density” Remote sensing of environment, vol.76 pp.156–172, 2000.
snow transforms¶
functions of use for enhancement of snow and ice
- dhdt.preprocessing.snow_transforms.automated_glacier_extraction_index(Red, Near, Short, alpha=0.5)¶
transform red, near and shortwave infrared arrays to AGEI, based on [Zh19].
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Near (numpy.array, size=(m,n)) – near infrared band of satellite image
Short (numpy.array, size=(m,n)) – shortwave infrared band of satellite image
alpha (float) –
- Returns
AGEI – array with AGEI transform
- Return type
numpy.array, size=(m,n)
See also
Notes
Based on the bands of Sentinel-2:
\[AGEI = (B_{3}-B_{8}) / (B_{3}+B_{8})\]References
- Zh19
Zhang et al. “Automated glacier extraction index by optimization of red/SWIR and NIR/SWIR ratio index for glacier mapping using Landsat imagery” Water, vol.11 pp.1233, 2019.
- dhdt.preprocessing.snow_transforms.normalized_difference_snow_ice_index(Red, Short)¶
transform red and short wave infrared arrays NDSII-index, see also [Xi01]
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Short (numpy.array, size=(m,n)) – shortwave infrared band of satellite image
- Returns
NDSII – array with NDSII transform
- Return type
numpy.array, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[NDSII = (B_{4}-B_{11}) / (B_{4}+B_{11})\]The bands of Sentinel-2 correspond to the following spectral range:
\(B_{4}\) : red
\(B_{11}\) : shortwave infrared
References
- Xi01
Xiao et al. “Assessing the potential of Vegetation sensor data for mapping snow and ice cover: A normalized snow and ice index” International journal of remote sensing, vol.22(13) pp.2479-2487, 2001.
- dhdt.preprocessing.snow_transforms.normalized_difference_snow_index(Green, Short)¶
transform red and short wave infrared arrays NDSII-index, see also [Do89].
- Parameters
Green (numpy.array, size=(m,n)) – green band of satellite image
Short (numpy.array, size=(m,n)) – shortwave infrared band of satellite image
- Returns
NDSI – array with NDSI transform
- Return type
numpy.array, size=(m,n)
Notes
- Based on the bands of Sentinel-2:
- \[NDSI = (B_{3}-B_{8}) / (B_{3}+B_{8})\]
- The bands of Sentinel-2 correspond to the following spectral range:
\(B_{3}\) : green
\(B_{8}\) : near infrared
See also
References
- Do89
Dozier. “Spectral signature of alpine snow cover from the Landsat Thematic Mapper” Remote sensing of environment, vol.28 pp.9-22, 1989.
- dhdt.preprocessing.snow_transforms.s3(Red, Near, Short)¶
transform red, near and shortwave infrared arrays to s3, see [SY99].
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Near (numpy.array, size=(m,n)) – near infrared band of satellite image
Short (numpy.array, size=(m,n)) – shortwave infrared band of satellite image
- Returns
S3 – array with S3 transform
- Return type
numpy.array, size=(m,n)
See also
Notes
Based on the bands of Sentinel-2:
\[S3 = B_{8}(B_{4}-B_{11}) / ((B_{4} + B_{8})(B_{4} + B_{11}))\]The bands of Sentinel-2 correspond to the following spectral range:
\(B_{4}\) : red
\(B_{8}\) : near infrared
\(B_{11}\) : shortwave infrared
See also
References
- SY99
Saito & Yamazaki. “Characteristics of spectral reflectance for vegetation ground surfaces with snowcover; Vegetation indeces and snow indices” Journal of Japan society of hydrology and water resources, vol.12(1) pp.28-38, 1999.
- dhdt.preprocessing.snow_transforms.snow_cover_fraction(Near, Short)¶
transform near and short wave infrared arrays to get a snow-index, from [Pa02].
- Parameters
Near (numpy.array, size=(m,n)) – near infrared band of satellite image
Short (numpy.array, size=(m,n)) – shortwave infrared band of satellite image
- Returns
NDSII – array with NDSII transform
- Return type
numpy.array, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[SCF= B_{8} / B_{11}\]The bands of Sentinel-2 correspond to the following spectral range:
\(B_{8}\) : near infrared
\(B_{11}\) : shortwave infrared
References
- Pa02
Paul et al. “The new remote-sensing-derived Swiss glacier inventory: I methods” Annuals of glaciology, vol.34 pp.355-361, 2002.
- dhdt.preprocessing.snow_transforms.snow_fraction(Red, Short)¶
transform red and short wave infrared arrays to get a snow-index, see also [PK05].
- Parameters
Red (numpy.array, size=(m,n)) – red band of satellite image
Short (numpy.array, size=(m,n)) – shortwave infrared band of satellite image
- Returns
SF – array with snow fraction transform
- Return type
numpy.array, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[SF = B_{4} / B_{11}\]The bands of Sentinel-2 correspond to the following spectral range:
\(B_{4}\) : red
\(B_{11}\) : shortwave infrared
References
- PK05
Paul & Kääb “Perspectives on the production of a glacier inventory from multispectral satellite data in Arctic Canada: Cumberland peninsula, Baffin Island” Annuals of glaciology, vol.42 pp.59-66, 2005.
- dhdt.preprocessing.snow_transforms.snow_water_index(Green, Near, Short)¶
transform green, near and shortwave infrared arrays to snow water index, as described in [Di19].
- Parameters
Green (numpy.array, size=(m,n)) – green band of satellite image
Near (numpy.array, size=(m,n)) – near infrared band of satellite image
Short (numpy.array, size=(m,n)) – shortwave infrared band of satellite image
- Returns
SWI – array with SWI transform
- Return type
numpy.array, size=(m,n)
See also
Notes
Based on the bands of Sentinel-2:
\[SWI = B_{8}(B_{3}-B_{11}) / ((B_{3} + B_{8})(B_{3} + B_{11}))\]The bands of Sentinel-2 correspond to the following spectral range:
\(B_{3}\) : green
\(B_{8}\) : near infrared
\(B_{11}\) : shortwave infrared
References
- Di19
Dixit et al. “Development and evaluation of a new snow water index (SWI) for accurate snow cover delineation” Remote sensing, vol.11(23) pp.2774, 2019.
shadow transforms¶
functions of use for enhancement of intensities, making use of spectral properties
- dhdt.preprocessing.shadow_transforms.apply_shadow_transform(method, Blue, Green, Red, RedEdge, Near, Shw, **kwargs)¶
Given a specific method, employ shadow transform
- Parameters
method ({‘siz’,’ isi’,’ sei’,’ fcsdi’, ‘nsvi’, ‘nsvdi’, ‘sil’, ‘sr’, ‘sdi’, ‘sisdi’, ‘mixed’, ‘c3’, ‘entropy’, ‘shi’,’ sp’, 'nri'}) –
method name to be implemented,can be one of the following:
’siz’ : shadow index first version
’isi’ : improved shadow index
’sei’ : shadow enhancement index
’fcsdi’ : false color shadow difference index
’csi’ : combinational shadow index
’nsvdi’ : normalized saturation value difference index
’sil’ : shadow index second
’sr’ : specthem ratio
’sdi’ : shadow detector index
’sisdi’ : saturation intensity shadow detector index
’mixed’ : mixed property based shadow index
’msf’ : modified shadow fraction
’c3’ : color invariant
’entropy’ : entropy shade removal
’shi’ : shade index
’sp’ : shadow probabilities
’nri’ : normalized range shadow index
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
RedEdge (numpy.ndarray, size=(m,n)) – red edge band of satellite image
Near (numpy.ndarray, size=(m,n)) – near-infrared band of satellite image
Shw (numpy.ndarray, size=(m,n)) – synthetic shading from for example a digital elevation model
- Returns
M – shadow enhanced satellite image
- Return type
numpy.ndarray, size=(m,n)
See also
handler_multispec.get_shadow_bandsprovides Blue, Green, Red bandnumbers for a specific satellite
handler_multispec.read_shadow_bandsreads the specific bands given
- dhdt.preprocessing.shadow_transforms.color_invariant(Blue, Green, Red)¶
transform red, green, blue arrays to color invariant (c3), see also [GS99] and [AA08].
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Red (numpy.array, size=(m,n)) – red band of satellite image
- Returns
c3 – array with shadow enhancement
- Return type
numpy.array, size=(m,n)
References
- GS99
Gevers & Smeulders. “Color-based object recognition” Pattern recognition, vol.32 pp.453–464, 1999
- AA08
Arévalo González & Ambrosio. “Shadow detection in colour high‐resolution satellite images” International journal of remote sensing, vol.29(7) pp.1945–1963, 2008
- dhdt.preprocessing.shadow_transforms.combinational_shadow_index(Blue, Green, Red, Near)¶
transform red, green, blue arrays to combined shadow index. See also [Su19].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
CSI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
See also
References
- Su19
Sun et al. “Combinational shadow index for building shadow extraction in urban areas from Sentinel-2A MSI imagery” International journal of applied earth observation and geoinformation, vol.78 pp.53–65, 2019.
- dhdt.preprocessing.shadow_transforms.entropy_shade_removal(Ia, Ib, Ic, a=None)¶
Make use of Wiens’ law to get illumination component, see also [Fi09].
- Parameters
Ia (numpy.array, size=(m,n)) – highest frequency band of satellite image
Ib (numpy.array, size=(m,n)) – lower frequency band of satellite image
Ic (numpy.array, size=(m,n)) – lowest frequency band of satellite image
a (float, default=-1, unit=degrees) – angle of the coordinate rotation
- Returns
S – shadow enhanced band
- Return type
numpy.array, size=(m,n)
References
- Fi09
Finlayson et al. “Entropy minimization for shadow removal” International journal of computer vision vol.85(1) pp.35-57 2009.
- dhdt.preprocessing.shadow_transforms.false_color_shadow_difference_index(Green, Red, Near)¶
transform red and near infrared arrays to shadow index. See also [Te11].
- Parameters
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
FCSI – array with shadow enhancement
- Return type
numpy.ndarray, size=(m,n)
See also
Notes
\[FCSI = (I_{sat} - I_{int}) / (I_{sat} + I_{int})\]References
- Te11
Teke et al. “Multi-spectral false color shadow detection”, Proceedings of the ISPRS conference on photogrammetric image analysis, pp.109-119, 2011.
- dhdt.preprocessing.shadow_transforms.fractional_range_shadow_index(*args)¶
transform multi-spectral arrays to shadow index
- Parameters
*args (numpy.array, size=(m,n)) – different bands of a satellite image, all of the same extent
- Returns
SI – shadow enhanced band
- Return type
numpy.array, size=(m,n)
Notes
Schematically this looks like:
I1 ┌-----┐ -┬---┤ | I2|┌--┤ | Imin -┼┼┬-┤ min ├------┐ |||┌┤ | |┌-------┐ ||||└-----┘ └┤ a / b | SI I3||||┌-----┐ ┌┤ ├------ -┼┴┼┼┤ | Imax |└-------┘ I4└-┼┼┤ ├------┘ ----┼┴┤ max | └-┤ | └-----┘
- dhdt.preprocessing.shadow_transforms.improved_shadow_index(Blue, Green, Red, Near)¶
transform red, green, blue arrays to improved shadow index. Based upon [Zh21].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
ISI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
See also
Notes
\[ISI = (SI + (1-I_{near}))/(SI + (1+I_{near}))\]References
- Zh21
Zhou et al. “Shadow detection and compensation from remote sensing images under complex urban conditions”, 2021.
- dhdt.preprocessing.shadow_transforms.mixed_property_based_shadow_index(Blue, Green, Red)¶
transform red, green, blue arrays to Mixed property-based shadow index, see also [Ha18].
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Red (numpy.array, size=(m,n)) – red band of satellite image
- Returns
MPSI – array with shadow enhancement
- Return type
numpy.array, size=(m,n)
References
- Ha18
Han et al. “A mixed property-based automatic shadow detection approach for VHR multispectral remote sensing images” Applied sciences vol.8(10) pp.1883, 2018.
- dhdt.preprocessing.shadow_transforms.modified_color_invariant(Blue, Green, Red, Near)¶
transform red, green, blue arrays to color invariant (c3), see also [GS99] and [BA15].
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Red (numpy.array, size=(m,n)) – red band of satellite image
Near (numpy.array, size=(m,n)) – near infrared band of satellite image
- Returns
c3 – array with shadow enhancement
- Return type
numpy.array, size=(m,n)
References
- GS99
Gevers & Smeulders. “Color-based object recognition” Pattern recognition, vol.32 pp.453–464, 1999
- BA15
Besheer & Abdelhafiz. “Modified invariant colour model for shadow detection” International journal of remote sensing, vol.36(24) pp.6214–6223, 2015.
- dhdt.preprocessing.shadow_transforms.modified_normalized_difference_water_index(Green, Shortwave)¶
transform green and shortwave infrared arrays to modified NDW-index. See also [Xu06].
- Parameters
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Shortwave (numpy.ndarray, size=(m,n)) – shortwave infrared band of satellite image, for Sentinel-2 this is B11, originally the Landsat TM band 5 was used in [Xu06].
- Returns
NDWI – array with NDWI transform
- Return type
numpy.ndarray, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[NDWI = [B_{3}-B_{11}]/[B_{3}+B_{11}]\]References
- dhdt.preprocessing.shadow_transforms.modified_shadow_fraction(Blue, Green, Red, P_S=0.95)¶
transform red, green, blue arrays to shadow index. Based upon [Ch09].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
P_S (float, range=0...1) – threshold value for classification, based on the quantile
- Returns
SI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
See also
References
- Ch09
Chung, et al. “Efficient shadow detection of color aerial images based on successive thresholding scheme.” IEEE transactions on geoscience and remote sensing vol.47, pp.671–682, 2009.
- dhdt.preprocessing.shadow_transforms.normalized_difference_blue_water_index(Blue, Near)¶
transform green and near infrared arrays NDW-index. See also [Qu11].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
NDWI – array with NDWI transform
- Return type
numpy.ndarray, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[NDWI = [B_{3}-B_{8}]/[B_{3}+B_{8}]\]References
- Qu11
Qu et al, “Research on automatic extraction of water bodies and wetlands on HJU satellite CCD images” Remote sensing information, vol.4 pp.28–33, 2011
- dhdt.preprocessing.shadow_transforms.normalized_difference_moisture_index(Near, Shortwave)¶
transform near and shortwave infrared arrays to NDM-index. See also [Wi02].
- Parameters
Near (numpy.ndarray, size=(m,n)) – green band of satellite image
Shortwave (numpy.ndarray, size=(m,n)) – shortwave infrared band of satellite image
- Returns
NDMI – array with NDMI transform
- Return type
numpy.ndarray, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[NDMI = [B_{8}-B_{11}]/[B_{8}+B_{11}]\]References
- Wi02
Wilson, “Detection of forest harvest type using multiple dates of Landsat TM imagery” Remote sensing of environment, vol.80 pp.385–396, 2002.
- dhdt.preprocessing.shadow_transforms.normalized_difference_water_index(Green, Near)¶
transform green and near infrared arrays NDW-index. See also [Ga96].
- Parameters
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
NDWI – array with NDWI transform
- Return type
numpy.ndarray, size=(m,n)
Notes
Based on the bands of Sentinel-2:
\[NDWI = [B_{3}-B_{8}]/[B_{3}+B_{8}]\]References
- Ga96
Gao, “NDWI - a normalized difference water index for remote sensing of vegetation liquid water from space” Remote sensing of environonment, vol.58 pp.257–266, 1996
- dhdt.preprocessing.shadow_transforms.normalized_range_shadow_index(*args)¶
transform multi-spectral arrays to shadow index
- Parameters
*args (numpy.array, size=(m,n)) – different bands of a satellite image, all of the same extent
- Returns
SI – shadow enhanced band
- Return type
numpy.array, size=(m,n)
Notes
Schematically this looks like:
I1 ┌-----┐ -┬---┤ | I2|┌--┤ | Imin -┼┼┬-┤ min ├------┐ |||┌┤ | |┌-------┐ ||||└-----┘ └┤(a-b)/ | SI I3||||┌-----┐ ┌┤ (a+b)├------ -┼┴┼┼┤ | Imax |└-------┘ I4└-┼┼┤ ├------┘ ----┼┴┤ max | └-┤ | └-----┘
- dhdt.preprocessing.shadow_transforms.normalized_sat_value_difference_index(Blue, Green, Red)¶
transform red, green, blue arrays to normalized sat-value difference. See also [Ma08].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
NSVDI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
See also
false_color_shadow_difference_index,rgb2hsiNotes
\[NSVDI = (I_{sat} - I_{int}) / (I_{sat} + I_{int})\]References
- Ma08
Ma et al. “Shadow segmentation and compensation in high resolution satellite images”, Proceedings of IEEE IGARSS, pp.II-1036–II-1039, 2008.
- dhdt.preprocessing.shadow_transforms.reinhard(Blue, Green, Red)¶
transform red, green, blue arrays to luminance, see also [Re01].
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Red (numpy.array, size=(m,n)) – red band of satellite image
- Returns
l – shadow enhanced band
- Return type
numpy.array, size=(m,n)
References
- Re01
Reinhard et al. “Color transfer between images” IEEE Computer graphics and applications vol.21(5) pp.34-41, 2001.
- dhdt.preprocessing.shadow_transforms.sat_int_shadow_detector_index(Blue, RedEdge, Near)¶
transform red, green, blue arrays to sat-int shadow detector index, see also [MA18].
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
RedEdge (numpy.array, size=(m,n)) – rededge band of satellite image
Near (numpy.array, size=(m,n)) – near infrared band of satellite image
- Returns
SISDI – array with shadow enhancement
- Return type
numpy.array, size=(m,n)
References
- MA18
Mustafa & Abedelwahab. “Corresponding regions for shadow restoration in satellite high-resolution images” International journal of remote sensing, vol.39(20) pp.7014–7028, 2018.
- dhdt.preprocessing.shadow_transforms.shade_index(*args)¶
transform multi-spectral arrays to shade index, see also [HS10] and [Al12].
- Parameters
*args (numpy.array, size=(m,n)) – different bands of a satellite image, all of the same extent
- Returns
SI – shadow enhanced band
- Return type
numpy.array, size=(m,n)
Notes
Amplify dark regions by simple multiplication:
\[SI = Pi_{i}^{N}{[1 - B_{i}]}\]References
- HS10
Hogan & Smith, “Refinement of digital elvation models from shadowing cues” IEEE computer society conference on computer vision and pattern recognition, 2010.
- Al12
Altena, “Filling the white gap on the map: Photoclinometry for glacier elevation modelling” MSc thesis TU Delft, 2012.
- dhdt.preprocessing.shadow_transforms.shadow_detector_index(Blue, Green, Red)¶
transform red, green, blue arrays to shadow detector index, see [MA17].
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Red (numpy.array, size=(m,n)) – red band of satellite image
- Returns
SDI – array with shadow enhancement
- Return type
numpy.array, size=(m,n)
References
- MA17
Mustafa & Abedehafez. “Accurate shadow detection from high-resolution satellite images” IEEE geoscience and remote sensing letters, vol.14(4) pp.494–498, 2017.
- dhdt.preprocessing.shadow_transforms.shadow_enhancement_index(Blue, Green, Red, Near)¶
transform red, green, blue arrays to SEI index. See also [Su19].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
Near (numpy.ndarray, size=(m,n)) – near infrared band of satellite image
- Returns
SEI – array with sei transform
- Return type
numpy.ndarray, size=(m,n)
See also
Notes
Based on the bands of Sentinel-2:
\[SEI = [(B_{1}+B_{9})-(B_{3}+B_{8})]/[(B_{1}+B_{9})+(B_{3}+B_{8})]\]References
- Su19
Sun et al. “Combinational shadow index for building shadow extraction in urban areas from Sentinel-2A MSI imagery” International journal of applied earth observation and geoinformation, vol.78 pp.53–65, 2019.
- dhdt.preprocessing.shadow_transforms.shadow_hcv_fraction(Blue, Green, Red)¶
transform red, green, blue arrays to shadow index. Based upon [Ts06].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
SI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
References
- Ts06
Tsai. “A comparative study on shadow compensation of color aerial images in invariant color models.” IEEE Transactions on geoscience and remote sensing vol.44 pp.1661–1671, 2006.
- dhdt.preprocessing.shadow_transforms.shadow_hsv_fraction(Blue, Green, Red)¶
transform red, green, blue arrays to shadow index. Based upon [Ts06].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
SI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
References
- Ts06
Tsai. “A comparative study on shadow compensation of color aerial images in invariant color models.” IEEE Transactions on geoscience and remote sensing vol.44 pp.1661–1671, 2006.
- dhdt.preprocessing.shadow_transforms.shadow_identification(Blue, Green, Red)¶
transform red, green, blue arrays to sat-value difference, following [Po03].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
SI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
See also
false_color_shadow_difference_index,rgb2hsiNotes
\[SI = I_{sat} - I_{int}\]References
- Po03
Polidorio et al. “Automatic shadow segmentation in aerial color images”, Proceedings of the 16th Brazilian symposium on computer graphics and image processing, pp.270-277, 2003.
- dhdt.preprocessing.shadow_transforms.shadow_index_liu(Blue, Green, Red)¶
transform red, green, blue arrays to shadow index, see [LX13] and [Su19].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
SI – array with shadow enhancement
- Return type
numpy.ndarray, size=(m,n)
References
- LX13
Liu & Xie. “Study on shadow detection in high resolution remote sensing image of PCA and HIS model”, Remote sensing technology and application, vol.28(1) pp. 78–84, 2013.
- Su19
Sun et al. “Combinational shadow index for building shadow extraction in urban areas from Sentinel-2A MSI imagery” International journal of applied earth observation and geoinformation, vol.78 pp.53–65, 2019.
- dhdt.preprocessing.shadow_transforms.shadow_index_zhou(Blue, Green, Red)¶
transform red, green, blue arrays to shadow index, based upon [Zh21].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
SI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
See also
References
- Zh21
Zhou et al. “Shadow detection and compensation from remote sensing images under complex urban conditions”, 2021.
- dhdt.preprocessing.shadow_transforms.shadow_probabilities(Blue, Green, Red, Near, **kwargs)¶
transform blue, green, red and near infrared to shade probabilities, see also [FS10] and [Rü13].
- Parameters
Blue (numpy.array, size=(m,n)) – blue band of satellite image
Green (numpy.array, size=(m,n)) – green band of satellite image
Red (numpy.array, size=(m,n)) – red band of satellite image
Near (numpy.array, size=(m,n)) – near infrared band of satellite image
- Returns
M – shadow enhanced band
- Return type
numpy.array, size=(m,n)
References
- FS10
Fredembach & Süsstrunk. “Automatic and accurate shadow detection from (potentially) a single image using near-infrared information” EPFL Tech Report 165527, 2010.
- Rü13
Rüfenacht et al. “Automatic and accurate shadow detection using near-infrared information” IEEE transactions on pattern analysis and machine intelligence, vol.36(8) pp.1672-1678, 2013.
- dhdt.preprocessing.shadow_transforms.shadow_ycbcr_fraction(Blue, Green, Red)¶
transform red, green, blue arrays to shadow index. Based upon [Ts06].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
SI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
References
- Ts06
Tsai. “A comparative study on shadow compensation of color aerial images in invariant color models.” IEEE Transactions on geoscience and remote sensing vol.44 pp.1661–1671, 2006.
- dhdt.preprocessing.shadow_transforms.shadow_yiq_fraction(Blue, Green, Red)¶
transform red, green, blue arrays to shadow index. Based upon [Ts06].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
SI – array with shadow transform
- Return type
numpy.ndarray, size=(m,n)
See also
References
- Ts06
Tsai. “A comparative study on shadow compensation of color aerial images in invariant color models.” IEEE Transactions on geoscience and remote sensing vol.44 pp.1661–1671, 2006.
- dhdt.preprocessing.shadow_transforms.shannon_entropy(X, band_width=100)¶
calculate shanon entropy from data sample
- Parameters
X (numpy.ndarray, size=(m,1)) – data array
band_width (integer, optional) – sample bandwith
- Returns
shannon – entropy value of the total data sample
- Return type
float
- dhdt.preprocessing.shadow_transforms.specthem_ratio(Blue, Green, Red)¶
transform red, green, blue arrays to specthem ratio, see [Si18].
- Parameters
Blue (numpy.ndarray, size=(m,n)) – blue band of satellite image
Green (numpy.ndarray, size=(m,n)) – green band of satellite image
Red (numpy.ndarray, size=(m,n)) – red band of satellite image
- Returns
Sr – array with shadow enhancement
- Return type
numpy.ndarray, size=(m,n)
Notes
In the paper the logarithmic is also used further to enhance the shadows,
\[Sr = log{[Sr+1]}\]References
- Si18
Silva et al. “Near real-time shadow detection and removal in aerial motion imagery application” ISPRS journal of photogrammetry and remote sensing, vol.140 pp.104–121, 2018.
shadow filters¶
functions of use for enhancement of intensities, making use of spatial properties
- dhdt.preprocessing.shadow_filters.L0_smoothing(Z, lamb=0.02, kappa=2.0, beta_max=100000.0)¶
- Parameters
Z (numpy.ndarray, size={(m,n), (m,n,b)}, dtype=float) – grid with intensity values
lamb (float, default=.002) –
kappa (float, default=2.) –
beta_max (float, default=10000) –
- Returns
Z – modified intensity image
- Return type
numpy.ndarray, size={(m,n), (m,n,b)}, dtype=float
References
- Xu11
Xu et al. “Image smoothing via L0 gradient minimization” ACM transactions on graphics, 2011.
- dhdt.preprocessing.shadow_filters.anistropic_diffusion_scalar(Z, iter=10, K=0.15, s=0.25, n=4)¶
non-linear anistropic diffusion filter of a scalar field
- Parameters
Z (np.array, size=(m,n), dtype={float,complex}, ndim={2,3}) – intensity array
iter (integer, {x ∈ ℕ | x ≥ 0}) – amount of iterations
K (float, default=.15) – parameter based upon the noise level
s (float, default=.25) – parameter for the integration constant
n ({4,8}) – amount of neighbors used
- Returns
I
- Return type
np.array, size=(m,n)
References
- PM87
Perona & Malik “Scale space and edge detection using anisotropic diffusion” Proceedings of the IEEE workshop on computer vision, pp.16-22, 1987
- Ge92
Gerig et al. “Nonlinear anisotropic filtering of MRI data” IEEE transactions on medical imaging, vol.11(2), pp.221-232, 1992
- dhdt.preprocessing.shadow_filters.diffusion_strength_1(Z, K)¶
first diffusion function, proposed by [PM87], when a complex array is provided the absolute magnitude is used, following [Ge92]
- Parameters
Z (np.array, size=(m,n), dtype={float,complex}, ndim={2,3}) – array with intensities
K (float) – parameter based upon the noise level
- Returns
g – array with diffusion parameters
- Return type
np.array, size=(m,n), dtype=float
References
- PM87
Perona & Malik “Scale space and edge detection using anisotropic diffusion” Proceedings of the IEEE workshop on computer vision, pp.16-22, 1987
- Ge92
Gerig et al. “Nonlinear anisotropic filtering of MRI data” IEEE transactions on medical imaging, vol.11(2), pp.221-232, 1992
- dhdt.preprocessing.shadow_filters.diffusion_strength_2(Z, K)¶
second diffusion function, proposed by [PM87], when a complex array is provided the absolute magnitude is used, following [Ge82].
- Parameters
Z (np.array, size=(m,n), dtype={float,complex}) – array with intensities
K (float) – parameter based upon the noise level
- Returns
g – array with diffusion parameters
- Return type
np.array, size=(m,n), dtype=float
References
- PM87
Perona & Malik “Scale space and edge detection using anisotropic diffusion” Proceedings of the IEEE workshop on computer vision, pp.16-22, 1987.
- Ge92
Gerig et al. “Nonlinear anisotropic filtering of MRI data” IEEE transactions on medical imaging, vol.11(2), pp.221-232, 1992.
- dhdt.preprocessing.shadow_filters.enhance_shadows(Shw, method, **kwargs)¶
Given a specific method, employ shadow transform
- Parameters
Shw (np.array, size=(m,n), dtype={float,integer}) – array with intensities of shading and shadowing
method ({‘mean’,’kuwahara’,’median’,’otsu’,'anistropic'}) –
method name to be implemented,can be one of the following:
’mean’ : mean shift filter
’kuwahara’ : kuwahara filter
’median’ : iterative median filter
’anistropic’ : anistropic diffusion filter
- Returns
M – shadow enhanced image, done through the given method
- Return type
np.array, size=(m,n), dtype={float,integer}
- dhdt.preprocessing.shadow_filters.fade_shadow_cast(Shw, az, t_size=9)¶
- Parameters
Shw (np.array, size=(m,n)) – image with mask of shadows
az (float, unit=degrees) – illumination orientation
t_size (integer, {x ∈ ℕ | x ≥ 1}) – buffer size of the Gaussian blur
- Returns
Shw – image with faded shadows on the occluding end
- Return type
np.array, size=(m,n), dtype=float, range=0…1
- dhdt.preprocessing.shadow_filters.iterative_median_filter(Z, tsize=5, loop=50)¶
Transform intensity to more clustered intensity, through iterative filtering with a median operation
- Parameters
Z (np.array, size=(m,n), dtype=float) – array with intensity values
tsize (integer, {x ∈ ℕ | x ≥ 1}) – dimension of the kernel
loop (integer, {x ∈ ℕ | x ≥ 0}) – amount of iterations
- Returns
I_new – array with stark edges
- Return type
np.array, size=(m,n), dtype=float
See also
- dhdt.preprocessing.shadow_filters.kuwahara_filter(Z, tsize=5)¶
Transform intensity to more clustered intensity
- Parameters
Z (np.array, size=(m,n), dtype=float) – array with intensity values
tsize (integer, {x ∈ ℕ | x ≥ 1}) – dimension of the kernel
- Returns
Z_new
- Return type
np.array, size=(m,n), dtype=float
Notes
The template is subdivided into four blocks, where the sampled mean and standard deviation are calculated. Then mean of the block with the lowest variance is assigned to the central pixel. The configuration of the blocks look like:
┌-----┬-┬-----┐ ^ | A | | B | | ├-----┼-┼-----┤ | tsize ├-----┼-┼-----┤ | | C | | D | | └-----┴-┴-----┘ v
See also
- dhdt.preprocessing.shadow_filters.mean_shift_filter(Z, quantile=0.1)¶
Transform intensity to more clustered intensity, through mean-shift
- Parameters
Z (np.array, size=(m,n), dtype=float) – array with intensity values
quantile (float, range=0...1) –
- Returns
labels – array with numbered labels
- Return type
np.array, size=(m,n), dtype=integer
composition and geometry of the atmosphere¶
functions of use for describing the refractive properties of the atmosphere
- dhdt.preprocessing.atmospheric_geometry.ciddor_eq5(ρ_a, ρ_w, n_axs, n_ws, ρ_axs, ρ_ws)¶
estimating the refractive index through summation
- Parameters
ρ_a ({numpy.array, float}, unit=kg m-3) – density of dry air
ρ_w ((numpy.array, float}, unit=kg m-3) – density of moist air
n_axs (numpy.array, size=(k,)) – refractive index for dry air at sea level at a specific wavelength
n_ws (numpy.array, size=(k,)) – refractive index for standard moist air at sea level for a specific wavelength
ρ_axs (float, unit=kg m-3) – density of standard dry air at sea level
ρ_ws (float, unit=kg m-3) – density of standard moist air at sea level
- Returns
n_prop – refractive index
- Return type
float
References
- Ow67
Owens, “Optical refractive index of air: Dependence on pressure, temperature and composition”, Applied optics, vol.6(1) pp.51-59, 1967.
- dhdt.preprocessing.atmospheric_geometry.ciddor_eq6(ρ_a, ρ_w, n_axs, n_ws, ρ_axs, ρ_ws)¶
estimating the refractive index through the Lorenz-Lorentz relation
- Parameters
ρ_a (float, unit=kg m-3) – density of dry air
ρ_w (float, unit=kg m-3) –
n_axs (float) – refractive index for dry air at sea level
n_ws (float) – refractive index for standard moist air at sea level
ρ_axs (float, unit=kg m-3) – density of standard dry air at sea level
ρ_ws (float, unit=kg m-3) – density of standard moist air at sea level
- Returns
n_LL – refractive index
- Return type
float
References
- Ow67
Owens, “Optical refractive index of air: Dependence on pressure, temperature and composition”, Applied optics, vol.6(1) pp.51-59, 1967.
- Ci96
Ciddor, “Refractive index of air: new equations for the visible and near infrared”, Applied optics, vol.35(9) pp.1566-1573, 1996.
- dhdt.preprocessing.atmospheric_geometry.compressability_moist_air(p, T, x_w)¶
- Parameters
p (float, unit=Pascal) – total pressure
T (float, unit=Kelvin) – temperature
x_w (float) – molar fraction of water vapor in moist air
- Returns
Z – compressability factor
- Return type
float
See also
References
- Ci96
Ciddor, “Refractive index of air: new equations for the visible and near infrared”, Applied optics, vol.35(9) pp.1566-1573, 1996.
- DaXX
Davis, “Equation for the determination of the density of moist air (1981/91)” Metrologica, vol.29 pp.67-70.
- dhdt.preprocessing.atmospheric_geometry.get_CO2_level(lat, date, m=3, nb=3)¶
Rough estimation of CO₂, based upon measurements of global ‘Station data’__
The annual trend is based upon the ‘Global trend’__ in CO₂ while the seasonal trend is based upon relations as shown in Figure 1 in [Ba16]. Where a latitudal component is present in the ampltitude, as well as, an asymmetric triangular wave.
- Parameters
lat (float, unit=degrees, range=-90...+90) – latitude of the location of interest
date (np.datetime64, range=1950...2030) – time of interest
m (float, default=3) – assymetry of the wave - m<=2: crest is leftwards leaning - m==2: symmetric wave - m>=2: the ramp is leftwards leaning
nb (integer, default=3) – number of base functions to use, higher numbers will give a more pointy wave function
- Returns
total_co2 – estimate of spatial temporal carbon dioxide (CO₂) value
- Return type
float, unit=ppm
Notes
Data where this fitting is based on, can be found here:
The following acronyms are used:
CO2 : carbon dioxide
ppm : parts per million
Examples
>>> import numpy as np >>> import matplotlib.pyplot as plt
Plot the annual trend at 60ºN
>>> t = np.arange("1950-01-01","2022-01-01",dtype="M8[Y]") >>> co2 = get_CO2_level(60., t) >>> plt.plot(t, co2)
Plot seasonal signal
>>> t = np.arange("2018-01-01","2022-01-01",dtype="M8[D]") >>> co2 = get_CO2_level(60., t) >>> plt.plot(t, co2)
References
- Ba16
Barnes et al. “Isentropic transport and the seasonal cycle amplitude of CO₂” Journal of geophysical research: atmosphere, vol.121 pp.8106-8124, 2016.
- dhdt.preprocessing.atmospheric_geometry.get_T_in_troposphere(h, T_0)¶
Using parameters of the first layer of the eight layer atmospheric model of ISO 2533:1970 to estimate the temperature at different altitudes
- Parameters
h ({float, numpy.array}, unit=m) – altitude above sea-level
T_0 (float, unit=Kelvin) – temperature at the surface (reference)
- Returns
T_h – temperature in the troposphere
- Return type
{float, numpy.array}, unit=Kelvin
References
- Ya16
Yan et al, “Correction of atmospheric refraction geolocation error of high resolution optical satellite pushbroom images” Photogrammetric engineering & remote sensing, vol.82(6) pp.427-435, 2016.
- dhdt.preprocessing.atmospheric_geometry.get_T_sealevel(lat, method='noerdlinger')¶
The temperature at sea level can be related to its latitude [Al73] and a least squares fit is given in [No99] and [Ya16].
- Parameters
lat ({float, numpy.ndarray}, unit=degrees, range=-90...+90) – latitude of interest
- Returns
T_sl ({float, numpy.array}, unit=Kelvin) – estimate of temperature at sea level
method ({‘noerdlinger’, ‘yan’}) –
‘noerdlinger’, from [2] has no dependency between hemispheres
’yan’, from [3] is hemisphere dependent
References
- Al73
Allen, “Astrophysical quantities”, pp.121-124, 1973.
- No99
Noerdlinger, “Atmospheric refraction effects in Earth remote sensing” ISPRS journal of photogrammetry and remote sensing, vol.54, pp.360-373, 1999.
- Ya16
Yan et al, “Correction of atmospheric refraction geolocation error of high resolution optical satellite pushbroom images” Photogrammetric engineering & remote sensing, vol.82(6) pp.427-435, 2016.
- dhdt.preprocessing.atmospheric_geometry.get_density_air(T, P, fH, moist=True, CO2=450, M_w=0.018015, R=8.31451, saturation=False)¶
- Parameters
T (float, unit=Kelvin) – temperature at location
P (float, unit=Pascal) – the total pressure in Pascals
fH (float, range=0...1) – fractional humidity
moist (bool) – estimate moist (True) or dry air (False)
CO2 (float, unit=ppm of CO₂) – parts per million of CO₂ in the atmosphere
M_w (float, unit=kg mol-1) – molar mass of water vapor
R (float, unit=J mol-1 K-1) – universal gas constant
- Returns
ρ – density of moist air
- Return type
float, unit=kg m3
Notes
The following parameters are used in this function:
M_a : unit=kg mol-1, molar mass of dry air
x_w : molecular fraction of water vapor
Z : compressibility of moist air
CO2 : unit=ppm, carbon dioxide
ppm : parts per million volume
References
- Ci96
Ciddor, “Refractive index of air: new equations for the visible and near infrared”, Applied optics, vol.35(9) pp.1566-1573, 1996.
- Da92
Davis, “Equation for the determination of the density of moist air (1981/91)” Metrologica, vol.29 pp.67-70, 1992.
- dhdt.preprocessing.atmospheric_geometry.get_density_fraction(φ, h_0, alpha=0.0065, R=8314.36, M_t=28.825)¶
- Parameters
ϕ (unit=degrees) – latitude of point of interest
h_0 (float, unit=meters) – initial altitude above geoid
alpha (float, unit=K km-1, default=0.0065) – lapse rate in the Troposphere
R (float, unit=J kmol-1 K-1, default=8314.36) – universal gas constant
M_t (float, unit=kg kmol-1, default=28.825) – mean moleculear mass in the Troposphere
- Returns
ρ_frac – fraction of density in relation to density as sea-level
- Return type
float
Notes
The following parameters are used
T_sl : unit=Kelvin, temperature at sea level T_t : unit=Kelvin, temperature at the Tropopause h_t : unit=meters, altitude of the Tropopause g : unit=m s-2, acceleration due to gravity
References
- No99
Noerdlinger, “Atmospheric refraction effects in Earth remote sensing” ISPRS journal of photogrammetry and remote sensing, vol.54, pp.360-373, 1999.
- dhdt.preprocessing.atmospheric_geometry.get_height_tropopause(φ)¶
The tropopause can be related to its latitude [Al73] and a least squares fit is given in [No99].
- Parameters
ϕ ({float, numpy.ndarray}, unit=degrees, range=-90...+90) – latitude of interest
- Returns
h_t – estimate of tropopause, that is the layer between the troposphere and the stratosphere
- Return type
{float, numpy.ndarray}, unit=meters
See also
get_mean_latNotes
The Strato- and Troposphere have different refractive behaviour. A simplified view of this can be like:
z_s Mesosphere / ┌-----------------/-┐ ↑ 80 km | * T_s |* | / | | | * |* | Stratosphere / | | | * |* | / | | | * |* ├--Tropopause-/-----┤ | ↑ 11 km |_*_________T_t |*________ | | | | | | ** | * | Tropo | | | | | ** | ** | sphere | z_h | | | | ** T_h | **** └------------#------┘ ↓ ↓ └-----------T_0 └--------- Earth surface Temperature Pressure
References
- Al73
Allen, “Astrophysical quantities”, pp.121-124, 1973.
- No99
Noerdlinger, “Atmospheric refraction effects in Earth remote sensing” ISPRS journal of photogrammetry and remote sensing, vol.54, pp.360-373, 1999.
- dhdt.preprocessing.atmospheric_geometry.get_sat_vapor_press(T, method='IAPWS')¶
uses either [1] or the international standard of IAPWS [2].
- Parameters
T ({float, numpy.array}, unit=Kelvin) – temperature
- Returns
svp – saturation vapor pressure
- Return type
{float, numpy.array}
References
- Gi82
Giacomo, “Equation for the determination of the density of moist air” Metrologica, vol.18, pp.33-40, 1982.
- WP02
Wagner and Pruß, “The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use.” Journal of physical and chemical reference data, vol.31(2), pp.387-535, 2002.
- dhdt.preprocessing.atmospheric_geometry.get_sat_vapor_press_ice(T)¶
following the international standard of IAPWS [2].
- Parameters
T ({float, numpy.array}, unit=Kelvin) – temperature
- Returns
svp – saturation vapor pressure
- Return type
{float, numpy.array}
References
- WP02
Wagner and Pruß, “The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use.” Journal of physical and chemical reference data, vol.31(2), pp.387-535, 2002.
- dhdt.preprocessing.atmospheric_geometry.get_simple_density_fraction(φ)¶
The density fraction is related to latitude through a least squares fit as is given in [No99].
- Parameters
ϕ ({float, numpy.array}, unit=degrees, range=-90...+90) – latitude of interest
- Returns
densFrac – estimate of the density fraction
- Return type
{float, numpy.array}
See also
get_mean_latReferences
- No99
Noerdlinger, “Atmospheric refraction effects in Earth remote sensing” ISPRS journal of photogrammetry and remote sensing, vol.54, pp.360-373, 1999.
- dhdt.preprocessing.atmospheric_geometry.get_water_vapor_enhancement(t, p)¶
estimate the enhancement factor for the humidity calculation
- Parameters
t ({float, numpy.array}, unit=Celsius) – temperature
p ({float, numpy.array}, unit=Pascal) – atmospheric pressure
- Returns
f – enhancement factor
- Return type
{float, numpy.array}
References
- Gi82
Giacomo, “Equation for the determination of the density of moist air” Metrologica, vol.18, pp.33-40, 1982.
- dhdt.preprocessing.atmospheric_geometry.get_water_vapor_press(T)¶
- Parameters
t ({float, numpy.array}, unit=Kelvin) – temperature
- Returns
P_w – water vapor pressure
- Return type
{float, numpy.array}, unit=Pascal
See also
References
- Ya16
Yan et al, “Correction of atmospheric refraction geolocation error of high resolution optical satellite pushbroom images” Photogrammetric engineering & remote sensing, vol.82(6) pp.427-435, 2016.
- dhdt.preprocessing.atmospheric_geometry.get_water_vapor_press_ice(T)¶
- Parameters
T ({float, numpy.array}, unit=Kelvin) – temperature
- Returns
P_w – water vapor pressure
- Return type
{float, numpy.array}, unit=Pascal
See also
References
- MM93
Marti & Mauersberger, “A survey and new measurements of ice vapor pressure at temperatures between 170 and 250K” Geophysical research letters, vol.20(5) pp.363-366, 1993.
- dhdt.preprocessing.atmospheric_geometry.refraction_angle_analytical(zn_0, n_0)¶
- Parameters
zn_0 ({float, numpy.array}, unit=degrees) – observation angle, when atmosphere is not taken into account
n_0 ({float, numpy.array}) – refractive index at the surface
- Returns
zn – corrected observation angle, taking refraction into account
- Return type
{float, numpy.array}, unit=degrees
References
- No99
Noerdlinger, “Atmospheric refraction effects in Earth remote sensing” ISPRS journal of photogrammetry and remote sensing, vol.54, pp.360-373, 1999.
- dhdt.preprocessing.atmospheric_geometry.refraction_spherical_symmetric(zn_0, lat=None, h_0=None)¶
- Parameters
zn_0 – observation angle, when atmosphere is not taken into account
float – observation angle, when atmosphere is not taken into account
lat (float, unit=degrees, range=-90...+90) – latitude of point of interest
h_0 (float, unit=meters) – initial altitude above geoid
- Returns
zn_hat – analytical refraction angle
- Return type
float
See also
get_mean_latNotes
The angles related to the sun are as follows:
surface normal * sun ↑ ↑ / | | / |-- zenith angle | / | / | /| |/ |/ | elevation angle +---- +------
References
- No99
Noerdlinger, “Atmospheric refraction effects in Earth remote sensing” ISPRS journal of photogrammetry and remote sensing, vol.54, pp.360-373, 1999.
- BJ94
Birch and Jones, “Correction to the updated Edlen equation for the refractive index of air”, Metrologica, vol.31(4) pp.315-316, 1994.
- dhdt.preprocessing.atmospheric_geometry.refractive_index_CO2(n_s, CO2)¶
correct standard refractivity for diffferent carbon dioxide concentrations, see [2] for more information.
- Parameters
n_s ({float, numpy.array}) – refractivity of standard air
CO2 ({float, numpy.array}, unit=ppm) – concentration of carbon dioxide (CO₂) in the atmosphere
- Returns
n_x – refractivity of standard air at a given CO₂ concentration
- Return type
{float, numpy.array}
Notes
The following acronyms are used:
CO2 : carbon dioxide
ppm : parts per million
References
- BD93
Birch & Downs, “An updated Edlen equation for the refractive index of air” Metrologica, vol.30(155) pp.155-162, 1993.
- Ed66
Edlén, “The refractive index of air”, Metrologica, vol.2(2) pp.71-80, 1966.
- dhdt.preprocessing.atmospheric_geometry.refractive_index_broadband(df, T_0, P_0, fH_0, CO2=450.0, T_a=15.0, P_a=101325.0, T_w=20.0, P_w=1333.0, LorentzLorenz=False)¶
Estimate the refractive index in for wavelengths in the range of 300…1700 nm, see also [1].
- Parameters
df ({dataframe, float}, unit=µm) – central wavelengths of the spectral bands
CO2 (float, unit=ppm of CO₂, default=450.) – parts per million of CO₂
T_0 (float, unit=Celsius) – temperature at location
P_0 (float, unit=Pascal) – the total pressure in Pascal
fH_0 (float, range=0...1) – fractional humidity
T_a (float, unit=Celcius) – temperature of standard dry air, see Notes
P_a (float, unit=Pascal) – pressure of standard dry air, see Notes
T_w (float, unit=Celsius) – temperature of standard moist air, see [1]
P_w (float, unit=Pascal) – pressure of standard moist air, see [1]
LorentzLorenz (bool) –
Two ways are possible to calculate the proportional refraction index, when
True : the Lorentz-Lorenz relationship is used.
False : component wise addition is used
- Returns
n – index of refraction
- Return type
{float, numpy.array}
Notes
Standard air [2] is defined to be dry air at a temperature of 15 degrees Celsius and a total pressure of 1013.25 mb, having the following composition by molar percentage:
78.09% nitrogen (N₂)
20.95% oxygen (O₂)
0.93% argon (Ar)
0.03% carbon dioxide (C0₂)
See also
ciddor_eq5calculation used for the Lorentz-Lorenz relationship
ciddor_eq6calculation used for component wise addition
References
- Ci96
Ciddor, “Refractive index of air: new equations for the visible and near infrared”, Applied optics, vol.35(9) pp.1566-1573, 1996.
- Ow67
Owens, “Optical refractive index of air: Dependence on pressure, temperature and composition”, Applied optics, vol.6(1) pp.51-59, 1967.
- dhdt.preprocessing.atmospheric_geometry.refractive_index_broadband_vapour(σ)¶
- Parameters
σ ({float, numpy.array}, unit=µm-1) – vacuum wavenumber, that is the reciprocal of the vacuum wavelength
- Returns
n_ws – refractive index for standard moist air at sea level
- Return type
{float, numpy.array}
References
- Ow67
Owens, “Optical refractive index of air: Dependence on pressure, temperature and composition”, Applied optics, vol.6(1) pp.51-59, 1967.
- Ci96
Ciddor, “Refractive index of air: new equations for the visible and near infrared”, Applied optics, vol.35(9) pp.1566-1573, 1996.
- dhdt.preprocessing.atmospheric_geometry.refractive_index_simple(p, RH, t)¶
simple shop-floor formula fro refraction, see appendix B in [1].
- Parameters
p ({numpy.array, float}, unit=kPa) – pressure
RH ({numpy.array, float}, unit=percent) – relative humidity
t ({numpy.array, float}, unit=Celsius) – temperature
- Returns
n – refractive index
- Return type
{numpy.array, float}
References
- dhdt.preprocessing.atmospheric_geometry.refractive_index_visible(df, T, P, p_w=None, CO2=None)¶
calculate refractive index, based on measurement around 633nm, see [1].
It applies to ambient atmospheric conditions over the range of wavelengths 350 nm to 650 nm.
- Parameters
df ({dataframe, float}, unit=µm) – central wavelengths of the spectral bands
T ({float, numpy.array}, unit=Celcius) – temperature at location
P ({float, numpy.array}, unit=Pascal) – atmospheric pressure at location
CO2 ({float, numpy.array}, unit=ppm) – concentration of carbon dioxide (CO₂) in the atmosphere
- Returns
n_0 – index of refraction
- Return type
numpy.array
See also
References
- BJ94
Birch and Jones, “Correction to the updated Edlén equation for the refractive index of air”, Metrologica, vol.31(4) pp.315-316, 1994.
- Ed66
Edlén, “The refractive index of air”, Metrologica, vol.2(2) pp.71-80, 1966.
- BJ94
Birch and Jones, “An updated Edlén equation for the refractive index of air”, Metrologica, vol.30 pp.155-162, 1993.
- dhdt.preprocessing.atmospheric_geometry.update_list_with_corr_zenith_pd(dh, file_dir, file_name='conn.txt')¶
- Parameters
dh (pandas.DataFrame) – array with photohypsometric information
file_dir (string) – directory where the file should be situated
file_name (string, default='conn.txt') – name of the file to export the data columns to
Notes
The nomenclature is as follows:
* sun \ <-- sun trace |\ <-- caster | | \ ^ Z | \ | ----┴----+ <-- casted └-> {X,Y}The angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / |-- zenith angle | / | / | /| |/ |/ | elevation angle └---- └------ {X,Y}The azimuth angle declared in the following coordinate frame:
^ North & Y | - <--|--> + | +----> East & XThe text file has at least the following columns:
* conn.txt ├ caster_X <- map location of the start of the shadow trace ├ caster_Y ├ casted_X ├ casted_Y <- map location of the end of the shadow trace ├ azimuth <- argument of the illumination direction ├ zenith <- off-normal angle of the sun without atmosphere └ zenith_refrac <- off-normal angle of the sun with atmosphere
See also
get_refraction_angle
- dhdt.preprocessing.atmospheric_geometry.water_vapor_frac(p, T, fH)¶
- Parameters
p (float, unit=Pascal) – total pressure
T (float, unit=Kelvin) – temperature
fH (float, range=0...1) – fractional humidity
- Returns
x_w – molar fraction of water vapor in moist air
- Return type
float
Notes
The following parameters are used in this function:
T : unit=Kelvin, temperature svp : unit=Pascal, saturation vapor pressure t : unit=Celcius, temperature
References
- Ci96
Ciddor, “Refractive index of air: new equations for the visible and near infrared”, Applied optics, vol.35(9) pp.1566-1573, 1996.
- MM93
Marti & Mauersberger, “A survey and new measurements of ice vapor pressure at temperatures between 170 and 250K” Geophysical research letters, vol.20(5) pp.363-366, 1993.
- Gi82
Giacomo, “Equation for the determination of the density of moist air” Metrologica, vol.18, pp.33-40, 1982.
processing¶
organization¶
- dhdt.processing.matching_tools_organization.estimate_precision(C, di, dj, method='gaussian')¶
given a similarity surface, estimate its matching dispersion. See for a full describtion [Al21]
- Parameters
C (numpy.ndarray, size=(m,n)) – cross correlation function
di (float) – locational estimate of the cross-correlation peak, this can also be negative, see Notes below
dj (float) – locational estimate of the cross-correlation peak, this can also be negative, see Notes below
method ({'gaussian', 'hessian'}) –
- Returns
cov_ii, cov_jj (float) – co-variance estimate
cov_ij (float) – off-diagonal co-variance estimate
Notes
It is important to know what type of coordinate systems exist, hence:
coordinate | +--------> collumns system 'ij'| | | | | j | image frame --------+--------> | | | | v | i rows vSee also
estimate_subpixel,match_translation_of_two_subsets,hessian_spread,gauss_spreadReferences
- Al21
Altena et al. “Correlation dispersion as a measure to better estimate uncertainty of remotely sensed glacier displacements” The cryosphere, vol.16(6) pp.2285-2300, 2021.
- dhdt.processing.matching_tools_organization.estimate_subpixel(QC, subpix, m0=array([[0., 0.]]), **kwargs)¶
estimate the sub-pixel translational displacement, present in a cross-correlation function or spectra.
- Parameters
QC (numpy.ndarray, size=(m,n), dtype={complex,float}, ndim=2) – cross-correlation spectra or function
subpix (string) – methodology to refine the to sub-pixel localization
m0 (numpy.ndarray, size=(1,2), dtype=integer) – row, collumn location of the peak of the cross-correlation function
- Returns
ddi,ddj – sub-pixel location of the peak/phase-plane
- Return type
float
- dhdt.processing.matching_tools_organization.estimate_translation_of_two_subsets(I1, I2, M1, M2, correlator='lucas_kan', **kwargs)¶
- Parameters
I1 (numpy.ndarray, size=(m,n), dtype={float,integer}) – grids with intensities
I2 (numpy.ndarray, size=(m,n), dtype={float,integer}) – grids with intensities
M1 (numpy.ndarray, size=(m,n), dtype=bool) – labelled array corresponding to I1,I2. Hence True denotes exclusion.
M2 (numpy.ndarray, size=(m,n), dtype=bool) – labelled array corresponding to I1,I2. Hence True denotes exclusion.
correlator (string, default='lucas_kanade') – methodology to use to correlate I1 and I2
- Returns
di, dj (float, unit=pixels) – relative displacement (translation) between I1 and I2
score (float) – (dis)-similarity score
- dhdt.processing.matching_tools_organization.list_differential_correlators()¶
” list the abbreviations of the different implemented differential correlators.
The following are implemented:
lucas_kan
lucas_aff
hough_opt_flw
- Returns
correlator_list
- Return type
list of strings
See also
- dhdt.processing.matching_tools_organization.list_frequency_correlators(unique=False)¶
list the abbreviations of the different implemented correlators
The following correlators are implemented:
cosi_corr : cosicorr
phas_only : phase only correlation
symm_phas : symmetric phase correlation
ampl_comp : amplitude compensation phase correlation
orie_corr : orientation correlation
grad_corr : gradient correlation
grad_norm : normalize gradient correlation
bina_phas : binary phase correlation
wind_corr : windrose correlation
gaus_phas : gaussian transformed phase correlation
cros_corr : cross correlation
robu_corr : robust phase correlation
proj_phas : projected phase correlation
phas_corr : phase correlation
- Returns
correlator_list
- Return type
list of strings
- dhdt.processing.matching_tools_organization.list_peak_estimators()¶
list the abbreviations of the different implemented for the peak fitting of the similarity function.
The following methods are implemented:
‘gauss_1’ : 1D Gaussian fitting
‘parab_1’ : 1D parabolic fitting
‘moment’ : 2D moment of the peak
‘mass’ : 2D center of mass fitting
‘centroid’ : 2D estimate the centroid
‘blais’ : 1D estimate through forth order filter
‘ren’ : 1D parabolic fitting
‘birch’ : 1D weighted sum
‘eqang’ : 1D equiangular line fitting
‘trian’ : 1D triangular fitting
‘esinc’ : 1D exponential esinc function
‘gauss_2’ : 2D Gaussian fitting
‘parab_2’ : 2D parabolic fitting
‘optical_flow’ : optical flow refinement
- Returns
subpix_list
- Return type
list of strings
- dhdt.processing.matching_tools_organization.list_phase_estimators()¶
list the abbreviations of the different implemented phase plane estimation procedures
The following methods are implemented:
‘tpss’ : two point step size
‘svd’ : single value decomposition
‘radon’ : radon transform
‘hough’ : hough transform
‘ransac’ : random sampling and consensus
‘wpca’ : weighted principle component analysis
‘pca’ : principle component analysis
‘lsq’ : least squares estimation
‘diff’ : phase difference
- Returns
subpix_list
- Return type
list of strings
- dhdt.processing.matching_tools_organization.list_spatial_correlators(unique=False)¶
list the abbreviations of the different implemented correlators.
The following methods are implemented:
norm_corr : normalized cross correlation
sq_diff : sum of squared differences
sad_diff : sum of absolute differences
max_like : maximum likelihood
wght_corr : weighted normalized cross correlation
- Returns
correlator_list
- Return type
list of strings
- dhdt.processing.matching_tools_organization.match_translation_of_two_subsets(I1_sub, I2_sub, correlator, subpix, M1_sub=array([], dtype=float64), M2_sub=array([], dtype=float64))¶
- Parameters
I1_sub (numpy.ndarray, size={(m,n),(k,l)}, dtype={float,int}, ndim={2,3}) – grid with intensities, a.k.a. template to locate
I2_sub (numpy.ndarray, size=(m,n), dtype={float,int}, ndim={2,3}) – grid with intensities, a.k.a. search space
correlator (string) – methodology to use to correlate I1 and I2
subpix (string) – methodology to refine the to sub-pixel localization
M1_sub (numpy.ndarray, size=(m,n), dtype=bool) – labels corresponding to I1_sub,I2_sub. Hence True denotes exclusion.
M2_sub (numpy.ndarray, size=(m,n), dtype=bool) – labels corresponding to I1_sub,I2_sub. Hence True denotes exclusion.
- Returns
{Q,C} – cross-correlation spectra or function
- Return type
numpy.ndarray, size={(m,n),(k,l)}, dtype={complex,float}, ndim=2
frequency filters¶
- dhdt.processing.matching_tools_frequency_filters.adaptive_masking(S, m=0.9)¶
mark significant intensities in spectrum, following [Le07].
- Parameters
S (numpy.array, size=(m,n), dtype=complex) – array with spectrum, i.e.: S = np.fft.fft2(I)
m (float, default=.9) – cut-off intensity in respect to maximum
- Returns
M – frequency mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
tpssReferences
- Le07
Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.
- dhdt.processing.matching_tools_frequency_filters.blackman_window(Z)¶
create two-dimensional Blackman filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
raised_cosine,cosine_bell,high_pass_circle,hamming_window,hanning_window
- dhdt.processing.matching_tools_frequency_filters.construct_phase_plane(Z, di, dj, indexing='ij')¶
given a displacement, create what its phase plane in Fourier space
- Parameters
Z (np.array, size=(m,n)) – image domain
di (float) – displacement along the vertical axis
dj (float) – displacement along the horizantal axis
indexing ({‘xy’, ‘ij’}) –
“xy” : using map coordinates
”ij” : using local image coordinates
- Returns
Q – array with phase angles
- Return type
np.array, size=(m,n), complex
Notes
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | i | x --------+--------> --------+--------> | | | | image | j map | based v based |
- dhdt.processing.matching_tools_frequency_filters.construct_phase_values(IJ, di, dj, indexing='ij', system='radians')¶
given a displacement, create what its phase plane in Fourier space
- Parameters
IJ (np.array, size=(_,2), dtype=float) – locations of phase values
di (float) – displacment along the vertical axis
dj (float) – displacment along the horizantal axis
indexing ({‘radians’ (default), ‘unit’, 'normalized'}) –
“xy” : using map coordinates
”ij” : using local image coordinates
indexing –
the extent of the cross-spectrum can span different ranges
”radians” : -pi..+pi
”unit” : -1…+1
”normalized” : -0.5…+0.5
- Returns
Q – array with phase angles
- Return type
np.array, size=(_,1), complex
- dhdt.processing.matching_tools_frequency_filters.cosine_bell(Z)¶
cosine bell filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=float
See also
- dhdt.processing.matching_tools_frequency_filters.cross_shading_filter(Q, az_1, az_2)¶
- Parameters
Q (numpy.array, size=(m,n)) – cross-power spectrum
az_1 (float, unit=degrees, range=-180...+180) – illumination direction of template
az_2 (float, unit=degrees, range=-180...+180) – illumination direction of search space
- Returns
W – weigthing grid
- Return type
numpy.array, size=(m,n), dtype=float
References
- Wa15
Wan, “Phase correlation-based illumination-insensitive image matching for terrain-related applications” PhD dissertation at Imperical College London, 2015.
- dhdt.processing.matching_tools_frequency_filters.cross_spectrum_to_coordinate_list(data, W=array([], dtype=float64))¶
if data is given in array for, then transform it to a coordinate list
- Parameters
data (np.array, size=(m,n), dtype=complex) – cross-spectrum
W (np.array, size=(m,n), dtype=boolean) – weigthing matrix of the cross-spectrum
- Returns
data_list – coordinate list with angles, in normalized ranges, i.e: -1 … +1
- Return type
np.array, size=(m*n,3), dtype=float
- dhdt.processing.matching_tools_frequency_filters.gaussian_mask(S)¶
mask significant intensities in spectrum, following [Ec08].
- Parameters
S (numpy.array, size=(m,n), dtype=complex) – array with spectrum, i.e.: S = np.fft.fft2(I)
- Returns
M – frequency mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
tpssReferences
- Ec08
Eckstein et al. “Phase correlation processing for DPIV measurements”, Experiments in fluids, vol.45 pp.485-500, 2008.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1) >>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2) >>> Q = spec1 * np.conjugate(spec2) # Fourier based image matching >>> Qn = normalize_spectrum(Q)
>>> W = gaussian_mask(Q) >>> C = np.fft.ifft2(W*Q)
- dhdt.processing.matching_tools_frequency_filters.gradient_fourier(Z)¶
spectral derivative estimation
- Parameters
Z (numpy.ndarray, size=(m,n), dtype=float) – intensity array
- Returns
dZdx_1, dZdx_2 – first order derivative along both axis
- Return type
numpy.ndarray, size=(m,n), dtype=float
- dhdt.processing.matching_tools_frequency_filters.hamming_window(Z)¶
create two-dimensional Hamming filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
- dhdt.processing.matching_tools_frequency_filters.hanning_window(Z)¶
create two-dimensional Hanning filter, also known as Cosine Bell
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
- dhdt.processing.matching_tools_frequency_filters.high_pass_circle(Z, r=0.5)¶
create hard high-pass filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
r (float, default=0.5) – radius of the circle, r=.5 is same as its width
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
- dhdt.processing.matching_tools_frequency_filters.kaiser_window(Z, beta=14.0)¶
create two dimensional Kaiser filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
beta (float) – 0.0 - rectangular window 5.0 - similar to Hamming window 6.0 - similar to Hanning window 8.6 - similar to Blackman window
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
raised_cosine,cosine_bell,high_pass_circle,hamming_window,hanning_window
- dhdt.processing.matching_tools_frequency_filters.local_coherence(Q, ds=1)¶
estimate the local coherence of a spectrum
- Parameters
Q (numpy.ndarray, size=(m,n), dtype=complex) – array with cross-spectrum, with centered coordinate frame
ds (integer, default=1) – kernel radius to describe the neighborhood
- Returns
M – vector coherence from no to ideal, i.e.: 0…1
- Return type
numpy.ndarray, size=(m,n), dtype=float
See also
Examples
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> # create cross-spectrum with random displacement >>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1) >>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2) >>> Q = spec1 * np.conjugate(spec2) >>> Q = normalize_spectrum(Q) >>> Q = np.fft.fftshift(Q) # transform to centered grid
>>> C = local_coherence(Q)
>>> plt.figure(), plt.imshow(C, cmap='OrRd'), plt.colorbar(), plt.show() >>> plt.figure(), plt.imshow(np.angle(Q), cmap='twilight'), >>> plt.colorbar(), plt.show()
- dhdt.processing.matching_tools_frequency_filters.low_pass_bell(Z, r=0.5)¶
create low-pass two-dimensional filter with a bell shape, see also [Ta03].
- Parameters
Z (numpy.ndarray, size=(m,n)) – array with intensities
r (float, default=0.5) – radius of the mother rectangle, r=.5 is same as its width
- Returns
W – weighting mask
- Return type
numpy.ndarray, size=(m,n), dtype=bool
See also
References
- Ta03
Takita et al. “High-accuracy subpixel image registration based on phase-only correlation” IEICE transactions on fundamentals of electronics, communications and computer sciences, vol.86(8) pp.1925-1934, 2003.
- dhdt.processing.matching_tools_frequency_filters.low_pass_circle(Z, r=0.5)¶
create hard two-dimensional low-pass filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
r (float, default=0.5) – radius of the circle, r=.5 is same as its width
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
- dhdt.processing.matching_tools_frequency_filters.low_pass_ellipse(Z, r1=0.5, r2=0.5)¶
create hard low-pass filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
r1 (float, default=0.5) – radius of the ellipse along the first axis, r=.5 is same as its width
r2 (float, default=0.5) – radius of the ellipse along the second axis
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
raised_cosine,cosine_bell,high_pass_circle,low_pass_circle
- dhdt.processing.matching_tools_frequency_filters.low_pass_pyramid(Z, r=0.5)¶
create low-pass two-dimensional filter with pyramid shape, see also [Ta03].
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
r (float, default=0.5) – radius of the mother rectangle, r=.5 is same as its width
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
References
- Ta03
Takita et al. “High-accuracy subpixel image registration based on phase-only correlation” IEICE transactions on fundamentals of electronics, communications and computer sciences, vol.86(8) pp.1925-1934, 2003.
- dhdt.processing.matching_tools_frequency_filters.low_pass_rectancle(Z, r=0.5)¶
create hard two dimensional low-pass filter
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
r (float, default=0.5) – radius of the rectangle, r=.5 is same as its width
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
References
- Ta03
Takita et al. “High-accuracy subpixel image registration based on phase-only correlation” IEICE transactions on fundamentals of electronics, communications and computer sciences, vol.86(8) pp.1925-1934, 2003.
- dhdt.processing.matching_tools_frequency_filters.make_fourier_grid(Q, indexing='ij', system='radians', shift=True, axis='center')¶
The four quadrants of the coordinate system of the discrete Fourier transform are flipped. This function gives its coordinate system as it would be in a map (xy) or pixel based (ij) system.
- Parameters
Q (numpy.ndarray, size=(m,n), dtype=complex) – Fourier based (cross-)spectrum.
indexing ({‘xy’, ‘ij’}) –
“xy” : using map coordinates
”ij” : using local image coordinates
system ({‘radians’, ‘unit’, 'normalized'}) –
- the extent of the cross-spectrum can span from
”radians” : -pi..+pi (default)
”unit” : -1…+1
”normalized” : -0.5…+0.5
”pixel” : -m/2…+m/2
shift (boolean, default=True) – the coordinate frame of a Fourier transform is flipped
- Returns
F_1 (numpy.ndarray, size=(m,n), dtype=integer) – first coordinate index of the Fourier spectrum in a map system.
F_2 (numpy.ndarray, size=(m,n), dtype=integer) – second coordinate index of the Fourier spectrum in a map system.
Notes
metric system: Fourier-based flip y ┼------><------┼ ^ | | | | | | v v <------┼-------> x | ^ ^ | | | v ┼------><------┼It is important to know what type of coordinate systems exist, hence:
coordinate | coordinate ^ y system 'ij'| system 'xy' | | | | j | x --------┼--------> --------┼--------> | | | | | i | v |
- dhdt.processing.matching_tools_frequency_filters.make_template_float(Z)¶
templates for frequency matching should be float and ideally have limited border issues.
- Parameters
Z (numpy.ndarray, size={(m,n), (m,n,b)}) – grid with intensity values
- Returns
Z – grid with intensity values
- Return type
numpy.ndarray, size={(m,n), (m,n,b)}, dtype=float, range=0…1
- dhdt.processing.matching_tools_frequency_filters.normalize_power_spectrum(Q)¶
transform spectrum to complex vectors with unit length
- Parameters
Q (numpy.ndarray, size=(m,n), dtype=complex) – cross-spectrum
- Returns
Qn – normalized cross-spectrum, that is elements with unit length
- Return type
numpy.ndarray, size=(m,n), dtype=complex
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1) >>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2) >>> Q = spec1 * np.conjugate(spec2) # fourier based image matching >>> Qn = normalize_spectrum(Q)
- dhdt.processing.matching_tools_frequency_filters.perdecomp(img)¶
calculate the periodic and smooth components of an image, based upon [Mo11].
- Parameters
img (numpy.ndarray, size=(m,n)) – array with intensities
- Returns
per (numpy.ndarray, size=(m,n)) – periodic component
cor (numpy.ndarray, size=(m,n)) – smooth component
References
- Mo11
Moisan, L. “Periodic plus smooth image decomposition”, Journal of mathematical imaging and vision vol. 39.2 pp. 161-179, 2011.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,_,_,_,_ = create_sample_image_pair(d=2**7, max_range=1) >>> per,cor = perdecomp(im1)
>>> spec1 = np.fft.fft2(per)
- dhdt.processing.matching_tools_frequency_filters.raised_cosine(Z, beta=0.35)¶
raised cosine filter, based on [St01] and used by [Le07].
- Parameters
Z (numpy.array, size=(m,n)) – array with intensities
beta (float, default=0.35) – roll-off factor
- Returns
W – weighting mask
- Return type
numpy.array, size=(m,n), dtype=float
See also
tpssReferences
- St01
Stone et al. “A fast direct Fourier-based algorithm for subpixel registration of images.” IEEE Transactions on geoscience and remote sensing. vol. 39(10) pp. 2235-2243, 2001.
- Le07
Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1) >>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2)
>>> rc1 = raised_cosine(spec1, beta=0.35) >>> rc2 = raised_cosine(spec2, beta=0.50)
>>> Q = (rc1*spec1) * np.conjugate((rc2*spec2)) # Fourier based image matching >>> Qn = normalize_spectrum(Q)
- dhdt.processing.matching_tools_frequency_filters.thresh_masking(S, m=0.0001, s=10)¶
Mask significant intensities in spectrum, following [St01] and [Le07].
- Parameters
S (numpy.array, size=(m,n), dtype=complex) – array with spectrum, i.e.: S = np.fft.fft2(I)
m (float, default=1e-3) – cut-off intensity in respect to maximum
s (integer, default=10) – kernel size of the median filter
- Returns
M – frequency mask
- Return type
numpy.array, size=(m,n), dtype=bool
See also
tpssReferences
- St01
Stone et al. “A fast direct Fourier-based algorithm for subpixel registration of images.” IEEE Transactions on geoscience and remote sensing vol. 39(10) pp. 2235-2243, 2001.
- Le07
Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.
frequency correlators¶
- dhdt.processing.matching_tools_frequency_correlators.amplitude_comp_corr(I1, I2, F_0=0.04)¶
match two imagery through amplitude compensated phase correlation
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
F_0 (float, default=4e-2) – cut-off intensity in respect to maximum
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
References
- Mu88
Mu et al. “Amplitude-compensated matched filtering”, Applied optics, vol. 27(16) pp. 3461-3463, 1988.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = amplitude_comp_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.binary_orientation_corr(I1, I2)¶
match two imagery through binary phase only correlation, see [KJ91].
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.orientation_corr,dhdt.processing.matching_tools_frequency_correlators.phase_only_corrNotes
The matching equations are as follows:
\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]\[\mathbf{W} = sign(\mathfrak{Re} [ \mathbf{S}_2 ] )\]\[\mathbf{Q}_{12} = \mathbf{S}_1 [\mathbf{W}\mathbf{S}_2]^{\star}\]where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation
References
- KJ91
Kumar & Juday, “Design of phase-only, binary phase-only, and complex ternary matched filters with increased signal-to-noise ratios for colored noise”, Optics letters, vol.16(13) pp.1025-1027, 1991.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = binary_orientation_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.cosine_corr(I1, I2)¶
match two imagery through discrete cosine transformation, see [Li04].
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}, dtype=float) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}, dtype=float) – array with intensities
- Returns
C – displacement surface
- Return type
numpy.array, size=(m,n), dtype=real
See also
dhdt.processing.matching_tools_frequency_correlators.create_complex_DCT,dhdt.processing.matching_tools_frequency_correlators.sign_only_corrReferences
- Li04
Li, et al. “DCT-based phase correlation motion estimation”, IEEE international conference on image processing, vol. 1, 2004.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
- dhdt.processing.matching_tools_frequency_correlators.cross_corr(I1, I2)¶
match two imagery through cross correlation in FFT
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.phase_corr,dhdt.processing.matching_tools_frequency_correlators.upsampled_cross_corrNotes
The matching equations are as follows:
\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]\[\mathbf{Q}_{12} = \mathbf{S}_1 {\mathbf{S}_2}^{\star}\]where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation.
Schematically this looks like:
I1 ┌----┐ S1 ---┤ F ├-----┐ └----┘ | Q12 ┌------┐ C12 ×-----┤ F^-1 ├--- I2 ┌----┐ S2 | └------┘ ---┤ F ├-----┘ └----┘References
- HK12
Heid & Kääb. “Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery”, Remote sensing of environment, vol.118 pp.339-355, 2012.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = cross_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.gaussian_transformed_phase_corr(I1, I2)¶
match two imagery through Gaussian transformed phase correlation, see [Ec08].
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
Notes
The matching equations are as follows:
\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]\[\mathbf{Q}_{12} = \mathbf{S}_1 {\mathbf{S}_2}^{\star}\]\[\mathbf{Q}_{12} / \| \mathbf{Q}_{12}\|\]\[\mathbf{Q}_{12} = \mathbf{W} \mathbf{Q}_{12}\]where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation.
Schematically this looks like:
I1 ┌----┐ S1 ┌----┐ ---┤ F ├-----┐ | W | └----┘ | Q12 ┌-------┐ └-┬--┘ ┌------┐ C12 ×-----┤ 1/|x| ├---┴----┤ F^-1 ├---- I2 ┌----┐ S2 | └-------┘ └------┘ ---┤ F ├-----┘ └----┘References
- Ec08
Eckstein et al. “Phase correlation processing for DPIV measurements”, Experiments in fluids, vol.45 pp.485-500, 2008.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = gaussian_transformed_phase_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.gradient_corr(I1, I2)¶
match two imagery through gradient correlation
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.normalized_gradient_corr,dhdt.processing.matching_tools_frequency_correlators.phase_corr,dhdt.processing.matching_tools_frequency_correlators.orientation_corrNotes
The matching equations are as follows:
\[\mathbf{G}_1, \mathbf{G}_2 = \partial_{x}\mathbf{I}_1 + i \partial_{y}\mathbf{I}_1, \quad \partial_{x}\mathbf{I}_2 + i \partial_{y}\mathbf{I}_2\]\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{G}_1], \quad \mathcal{F}[\mathbf{G}_2]\]\[\mathbf{Q}_{12} = \mathbf{S}_1 \mathbf{S}_2^{\star}\]where \(\partial_{x}\) is the spatial derivate in horizontal direction, \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation.
Schematically this looks like:
∂x I1 I1 ┌----┬-----┬----┐ S1 ---┤ ∂ | | F ├------┐ └----┴-----┴----┘ | ∂y I1 | | Q12 ┌------┐ C12 ∂x I2 ×-----┤ F^-1 ├--- I2 ┌----┬-----┬----┐ S2 | └------┘ ---┤ ∂ | | F ├------┘ └----┴-----┴----┘ ∂y I2References
- AV03
Argyriou & Vlachos. “Estimation of sub-pixel motion using gradient cross-correlation”, Electronics letters, vol.39(13) pp.980–982, 2003.
- Tz10
Tzimiropoulos et al. “Subpixel registration with gradient correlation” IEEE transactions on image processing, vol.20(6) pp.1761–1767, 2010.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = gradient_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.masked_corr(I1, I2, M1=array([], dtype=float64), M2=array([], dtype=float64))¶
match two imagery through masked normalized cross-correlation in FFT, see [Pa11].
- Parameters
I1 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities
I2 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities
M1 (numpy.array, size=(m,n)) – array with mask, hence False’s are neglected
M2 (numpy.array, size=(m,n)) – array with mask, hence False’s are neglected
- Returns
NCC – correlation surface
- Return type
numpy.array, size=(m,n)
See also
dhdt.processing.matching_tools_spatial_correlators.normalized_cross_corr,dhdt.processing.matching_tools_spatial_correlators.weighted_normalized_cross_corrReferences
- Pa11
Padfield. “Masked object registration in the Fourier domain”, IEEE transactions on image processing, vol. 21(5) pp. 2706-2718, 2011.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> msk1,msk2 = np.ones_like(im1), np.ones_like(im2) >>> Q = masked_corr(im1, im2, msk1, msk2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.masked_cosine_corr(I1, I2, M1, M2)¶
- Parameters
I1 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities
I2 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities
M1 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected
M2 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected
- dhdt.processing.matching_tools_frequency_correlators.normalized_gradient_corr(I1, I2)¶
match two imagery through gradient correlation
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.gradient_corr,dhdt.processing.matching_tools_frequency_correlators.phase_corr,dhdt.processing.matching_tools_frequency_correlators.windrose_corrNotes
The matching equations are as follows:
\[\mathbf{G}_1, \mathbf{G}_2 = \partial_{x}\mathbf{I}_1 + i \partial_{y}\mathbf{I}_1, \quad \partial_{x}\mathbf{I}_2 + i \partial_{y}\mathbf{I}_2\]\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\dot{\mathbf{G}}_1], \quad \mathcal{F}[\dot{\mathbf{G}}_2]\]\[\mathbf{Q}_{12} = \mathbf{S}_1 \mathbf{S}_2^{\star}\]where \(\partial_{x}\) is the spatial derivate in horizontal direction, \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation and \(\dot{x}\) represents a statistical normalization through the sampled mean and standard deviation.
Schematically this looks like:
∂x I1 I1 ┌----┬-----┬----┐ S1 ┌------┐ ---┤ ∂ | | F ├----┤ norm ├--┐ └----┴-----┴----┘ └------┘ | ∂y I1 | | Q12 ┌------┐ C12 ∂x I2 ×-----┤ F^-1 ├--- I2 ┌----┬-----┬----┐ S2 ┌------┐ | └------┘ ---┤ ∂ | | F ├----┤ norm ├--┘ └----┴-----┴----┘ └------┘ ∂y I2References
- Tz10
Tzimiropoulos et al. “Subpixel registration with gradient correlation” IEEE transactions on image processing, vol.20(6) pp.1761–1767, 2010.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = normalized_gradient_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.orientation_corr(I1, I2)¶
match two imagery through orientation correlation, developed by [Fi02] while demonstrated its benefits over glaciated terrain by [HK12].
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.phase_corr,dhdt.processing.matching_tools_frequency_correlators.windrose_corr,dhdt.processing.matching_tools_spatial_correlators.cosine_similarityNotes
The matching equations are as follows:
\[\mathbf{G}_1, \mathbf{G}_2 = \partial_{x}\mathbf{I}_1 + i \partial_{y}\mathbf{I}_1,\quad \partial_{x}\mathbf{I}_2 + i \partial_{y}\mathbf{I}_2\]\[\hat{\mathbf{G}_1}, \hat{\mathbf{G}_2} = \mathbf{G}_1 / \| \mathbf{G}_1\|,\quad \mathbf{G}_2 / \| \mathbf{G}_2\|\]\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\hat{\mathbf{G}_1}],\quad \mathcal{F}[\hat{\mathbf{G}_2}]\]\[\mathbf{Q}_{12} = \mathbf{S}_1 \mathbf{S}_2^{\star}\]where \(\partial_{x}\) is the spatial derivate in horizontal direction, \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation and \(\hat{x}\) is a normalizing operator.
Schematically this looks like:
∂x I1┌-------┐ I1 ┌----┬--┤ 1/|x| ├--┬----┐ S1 →--┤ ∂ | ├-------┤ | F ├--→-┐ └----┴--┤ 1/|x| ├--┴----┘ | ∂y I1└-------┘ | | Q12 ┌------┐ C12 ×--→--┤ F^-1 ├--- ∂x I2┌-------┐ | └------┘ I2 ┌----┬--┤ 1/|x| ├--┬----┐ S2 | →--┤ ∂ | ├-------┤ | F ├--→-┘ └----┴--┤ 1/|x| ├--┴----┘ ∂y I2└-------┘References
- Fi02
Fitch et al. “Orientation correlation”, Proceeding of the Britisch machine vison conference, pp. 1–10, 2002.
- HK12
Heid & Kääb. “Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery”, Remote sensing of environment, vol.118 pp.339-355, 2012.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = orientation_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.phase_corr(I1, I2)¶
match two imagery through phase correlation
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.orientation_corr,dhdt.processing.matching_tools_frequency_correlators.cross_corrNotes
The matching equations are as follows:
\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]\[\hat{\mathbf{S}}_1, \hat{\mathbf{S}}_2 = \mathbf{S}_1 / \| \mathbf{S}_1\|, \quad \mathbf{S}_2 / \| \mathbf{S}_2\|\]\[\mathbf{Q}_{12} = \hat{\mathbf{S}}_1 {\hat{\mathbf{S}}_2}^{\star}\]where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation and \(\hat{x}\) is a normalizing operator.
Schematically this looks like:
I1 ┌----┐ S1 ┌-------┐ ---┤ F ├----┤ 1/|x| ├--┐ └----┘ └-------┘ | Q12 ┌------┐ C12 ×-----┤ F^-1 ├--- I2 ┌----┐ S2 ┌-------┐ | └------┘ ---┤ F ├----┤ 1/|x| ├--┘ └----┘ └-------┘References
- KH75
Kuglin & Hines. “The phase correlation image alignment method”, proceedings of the IEEE international conference on cybernetics and society, pp. 163-165, 1975.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = phase_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.phase_only_corr(I1, I2)¶
match two imagery through phase only correlation, based upon [HG84] and [KJ91].
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.phase_corr,dhdt.processing.matching_tools_frequency_correlators.symmetric_phase_corr,dhdt.processing.matching_tools_frequency_correlators.amplitude_comp_corrNotes
The matching equations are as follows:
\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]\[\mathbf{W} = 1 / \|\mathbf{S}_2 \|\]\[\mathbf{Q}_{12} = \mathbf{S}_1 [\mathbf{W}\mathbf{S}_2]^{\star}\]where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation
Schematically this looks like:
I2 ┌----┐ S2 ---┤ F ├--┬------------┐ └----┘ | ┌-------┐W | └-┤ 1/|x| ├--× └-------┘ | Q12 ┌------┐ C12 ×-----┤ F^-1 ├--- I1 ┌----┐ S1 | └------┘ ---┤ F ├---------------┘ └----┘If multiple bands are given, cross-spectral stacking is applied,see [AL22].
References
- HG84
Horner & Gianino, “Phase-only matched filtering”, Applied optics, vol. 23(6) pp.812–816, 1984.
- KJ91
Kumar & Juday, “Design of phase-only, binary phase-only, and complex ternary matched filters with increased signal-to-noise ratios for colored noise”, Optics letters, vol.16(13) pp.1025–1027, 1991.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair >>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = phase_only_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.projected_phase_corr(I1, I2, M1=array([], dtype=float64), M2=array([], dtype=float64))¶
match two imagery through separated phase correlation
- Parameters
I1 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities
I2 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities
M1 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected
M2 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected
- Returns
C – displacement surface
- Return type
numpy.array, size=(m,n), dtype=real
Notes
Schematically this looks like:
I1 ┌----┐ I1x ┌----┐ S1x -┬-┤ Σx ├------┤ F ├------┐ Q12x ┌------┐ C12x | └----┘ └----┘ ×------┤ F^-1 ├------┐ | ┌----┐ I2x ┌----┐ S2x | └------┘ | |┌┤ Σx ├------┤ F ├------┘ | ┌---┐ C12 ||└----┘ └----┘ ×--┤ √ ├---- ||┌----┐ I1y ┌----┐ S1y | └---┘ └┼┤ Σy ├------┤ F ├------┐ Q12y ┌------┐ C12y | |└----┘ └----┘ ×------┤ F^-1 ├------┘ I2 |┌----┐ I2y ┌----┐ S2y | └------┘ --┴┤ Σy ├------┤ F ├------┘ └----┘ └----┘References
- Zh10
Zhang et al. “An efficient subpixel image registration based on the phase-only correlations of image projections”, IEEE proceedings of the 10th international symposium on communications and information technologies, 2010.
- dhdt.processing.matching_tools_frequency_correlators.robust_corr(I1, I2)¶
match two imagery through fast robust correlation
- Parameters
I1 (numpy.array, size=(m,n), ndim=2) – array with intensities
I2 (numpy.array, size=(m,n), ndim=2) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
References
- Fi05
Fitch et al. “Fast robust correlation”, IEEE transactions on image processing vol. 14(8) pp. 1063-1073, 2005.
- Es07
Essannouni et al. “Adjustable SAD matching algorithm using frequency domain” Journal of real-time image processing, vol.1 pp.257-265, 2007.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = robust_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.sign_only_corr(I1, I2)¶
match two imagery through phase only correlation
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
C – displacement surface
- Return type
numpy.array, size=(m,n), real
References
- IK07
Ito & Kiya, “DCT sign-only correlation with application to image matching and the relationship with phase-only correlation”, IEEE international conference on acoustics, speech and signal processing, vol. 1, 2007.
- dhdt.processing.matching_tools_frequency_correlators.symmetric_phase_corr(I1, I2)¶
match two imagery through symmetric phase only correlation (SPOF) also known as Smoothed Coherence Transform (SCOT)
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [3].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
Notes
The matching equations are as follows:
\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1], \quad \mathcal{F}[\mathbf{I}_2]\]\[\mathbf{W} = 1 / \sqrt{||\mathbf{S}_1||||\mathbf{S}_2||}\]\[\mathbf{Q}_{12} = \mathbf{S}_1 [\mathbf{W}\mathbf{S}_2]^{\star}\]where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation
Schematically this looks like:
I2 ┌----┐ S2 ---┤ F ├---┬-------------------┐ └----┘ | ┌--------------┐W | ├-┤ 1/(|x||x|)^½ ├--× | └--------------┘ | Q12 ┌------┐ C12 | ×-----┤ F^-1 ├--- I1 ┌----┐ S1| | └------┘ ---┤ F ├---┴-------------------┘ └----┘If multiple bands are given, cross-spectral stacking is applied, see [AL22].
References
- NP93
Nikias & Petropoulou. “Higher order spectral analysis: a nonlinear signal processing framework”, Prentice hall. pp.313-322, 1993.
- We05
Wernet. “Symmetric phase only filtering: a new paradigm for DPIV data processing”, Measurement science and technology, vol.16 pp.601-618, 2005.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair >>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = symmetric_phase_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.upsampled_cross_corr(S1, S2, upsampling=2)¶
apply cros correlation, and upsample the correlation peak
- Parameters
S1 (numpy.array, size=(m,n), dtype=complex, ndim=2) – array with intensities
S2 (numpy.array, size=(m,n), dtype=complex, ndim=2) – array with intensities
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
dhdt.processing.matching_tools_frequency_correlators.pad_dft,dhdt.processing.matching_tools_frequency_correlators.upsample_dftReferences
- GS08
Guizar-Sicairo, et al. “Efficient subpixel image registration algorithms”, Applied optics, vol. 33 pp.156–158, 2008.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> di,dj = upsampled_cross_corr(im1, im2)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
- dhdt.processing.matching_tools_frequency_correlators.windrose_corr(I1, I2)¶
match two imagery through windrose phase correlation, see [KJ91].
- Parameters
I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].
I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities
- Returns
Q – cross-power spectrum
- Return type
numpy.array, size=(m,n), dtype=complex
See also
dhdt.processing.matching_tools_frequency_correlators.binary_orientation_corr,dhdt.processing.matching_tools_frequency_correlators.orientation_corr,dhdt.processing.matching_tools_frequency_correlators.phase_only_corrNotes
The matching equations are as follows:
\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]\[\mathbf{Q}_{12} = sign(\mathbf{S}_1) sign(\mathbf{S}_2)^{\star}\]where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation
References
- KJ91
Kumar & Juday, “Design of phase-only, binary phase-only, and complex ternary matched filters with increased signal-to-noise ratios for colored noise”, Optics letters, vol. 16(13) pp. 1025–1027, 1991.
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery” Science of Remote Sensing, vol.6 pp.100070, 2022.
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = windrose_corr(im1, im2) >>> C = np.fft.ifft2(Q) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
spatial correlators¶
- dhdt.processing.matching_tools_spatial_correlators.cosine_similarity(I1, I2)¶
estimate the similarity of orientation imagery, via their dot product. Going along via different names, such as Gradient Inner Product (GIP)[Gl08]_ or Cosine Similarity [De21].
- Parameters
I1 (numpy.array, size=(m,n), dtype=complex, {x ∈ ℂ | ║x║ = 1}) – image with intensities (template), given as complex orientation image
I2 (numpy.array, size=(k,l), dtype=complex, {x ∈ ℂ | ║x║ = 1}) – image with intensities (search space), as complex orientation image
- Returns
score
- Return type
numpy.array, dtype=float, range=0…+1
Examples
>>> import numpy as np >>> from scipy import ndimage >>> from dhdt.testing.matching_tools import create_sample_image_pair >>> from dhdt.generic.handler_im import get_grad_filters >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.processing.matching_tools_frequency_filters import normalize_power_spectrum
>>> _,im2,ti,tj,im1 = create_sample_image_pair(d=2**4, max_range=1)
complex orientation imagery are used here, so these need to be transformed
>>> H_x = get_grad_filters(ftype='kroon', tsize=3, order=1)[0] >>> H_y = np.transpose(H_x) >>> G1 = ndimage.convolve(im2, H_x) + 1j*ndimage.convolve(im2, H_y) >>> G2 = ndimage.convolve(im1_big, H_x) + 1j*ndimage.convolve(im1_big, H_y) >>> G1, G2 = normalize_power_spectrum(G1), normalize_power_spectrum(G2)
now processing that actual processing can be done
>>> C = normalized_cross_corr(G1,G2) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
References
- Gl08
Glocker et al. “Optical flow estimation with uncertainties through dynamic MRFs” IEEE conference on computer vision and pattern recognition, 2008.
- De21
Dematteis et al., “Comparison of digital image correlation methods and the impact of noise in geoscience applications” Remote sensing, vol.13(2), pp.327, 2021.
- dhdt.processing.matching_tools_spatial_correlators.cumulative_cross_corr(I1, I2)¶
doing normalized cross correlation on distance imagery
- Parameters
I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)
I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)
- Returns
ccs – similarity surface
- Return type
numpy.array
Notes
The following nomenclature is used:
ccs: cross correlation surface
- dhdt.processing.matching_tools_spatial_correlators.maximum_likelihood(I1, I2)¶
maximum likelihood texture tracking
- Parameters
I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)
I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)
- Returns
C – texture correspondence surface
- Return type
float
References
- Er09
Erten et al. “Glacier velocity monitoring by maximum likelihood texture tracking” IEEE transactions on geosciences and remote sensing, vol.47(2) pp.394–405, 2009.
- dhdt.processing.matching_tools_spatial_correlators.normalized_cross_corr(I1, I2)¶
simple normalized cross correlation
- Parameters
I1 (numpy.ndarray, ndim=2, size=(m,n)) – grid with intensities (template)
I2 (numpy.ndarray, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)
- Returns
ccs – similarity surface,
- Return type
numpy.array
Examples
>>> import numpy as np >>> from dhdt.processing.matching_tools import get_integer_peak_location >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> _,im2,ti,tj,im1 = create_sample_image_pair(d=2**4, max_range=1) >>> C = normalized_cross_corr(im1,im2) >>> di,dj,_,_ = get_integer_peak_location(C)
>>> assert(np.isclose(ti, di, atol=1)) >>> assert(np.isclose(ti, di, atol=1))
Notes
The following nomenclature is used:
ccs: cross correlation surface
- dhdt.processing.matching_tools_spatial_correlators.sum_sad_diff(I1, I2)¶
sum of absolute difference correlation
- Parameters
I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)
I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)
- Returns
sad – dissimilarity surface
- Return type
numpy.array, ndim=2
Notes
The following nomenclature is used:
sad: sum of absolute difference
- dhdt.processing.matching_tools_spatial_correlators.sum_sq_diff(I1, I2)¶
sum of squared difference correlation
- Parameters
I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)
I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)
- Returns
ssd – dissimilarity surface
- Return type
numpy.array, ndim=2
Notes
The following nomenclature is used:
ssd: sum of squared differences
- dhdt.processing.matching_tools_spatial_correlators.total_gradients(I1, I2)¶
calculate total gradients and its kurtosis, inspired by [Ch17]
- Parameters
I1 (numpy.ndarray, ndim=2, size=(m,n)) – grid with intensities (template)
I2 (numpy.ndarray, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)
- Returns
tg – similarity score surface
- Return type
numpy.array, ndim=2
References
- Ch17
Chen et al. “Normalized total gradients: a new measure for multispectral image registration” IEEE transactions on image processing, vol.27(3) pp.1297–1310, 2017.
- dhdt.processing.matching_tools_spatial_correlators.weighted_normalized_cross_correlation(I1, I2, W1=None, W2=None)¶
- Parameters
I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)
I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)
W1 (numpy.array, size=(m,n), dtype={boolean,float}) – weighting matrix for the template
W2 (numpy.array, size=(m,n), dtype={boolean,float}) – weighting matrix for the search range
- Returns
wncc (numpy.array) – weighted normalized cross-correlation
support (numpy.array, range=0…1) – amount of overlap present and used in both imagery
References
- AK21
Altena & Kääb “Quantifying river ice movement through a combination of European satellite monitoring services” International journal of applied Earth observation and geoinformation, vol.98 pp.102315, 2021.
phase angle estimators¶
- class dhdt.processing.matching_tools_frequency_subpixel.PlaneModel¶
Least squares estimator for phase plane.
Vectors/lines are parameterized using polar coordinates as functional model:
z = x * dx + y * dy
This estimator minimizes the squared distances from all points to the plane, independent of distance.
- estimate(data)¶
Estimate plane from data using least squares. :param data: N points with
(x, y)coordinates of vector, respectively. :type data: (N, 3) array- Returns
success – True, if model estimation succeeds.
- Return type
bool
- predict_xy(xy, params=None)¶
Predict vector using the estimated heading. :param xy: x,y-coordinates. :type xy: numpy.array :param params: Optional custom parameter set. :type params: (2, ) array, optional
- Returns
Q_hat – Predicted plane height at x,y-coordinates.
- Return type
numpy.array
- residuals(data)¶
Determine residuals of data to model For each point the shortest distance to the plane is returned.
- Parameters
data ((N, 3) array) – N points with x, y coordinates and z values, respectively.
- Returns
residuals – Residual for each data point.
- Return type
(N, ) array
- class dhdt.processing.matching_tools_frequency_subpixel.SawtoothModel¶
Least squares estimator for phase sawtooth.
- estimate(data, params_bound=0)¶
Estimate plane from data using least squares. :param data: N points with
(x, y)coordinates of vector, respectively. :type data: (N, 3) array :param params_bound: bound of the parameter space :type params_bound: float- Returns
success – True, if model estimation succeeds.
- Return type
bool
- predict_xy(xy, params=None)¶
Predict vector using the estimated heading. :param xy: x,y-coordinates. :type xy: numpy.array :param params: Optional custom parameter set. :type params: (2, ) array, optional
- Returns
Q_hat – Predicted plane height at x,y-coordinates.
- Return type
numpy.array
- residuals(data)¶
Determine residuals of data to model For each point the shortest distance to the plane is returned.
- Parameters
data ((N, 3) array) – N points with x, y coordinates and z values, respectively.
- Returns
residuals – Residual for each data point.
- Return type
(N, ) array
- dhdt.processing.matching_tools_frequency_subpixel.phase_difference(Q, W=array([], dtype=float64))¶
get displacement from phase plane through neighbouring vector difference
find slope of the phase plane through local difference of the phase angles
- Parameters
Q (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum
W (numpy.array, size=(m,n), dtype=boolean) – weigthing matrix
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
References
- Ka89
Kay, S. “A fast and accurate frequency estimator”, IEEE transactions on acoustics, speech and signal processing, vol.37(12) pp.1987-1990, 1989.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj,_,_ = phase_svd(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_difference_1d(Q, W=array([], dtype=float64), axis=0)¶
get displacement from phase plane along one axis through differencing
find slope of the phase plane through local difference of the phase angles, see [Ka89]
- Parameters
Q (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum
W (numpy.array, size=(m,n), dtype=boolean) – weigthing matrix
- Returns
dj – sub-pixel displacement
- Return type
float
See also
References
- Ka89
Kay, S. “A fast and accurate frequency estimator”, IEEE transactions on acoustics, speech and signal processing, vol.37(12) pp.1987-1990, 1989.
- dhdt.processing.matching_tools_frequency_subpixel.phase_gradient_descend(data, W=array([], dtype=float64), x_0=array([0., 0.]), learning_rate=1, n_iters=50)¶
get phase plane of cross-spectrum through principle component analysis
find slope of the phase plane through principle component analysis
- Parameters
data (numpy.array, size=(m*n,3), dtype=complex) – normalized cross spectrum
data – coordinate list with complex cross-sprectum at last collumn
W (numpy.array, size=(m*n,1), dtype=boolean) – index of data that is correct
W – list with classification of correct data
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj,_,_ = phase_gradient_descend(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_jac(Q, m, W=array([], dtype=float64), F1=array([], dtype=float64), F2=array([], dtype=float64), rank=2)¶
- Parameters
Q (numpy.array, size=(_,_), dtype=complex) – cross spectrum
m (numpy.array, size=(2,1), dtype=float) – displacement estimate, in pixel coordinate system
W (numpy.array, size=(m,n), dtype=float | boolean) – weigthing matrix, in a range of 0…1
F1 (np,array, size=(m,n), dtype=integer) – coordinate of the first axis from the Fourier spectrum.
F2 (np,array, size=(m,n), dtype=integer) – coordinate of the second axis from the Fourier spectrum
rank (TYPE, optional) – DESCRIPTION. The default is 2.
- Returns
dQdm – Jacobian of phase estimate
- Return type
numpy.array, size=(m,n)
- dhdt.processing.matching_tools_frequency_subpixel.phase_lsq(data, W=array([], dtype=float64))¶
get phase plane of cross-spectrum through least squares plane fitting
find slope of the phase plane through principle component analysis
- Parameters
data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum
numpy.array (or) – coordinate list with complex cross-sprectum at last
size=(m*n – coordinate list with complex cross-sprectum at last
3) – coordinate list with complex cross-sprectum at last
dtype=complex – coordinate list with complex cross-sprectum at last
W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct
numpy.array – list with classification of correct data
size=(m*n – list with classification of correct data
1) – list with classification of correct data
dtype=boolean – list with classification of correct data
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
phase_pca,phase_ransac,phase_houghExamples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj,_,_ = phase_lsq(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_pca(data, W=array([], dtype=float64))¶
get phase plane of cross-spectrum through principle component analysis
find slope of the phase plane through principle component analysis
- Parameters
data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum
numpy.array (or) – coordinate list with complex cross-sprectum at last
size=(m*n – coordinate list with complex cross-sprectum at last
3) – coordinate list with complex cross-sprectum at last
dtype=complex – coordinate list with complex cross-sprectum at last
W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct
numpy.array – list with classification of correct data
size=(m*n – list with classification of correct data
1) – list with classification of correct data
dtype=boolean – list with classification of correct data
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj = phase_pca(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_radon(Q, coord_system='ij')¶
get direction and magnitude from phase plane through Radon transform
find slope of the phase plane through single value decomposition
- Parameters
Q (numpy.array, size=(m,n), dtype=complex) – cross spectrum
- Returns
* di,dj (float, default) – cartesian displacement
* θ,ρ (float) – magnitude and direction of displacement
See also
References
- BF06
Balci & Foroosh. “Subpixel registration directly from the phase difference” EURASIP journal on advances in signal processing, pp.1-11, 2006.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj,_,_ = phase_radon(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_ransac(data, max_displacement=1, precision_threshold=0.05)¶
robustly fit plane using RANSAC algorithm
find slope of the phase plane through random sampling and consensus
- Parameters
data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum
numpy.array (or) – coordinate list with complex cross-sprectum at last
size=(m*n – coordinate list with complex cross-sprectum at last
3) – coordinate list with complex cross-sprectum at last
dtype=complex – coordinate list with complex cross-sprectum at last
- Returns
di (float) – sub-pixel displacement along vertical axis
dj (float) – sub-pixel displacement along horizontal axis
References
- FB81
Fischler & Bolles. “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography” Communications of the ACM vol.24(6) pp.381-395, 1981.
- To15
Tong et al. “A novel subpixel phase correlation method using singular value decomposition and unified random sample consensus” IEEE transactions on geoscience and remote sensing vol.53(8) pp.4143-4156, 2015.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj,_,_ = phase_ransac(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_secant(data, W=array([], dtype=float64), x_0=array([0., 0.]))¶
get phase plane of cross-spectrum through secant
find slope of the phase plane through secant method (or Newton’s method) in multiple dimensions it is known as the Broyden’s method.
- Parameters
data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum
numpy.array (or) – coordinate list with complex cross-sprectum at last collumn
size=(m*n – coordinate list with complex cross-sprectum at last collumn
3) – coordinate list with complex cross-sprectum at last collumn
dtype=complex – coordinate list with complex cross-sprectum at last collumn
W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct
numpy.array – list with classification of correct data
size=(m*n – list with classification of correct data
1) – list with classification of correct data
dtype=boolean – list with classification of correct data
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
References
- Br65
Broyden, “A class of methods for solving nonlinear simultaneous equations” Mathematics and computation. vol.19(92) pp.577–593, 1965.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj,_,_ = phase_secant(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_slope_1d(t, rad=0.1)¶
estimate the slope and intercept for one-dimensional signal
- Parameters
t (numpy.array, size=(m,1), dtype=complex) – angle values.
rad (float, range=(0.0,0.5)) – radial inclusion, seen from the center
- Returns
x_hat – estimated slope and intercept.
- Return type
numpy.array, size=(2,1)
See also
- dhdt.processing.matching_tools_frequency_subpixel.phase_svd(Q, W, rad=0.1)¶
get phase plane of cross-spectrum through single value decomposition
find slope of the phase plane through single value decomposition, see [Ho03]
- Parameters
Q (numpy.array, size=(m,n), dtype=complex) – cross spectrum
W (numpy.array, size=(m,n), dtype=float) – weigthing matrix
rad (float, range=(0.0,0.5)) – radial inclusion, seen from the center
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
References
- Ho03
Hoge, W.S. “A subspace identification extension to the phase correlation method”, IEEE transactions on medical imaging, vol.22(2) pp.277-280, 2003.
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> di,dj,_,_ = phase_svd(Q)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.phase_tpss(Q, W, m, p=0.0001, k=4, j=5, n=3)¶
get phase plane of cross-spectrum through two point step size iteration
find slope of the phase plane through two point step size for phase correlation minimization [BB88] which is implemented by [Le07]
- Parameters
Q (numpy.ndarray, size=(_,_), dtype=complex) – cross spectrum
m0 (numpy.ndarray, size=(2,1)) – initial displacement estimate
p (float, default=1e4) – closing error threshold
k (integer, default=4) – number of refinements in iteration
j (integer, default=5) – number of sub routines during an estimation
n (integer, default=3) – mask convergence factor
- Returns
di,dj (float, size=(2,1)) – sub-pixel displacement
snr (float) – signal-to-noise ratio
See also
References
- BB88
Barzilai & Borwein. “Two-point step size gradient methods”, IMA journal of numerical analysis. vol.8 pp.141–148, 1988.
- Le07
Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.
- dhdt.processing.matching_tools_frequency_subpixel.phase_weighted_pca(Q, W)¶
get phase plane of cross-spectrum through principle component analysis
find slope of the phase plane through principle component analysis
- Parameters
data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum
numpy.array (or) – coordinate list with complex cross-sprectum at last
size=(m*n – coordinate list with complex cross-sprectum at last
3) – coordinate list with complex cross-sprectum at last
dtype=complex – coordinate list with complex cross-sprectum at last
W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct
numpy.array – list with classification of correct data
size=(m*n – list with classification of correct data
1) – list with classification of correct data
dtype=boolean – list with classification of correct data
- Returns
di,dj – sub-pixel displacement
- Return type
float
See also
Examples
>>> import numpy as np >>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1) >>> Q = phase_corr(im1, im2) >>> W = gaussian_mask(Q) >>> di,dj,_,_ = phase_weighted_pca(Q, W)
>>> assert(np.isclose(ti, di, atol=.2)) >>> assert(np.isclose(tj, dj, atol=.2))
- dhdt.processing.matching_tools_frequency_subpixel.ransac(data, model_class, min_samples, residual_threshold, is_data_valid=None, is_model_valid=None, max_trials=100, stop_sample_num=inf, stop_residuals_sum=0, stop_probability=1, random_state=None, initial_inliers=None, params_bounds=0)¶
Fit a model to data with the RANSAC (random sample consensus) algorithm. RANSAC is an iterative algorithm for the robust estimation of parameters from a subset of inliers from the complete data set. Each iteration performs the following tasks: 1. Select min_samples random samples from the original data and check
whether the set of data is valid (see is_data_valid).
Estimate a model to the random subset (model_cls.estimate(*data[random_subset]) and check whether the estimated model is valid (see is_model_valid).
Classify all data as inliers or outliers by calculating the residuals to the estimated model (model_cls.residuals(*data)) - all data samples with residuals smaller than the residual_threshold are considered as inliers.
Save estimated model as best model if number of inlier samples is maximal. In case the current estimated model has the same number of inliers, it is only considered as the best model if it has less sum of residuals.
These steps are performed either a maximum number of times or until one of the special stop criteria are met. The final model is estimated using all inlier samples of the previously determined best model. :param data: Data set to which the model is fitted, where N is the number of data
points and the remaining dimension are depending on model requirements. If the model class requires multiple input data arrays (e.g. source and destination coordinates of
skimage.transform.AffineTransform), they can be optionally passed as tuple or list. Note, that in this case the functionsestimate(*data),residuals(*data),is_model_valid(model, *random_data)andis_data_valid(*random_data)must all take each data array as separate arguments.- Parameters
model_class (object) –
- Object with the following object methods:
success = estimate(*data)residuals(*data)
where success indicates whether the model estimation succeeded (True or None for success, False for failure).
min_samples (int in range (0, N)) – The minimum number of data points to fit a model to.
residual_threshold (float larger than 0) – Maximum distance for a data point to be classified as an inlier.
is_data_valid (function, optional) – This function is called with the randomly selected data before the model is fitted to it: is_data_valid(*random_data).
is_model_valid (function, optional) – This function is called with the estimated model and the randomly selected data: is_model_valid(model, *random_data), .
max_trials (int, optional) – Maximum number of iterations for random sample selection.
stop_sample_num (int, optional) – Stop iteration if at least this number of inliers are found.
stop_residuals_sum (float, optional) – Stop iteration if sum of residuals is less than or equal to this threshold.
stop_probability (float in range [0, 1], optional) –
RANSAC iteration stops if at least one outlier-free set of the training data is sampled with
probability >= stop_probability, depending on the current best model’s inlier ratio and the number of trials. This requires to generate at least N samples (trials):N >= log(1 - probability) / log(1 - e**m)
where the probability (confidence) is typically set to a high value such as 0.99, e is the current fraction of inliers w.r.t. the total number of samples, and m is the min_samples value.
random_state ({None, int, numpy.random.Generator}, optional) – If random_state is None the numpy.random.Generator singleton is used. If random_state is an int, a new
Generatorinstance is used, seeded with random_state. If random_state is already aGeneratorinstance then that instance is used.initial_inliers (array-like of bool, shape (N,), optional) – Initial samples selection for model estimation
- Returns
model (object) – Best model with largest consensus set.
inliers ((N, ) array) – Boolean mask of inliers classified as
True.
References
- 1
“RANSAC”, Wikipedia, https://en.wikipedia.org/wiki/RANSAC
Examples
Generate ellipse data without tilt and add noise:
>>> t = np.linspace(0, 2 * np.pi, 50) >>> xc, yc = 20, 30 >>> a, b = 5, 10 >>> x = xc + a * np.cos(t) >>> y = yc + b * np.sin(t) >>> data = np.column_stack([x, y]) >>> rng = np.random.default_rng(203560) # do not copy this value >>> data += rng.normal(size=data.shape)
Add some faulty data:
>>> data[0] = (100, 100) >>> data[1] = (110, 120) >>> data[2] = (120, 130) >>> data[3] = (140, 130)
Estimate ellipse model using all available data:
>>> model = EllipseModel() >>> model.estimate(data) True >>> np.round(model.params) array([ 72., 75., 77., 14., 1.])
Estimate ellipse model using RANSAC:
>>> ransac_model, inliers = ransac( ... data, EllipseModel, 20, 3, max_trials=50 ... ) >>> abs(np.round(ransac_model.params)) array([20., 30., 10., 6., 2.])
>>> inliers array([False, False, False, False, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True], dtype=bool)
>>> sum(inliers) > 40 True
RANSAC can be used to robustly estimate a geometric transformation. In this section, we also show how to use a proportion of the total samples, rather than an absolute number.
>>> from skimage.transform import SimilarityTransform >>> rng = np.random.default_rng() >>> src = 100 * rng.random((50, 2)) >>> model0 = SimilarityTransform(scale=0.5, rotation=1, ... translation=(10, 20)) >>> dst = model0(src) >>> dst[0] = (10000, 10000) >>> dst[1] = (-100, 100) >>> dst[2] = (50, 50) >>> ratio = 0.5 # use half of the samples >>> min_samples = int(ratio * len(src)) >>> model, inliers = ransac((src, dst), SimilarityTransform, min_samples, ... 10, ... initial_inliers=np.ones(len(src), dtype=bool))
>>> inliers array([False, False, False, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True])
subpixel localization¶
- dhdt.processing.matching_tools_spatial_subpixel.get_top_2d_gaussian(C, top=array([], dtype=float64))¶
find location of highest score using 2D Gaussian
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- NH05
Nobach & Honkanen, “Two-dimensional Gaussian regression for sub-pixel displacement estimation in particle image velocimetry or particle position estimation in particle tracking velocimetry”, Experiments in fluids, vol.38 pp.511-515, 2005
- dhdt.processing.matching_tools_spatial_subpixel.get_top_birchfield(C, top=array([], dtype=float64))¶
find location of highest score along each axis
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- BT99
Birchfield & Tomasi. “Depth discontinuities by pixel-to-pixel stereo” International journal of computer vision, vol. 35(3) pp.269-293, 1999.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_blais(C, top=array([], dtype=float64))¶
find location of highest score through forth order filter
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- BR86
Blais & Rioux, “Real-time numerical peak detector” Signal processing vol.11 pp.145-155, 1986.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_blue(C, ds=1, top=array([], dtype=float64))¶
find top of correlation peak through best linear unbiased estimation
References
- 1
- dhdt.processing.matching_tools_spatial_subpixel.get_top_centroid(C, top=array([], dtype=float64))¶
find location of highest score through 1D centorid fit
- Parameters
C (numpy.ndarray, size=(_,_)) – similarity surface
top (numpy.ndarray, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak ddi : float estimated subpixel location on the vertical axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- Ra18
Raffel et al. “Particle Image Velocimetry” Ch.6 pp.184, 2018.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_equiangular(C, top=array([], dtype=float64))¶
find location of highest score along each axis by equiangular line
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the crossing
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- SO05
Shimizu & Okutomi. “Sub-pixel estimation error cancellation on area-based matching” International journal of computer vision, vol.63(3), pp.207–224, 2005.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_esinc(C, ds=1, top=array([], dtype=float64))¶
find location of highest score using exponential esinc function
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- AV06
Argyriou & Vlachos, “A study of sub-pixel motion estimation using phase correlation”, proceedings of the British machine vision conference, 2006
- dhdt.processing.matching_tools_spatial_subpixel.get_top_gaussian(C, top=array([], dtype=float64))¶
find location of highest score through 1D gaussian fit
- Parameters
C (numpy.ndarray, size=(_,_)) – similarity surface
top (numpy.ndarray, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak ddi : float estimated subpixel location on the vertical axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- WG91
Willert & Gharib, “Digital particle image velocimetry”, Experiments in fluids, vol.10 pp.181-193, 1991.
- AV06
Argyriou & Vlachos, “A Study of sub-pixel motion estimation using phase correlation” Proceeding of the British machine vision conference, pp.387-396, 2006.
- Ra18
Raffel et al. “Particle Image Velocimetry” Ch.6 pp.184 2018.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_mass(C, top=array([], dtype=float64))¶
find location of highest score through 1D center of mass
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- FN96
Fisher & Naidu, “A Comparison of algorithms for subpixel peak detection” in Image Technology - Advances in image processing, multimedia and machine vision pp.385-404, 1996.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_moment(C, ds=1, top=array([], dtype=float64))¶
find location of highest score through bicubic fitting
- Parameters
C (numpy.ndarray, size=(_,_)) – similarity surface
ds (integer, default=1) – size of the radius to use neighboring information
top (numpy.ndarray, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
References
- Fe12
Feng et al. “A subpixel registration algorithm for low PSNR images” IEEE international conference on advanced computational intelligence, pp.626-630, 2012.
- MG15
Messerli & Grinstad, “Image georectification and feature tracking toolbox: ImGRAFT” Geoscientific instrumentation, methods and data systems, vol.4(1) pp.23-34, 2015.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_parabolic(C, top=array([], dtype=float64), ds=1)¶
find location of highest score through 1D parabolic fit
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- AV06
Argyriou & Vlachos, “A Study of sub-pixel motion estimation using phase correlation” Proceeding of the British machine vision conference, pp.387-396), 2006.
- Ra18
Raffel et al. “Particle Image Velocimetry” Ch.6 pp.184 2018.
- Br21
Bradley, “Sub-pixel registration for low cost, low dosage, X-ray phase contrast imaging”, 2021.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_paraboloid(C, top=array([], dtype=float64))¶
find location of highest score using paraboloid
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- Pa20
Pallotta et al. “Subpixel SAR image registration through parabolic interpolation of the 2-D cross correlation”, IEEE transactions on geoscience and remote sensing, vol.58(6) pp.4132–4144, 2020.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_ren(C, top=array([], dtype=float64))¶
find location of highest score
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- Re10
Ren et al. “High-accuracy sub-pixel motion estimation from noisy images in Fourier domain.” IEEE transactions on image processing, vol.19(5) pp.1379-1384, 2010.
- dhdt.processing.matching_tools_spatial_subpixel.get_top_triangular(C, top=array([], dtype=float64))¶
find location of highest score through triangular fit
- Parameters
C (numpy.array, size=(_,_)) – similarity surface
top (numpy.array, size=(1,2)) – location of the maximum score
- Returns
ddi (float) – estimated subpixel location on the vertical axis of the peak
ddj (float) – estimated subpixel location on the horizontal axis of the peak
i_int (integer) – location of highest score on the vertical axis
j_int (integer) – location of highest score on the horizontal axis
Notes
- OC91
Olsen & Coombs, “Real-time vergence control for binocular robots” International journal of computer vision, vol.7(1), pp.67-89, 1991.
frequency metrics¶
- dhdt.processing.matching_tools_frequency_metrics.get_phase_metric(Q, di, dj, metric='phase_fit')¶
redistribution function, to get a metric from a matching score surface
- Parameters
Q (numpy.ndarray, size=(m,n), type=complex) – phase angles of the cross-correlation spectra
di ({float,integer}) – estimated displacement
dj ({float,integer}) – estimated displacement
metric (string) – abbreviation for the metric type to be calculated, for the options see
list_pahse_metrics()for the options
- Returns
score – metric of the matching phase angle in relation to its surface
- Return type
float
- dhdt.processing.matching_tools_frequency_metrics.list_phase_metrics()¶
list the abbreviations of the different implemented phase metrics, there are:
'phase_fit': the absolute score'phase_sup': the primary peak ratio i.r.t. the second peak
See also
- dhdt.processing.matching_tools_frequency_metrics.phase_fitness(Q, di, dj, norm=2)¶
estimate the goodness of fit of the phase diffence
- Parameters
Q (numpy.ndarray, size={(m,n), (k,3)}) – cross-spectrum
di ({float,integer}) – estimated displacement
dj ({float,integer}) – estimated displacement
norm (integer) – norm for the difference
- Returns
fitness – goodness of fit, see also [1]
- Return type
float, range=0…1
See also
References
- Le07
Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.
- dhdt.processing.matching_tools_frequency_metrics.phase_support(Q, di, dj, thres=1.4826)¶
estimate the support of the fit for the phase differences, following the implementation by [AL22].
- Parameters
Q (numpy.ndarray, size={(m,n), (k,3)}) – cross-spectrum
di ({float,integer}) – estimated displacement
dj ({float,integer}) – estimated displacement
norm (integer) – norm for the difference
- Returns
support – ammount of support, see also [1]
- Return type
float, range=0…1
See also
References
- AL22
Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.
- dhdt.processing.matching_tools_frequency_metrics.signal_to_noise(Q, C, norm=2)¶
calculate the signal to noise from a theoretical and an experimental cross-spectrum, see [Le07].
- Parameters
Q (numpy.ndarray, size=(m,n), dtype=complex) – cross-spectrum
C (numpy.ndarray, size=(m,n), dtype=complex) – phase plane
norm (integer) – norm for the difference
- Returns
snr – signal to noise ratio, see also [1]
- Return type
float, range=0…1
See also
References
- Le07
Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.
correlation metrics¶
- dhdt.processing.matching_tools_spatial_metrics.entropy_corr(C)¶
metric for uniqueness of the correlation estimate
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
- Returns
disorder – entropy of the correlation surface
- Return type
float
See also
dhdt.processing.matching_tools_spatial_metrics.peak_confidence,dhdt.processing.matching_tools_spatial_metrics.peak_to_noise,dhdt.processing.matching_tools_spatial_metrics.peak_corr_energy,dhdt.processing.matching_tools_spatial_metrics.peak_rms_ratio,dhdt.processing.matching_tools_spatial_metrics.num_of_peaks,dhdt.processing.matching_tools_spatial_metrics.peak_winner_margin,dhdt.processing.matching_tools_spatial_metrics.primary_peak_margin,dhdt.processing.matching_tools_spatial_metrics.primary_peak_ratioReferences
- SS98
Scharstein & Szeliski, “Stereo matching with nonlinear diffusion” International journal of computer vision, vol.28(2) pp.155-174, 1998.
- Xu14
Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.
- St26
Sturges, “The choice of a class interval”. Journal of the american statistical association. vol.21(153) pp.65–66, 1926.
- dhdt.processing.matching_tools_spatial_metrics.gauss_spread(C, intI, intJ, dI, dJ, est='dist')¶
estimate an oriented gaussian function through the vicinity of the correlation function. See for a full describtion [Al21].
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – array with correlation values
intI (integer, {x ∈ ℕ | x ≥ 0}) – location in rows of highest value in search space
intJ (integer, {x ∈ ℕ | x ≥ 0}) – location in collumns of highest value in search space
dI (float) – sub-pixel bias of top location along the row axis
dJ (float) – sub-pixel bias of top location along the collumn axis
est (string, default='dist') –
‘dist’ do weighted least sqaures dependent on the distance from the top location
otherwise : ordinary least squares
- Returns
cov_ii, cov_jj (float) – standard deviation of the vertical and horizontal axis
ρ (float) – orientation of the Gaussian
hess (numpy.ndarray, size=(4)) – estimate of the least squares computation
frac – scaling of the correlation peak
Notes
The image frame is used in this function, hence the difference is:
coordinate | +--------> collumns system 'ij'| | | | | j | image frame --------+--------> | | | | v | i rows vReferences
- Al21
Altena et al. “Correlation dispersion as a measure to better estimate uncertainty of remotely sensed glacier displacements” The cryosphere, vol.16(6) pp.2285-2300, 2021.
- dhdt.processing.matching_tools_spatial_metrics.get_correlation_metric(C, metric='peak_abs')¶
redistribution function, to get a metric from a matching score surface
- Parameters
C (numpy.ndarray, size=(m,n)) – grid with correlation scores
metric ({'peak_ratio', 'peak_rms', 'peak_ener', 'peak_nois', 'peak_conf',) – ‘peak_entr’, ‘peak_abs’, ‘peak_marg’, ‘peak_win’, ‘peak_num’} abbreviation for the metric type to be calculated, for the options see
list_matching_metrics()for the options
- Returns
score – metric of the matching peak in relation to its surface
- Return type
float
- dhdt.processing.matching_tools_spatial_metrics.hessian_spread(C, intI, intJ)¶
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – array with correlation values
intI (integer, {x ∈ ℕ | x ≥ 0}) – location in rows of highest value in search space
intJ (integer, {x ∈ ℕ | x ≥ 0}) – location in collumns of highest value in search space
- Returns
cov_ii (float, real, unit=pixel) – first diagonal element of the co-variance matrix
cov_jj (float, real, unit=pixel) – second diagonal element of the co-variance matrix
cov_ij (float, real, unit=pixel) – off-diagonal element of the co-variance matrix
Notes
The image frame is used in this function, hence the difference is:
coordinate | +--------> collumns system 'ij'| | | | | j | image frame --------+--------> | | | | v | i rows vReferences
- Ro04
Rosen et al. “Updated repeat orbit interferometry package released” EOS, vol.85(5) pp.47, 2004.
- Ca11
Casu et al. “Deformation time-series generation in areas characterized by large displacement dynamics: The SAR amplitude pixel-offset SBAS Technique” IEEE transactions in geoscience and remote sensing vol.49(7) pp.2752–2763, 2011.
- dhdt.processing.matching_tools_spatial_metrics.intensity_disparity(I1, I2)¶
- Parameters
I1 (numpy.ndarray, size=(m,n), dtype=float) – template with intensities
I2 (numpy.ndarray, size=(m,n), dtype=float) – template with intensities
- Returns
sigma_dI – temporal intensity variation
- Return type
float
References
- Sc13
Sciacchitano et al. “PIV uncertainty quantification by image matching” Measurement science and technology, vol.24 pp.045302, 2013.
- dhdt.processing.matching_tools_spatial_metrics.list_matching_metrics()¶
list the abbreviations of the different implemented correlation metrics
- The following metrics are implemented:
‘peak_abs’ : the absolute score
'peak_ratio': the primary peak ratio i.r.t. the second peak'peak_rms': the peak ratio i.r.t. the root mean square error'peak_ener': the peaks’ energy'peak_noise': the peak score i.r.t. to the noise level'peak_conf': the peak confidence'peak_entr': the peaks’ entropy'peak_num': the number of correlation peaks’'peak_marg': the dominance of the peak
See also
- dhdt.processing.matching_tools_spatial_metrics.num_of_peaks(C, filtering=True)¶
metric for the uniqueness of a match
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
filtering (boolean) – apply low-pass filtering, otherwise noise is also seen as a peak, see [2] for more motivation
- Returns
ppm – margin between primary and secondary peak
- Return type
float
See also
entropy_corr,peak_confidence,peak_to_noise,peak_corr_energy,peak_rms_ratio,peak_winner_margin,primary_peak_margin,primary_peak_ratioReferences
- Le07
Lefebvre et al., “A colour correlation-based stereo matching using 1D windows” IEEE conference on signal-image technologies and internet-based system, pp.702-710, 2007.
- HM12
Hu & Mordohai, “A quantitative evaluation of confidence measures for stereo vision” IEEE transactions on pattern analysis and machine intelligence, vol.34(11) pp.2121-2133, 2012.
- dhdt.processing.matching_tools_spatial_metrics.peak_confidence(C, radius=1)¶
metric for steepness of a correlation peak
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
radius (integer, {x ∈ ℕ | x ≥ 1}, default=1) – how far away should the sampling set be taken
- Returns
q – confidence measure
- Return type
float
See also
entropy_corr,peak_to_noise,peak_corr_energy,peak_rms_ratio,num_of_peaks,peak_winner_margin,primary_peak_margin,primary_peak_ratioReferences
- Er09
Erten et al. “Glacier velocity monitoring by maximum likelihood texture tracking” IEEE transactions on geosciences and remote sensing, vol.47(2) pp.394–405, 2009.
- Er13
Erten “Glacier velocity estimation by means of polarimetric similarity measure” IEEE transactions on geosciences and remote sensing, vol.51 pp.3319–3327, 2013.
- dhdt.processing.matching_tools_spatial_metrics.peak_corr_energy(C)¶
metric for uniqueness of the correlation estimate, also known as, “signal to noise”
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
- Returns
prmsr – peak to correlation energy
- Return type
float
See also
entropy_corr,peak_confidence,peak_to_noise,peak_rms_ratio,num_of_peaks,peak_winner_margin,primary_peak_margin,primary_peak_ratioReferences
- Xu14
Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.
- dhdt.processing.matching_tools_spatial_metrics.peak_rms_ratio(C)¶
metric for uniqueness of the correlation estimate this is a similar metric as implemented in ROI_PAC [RoXX]
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
- Returns
prmsr – peak to root mean square ratio
- Return type
float
See also
entropy_corr,peak_confidence,peak_to_noise,peak_corr_energy,num_of_peaks,peak_winner_margin,primary_peak_margin,primary_peak_ratioReferences
- Xu14
Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.
- Ro04
Rosen et al. “Updated repeat orbit interferometry package released” EOS, vol.85(5) pp.47, 2004.
- dhdt.processing.matching_tools_spatial_metrics.peak_to_noise(C)¶
metric for uniqueness of the correlation estimate
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
- Returns
snr – signal to noise ratio
- Return type
float
See also
entropy_corr,peak_confidence,peak_corr_energy,peak_rms_ratio,num_of_peaks,peak_winner_margin,primary_peak_margin,primary_peak_ratioReferences
- dL07
de Lange et al. “Improvement of satellite radar feature tracking for ice velocity derivation by spatial frequency filtering” IEEE transactions on geosciences and remote sensing, vol.45(7) pp.2309–2318, 2007.
- HK12
Heid & Kääb “Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery.” Remote sensing of environment vol.118 pp.339-355, 2012.
- dhdt.processing.matching_tools_spatial_metrics.peak_winner_margin(C)¶
metric for dominance of the correlation estimate in relation to other candidates and their surrounding scores. See also [SS98].
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
- Returns
ppm – margin between primary and secondary peak
- Return type
float
See also
entropy_corr,peak_confidence,peak_to_noise,peak_corr_energy,peak_rms_ratio,num_of_peaks,primary_peak_margin,primary_peak_ratioReferences
- SS98
Scharstein & Szeliski, “Stereo matching with nonlinear diffusion” International journal of computer vision, vol.28(2) pp.155-174, 1998.
- dhdt.processing.matching_tools_spatial_metrics.primary_peak_margin(C)¶
metric for dominance of the correlation estimate in relation to other candidates
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
- Returns
ppm – margin between primary and secondary peak
- Return type
float
See also
entropy_corr,peak_confidence,peak_to_noise,peak_corr_energy,peak_rms_ratio,num_of_peaks,peak_winner_margin,primary_peak_ratioReferences
- HM12
Hu & Mordohai, “A quantitative evaluation of confidence measures for stereo vision” IEEE transactions on pattern analysis and machine intelligence, vol.34(11) pp.2121-2133, 2012.
- dhdt.processing.matching_tools_spatial_metrics.primary_peak_ratio(C)¶
metric for uniqueness of the correlation estimate also known as, “Cross correlation peak ratio”
- Parameters
C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface
- Returns
ppr – peak ratio between primary and secondary peak
- Return type
float, range={0..1}
See also
entropy_corr,peak_confidence,peak_to_noise,peak_corr_energy,peak_rms_ratio,num_of_peaks,peak_winner_margin,primary_peak_marginReferences
- KA90
Keane & Adrian, “Optimization of particle image velocimeters. I. Double pulsed systems” Measurement science and technology. vol.1 pp.1202, 1990.
- CV13
Charonk & Vlachos “Estimation of uncertainty bounds for individual particle image velocimetry measurements from cross-correlation peak ratio” Measurement science and technology. vol.24 pp.065301, 2013.
- Xu14
Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.
postprocessing¶
displacement filters¶
functions of use for the extraction of information of displacement fields, or filtering thereof.
- dhdt.postprocessing.displacement_filters.local_infilling_filter(Z, tsize=5)¶
apply local inverse distance weighting, on no-data values in an array
- Parameters
Z (numpy.array, size=(m,n), dtype=float) – intensity array with no-data values
tsize (integer, {x ∈ ℕ | x ≥ 1}, unit=pixel) – size of the neighborhood
- Returns
Z_new – intensity array with interpolated values at the data gap locations
- Return type
numpy.array, size=(m,n), dtype=float
References
- Jo02
Joughin “Ice-sheet velocity mapping: a combined interferometric and speckle-tracking approach”, Annuals of glaciology vol.34 pp.195-201, 2002.
- Jo17
Joughin et al. “Greenland ice mapping project 2 (GIMP-2) algorithm theoretical basis document”, Making earth system data records for use in research environment (MEaSUREs) documentation, 2017.
- dhdt.postprocessing.displacement_filters.local_variance(V, tsize=5)¶
local non-linear variance calculation
- Parameters
V (numpy.array, size=(m,n), dtype=float) – array with one velocity component, all algorithms are indepent of their axis.
sig_V (numpy.array, size=(m,n), dtype=float) – statistical local variance, based on the procedure described in [JoXX], [JoYY]
References
- Jo02
Joughin “Ice-sheet velocity mapping: a combined interferometric and speckle-tracking approach”, Annuals of glaciology vol.34 pp.195-201, 2002.
- Jo17
Joughin et al. “Greenland ice mapping project 2 (GIMP-2) algorithm theoretical basis document”, Making earth system data records for use in research environment (MEaSUREs) documentation, 2017.
- dhdt.postprocessing.displacement_filters.nan_resistant_filter(Z, kernel)¶
operate zonal filter on image with missing values. Energy within the kernel that fals on missing instances is distributed over other elements. Hence, regions with no data are shrinking, instead of growing. For a more detailed describtion see [Al22]
- Parameters
Z ({numpy.ndarray, numpy.masked.array}, size=(m,n)) – grid with values
kernel (numpy.ndarray, size=(k,l)) – convolution kernel
Returns –
Z_new (numpy.ndarray, size=(m,n)) – convolved estimate
References
- Al22
Altena et al. “Correlation dispersion as a measure to better estimate uncertainty of remotely sensed glacier displacements”, The Crysophere, vol.16(6) pp.2285–2300, 2022.
solar tools¶
functions of use for illumination calculations
- dhdt.postprocessing.solar_tools.az_to_sun_vector(az, indexing='ij')¶
transform azimuth angle to 2D-unit vector
- Parameters
az (float, unit=degrees) – azimuth of sun.
indexing ({‘xy’, ‘ij’}) –
“xy” : using map coordinates
”ij” : using local image coordinates
- Returns
sun – unit vector in the direction of the sun.
- Return type
numpy.array, size=(2,1), range=0…1
See also
Notes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--┼--> + | ┼----> East & xThe angles related to the suns’ heading are as follows:
surface normal * sun ^ ^ / | | / ├-- zenith angle | / | / | /| |/ |/ | elevation angle └---- └--┴---
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | i | x --------┼--------> --------┼--------> | | | | image | j map | based v based |
- dhdt.postprocessing.solar_tools.make_doppler_range(Z, az, zn, Lambertian=True, spac=10)¶
- Parameters
Z (numpy.array, unit=meters) – array with elevation values
az (float, unit=degrees, range=-180...+180) – flight orientation of the satellite
zn ({float,array}, unit=degrees, range=0...+90) – illumination angle from the satellite
Notes
- dhdt.postprocessing.solar_tools.make_shading(Z, az, zn, spac=10)¶
create synthetic shading image from given sun angles
A simple Lambertian reflection model is used here.
- Parameters
Z (numpy.array, size=(m,n), dtype={integer,float}, unit=meter) – grid with elevation data
az (float, unit=degrees) – azimuth angle
zn (float, unit=degrees) – zenith angle
spac ({float, tuple}, default=10, unit=meter) – resolution of the square grid. Optionally the geoTransform is given.
- Returns
Sh – estimated shading grid
- Return type
numpy.array, size=(m,n), dtype=float, range=0…1
Notes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--┼--> + | ┼----> East & xThe angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / ├-- zenith angle | / | / | /| |/ |/ | elevation angle └---- └--┴---
- dhdt.postprocessing.solar_tools.make_shading_minnaert(Z, az, zn, k=1, spac=10)¶
create synthetic shading image from given sun angles
A simple Minnaert reflection model is used here.
- Parameters
Z (numpy.array, size=(m,n), dtype={integer,float}, unit=meter) – grid with elevation data
az (float, unit=degrees) – azimuth angle
zn (float, unit=degrees) – zenith angle
spac (float, default=10, unit=meter) – resolution of the square grid.
- Returns
Sh – estimated shading grid
- Return type
numpy.array, size=(m,n), dtype=float, range=0…1
Notes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--┼--> + | ┼----> East & xThe angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / ├-- zenith angle | / | / | /| |/ |/ | elevation angle └---- └--┴---
- dhdt.postprocessing.solar_tools.make_shadowing(Z, az, zn, spac=10, weights=None, radiation=False, border='copy')¶
create synthetic shadow image from given sun angles
- Parameters
Z (numpy.array, size=(m,n), dtype={integer,float}) – grid with elevation data
az (float, unit=degrees) – azimuth angle
zn ({float, numpy.array}, unit=degrees) – zenith angle
spac (float, optional) – resolution of the square grid. The default is 10.
weights (numpy.array) – if the zenith angle is a vector, it can be used to get a weighted sum
border (string, {'copy','zeros'}) – specify what to do at the border, since it is a zonal operation
- Returns
Sw – estimated shadow grid
- Return type
numpy.array, size=(m,n), dtype=bool
Notes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--┼--> + | ┼----> East & xThe angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / ├-- zenith angle | / | / | /| |/ |/ | elevation angle └---- └--┴---
- dhdt.postprocessing.solar_tools.sun_angles_to_vector(az, zn, indexing='ij')¶
transform azimuth and zenith angle to 3D-unit vector
- Parameters
az ({float, numpy.ndarray}, unit=degrees) – azimuth angle of sun.
zn ({float, numpy.ndarray}, unit=degrees) – zenith angle of sun.
indexing ({‘xy’, ‘ij’}) –
“xy” : using map coordinates
”ij” : using local image coordinates
- Returns
sun – unit vector in the direction of the sun.
- Return type
numpy.array, size=(3,1), dtype=float, range=0…1
See also
dhdt.postprocessing.solar_tools.az_to_sun_vector,dhdt.postprocessing.solar_tools.vector_to_sun_anglesNotes
The azimuth angle declared in the following coordinate frame:
^ North & y | - <--┼--> + | ┼----> East & xThe angles related to the sun are as follows:
surface normal * sun ^ ^ / | | / ├-- zenith angle | / | / | /| |/ |/ | elevation angle └---- surface ----- └--┴---
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | i | x --------┼--------> --------┼--------> | | | | image | j map | based v based |
- dhdt.postprocessing.solar_tools.vector_to_sun_angles(sun, indexing='ij')¶
transform 3D-unit vector to azimuth and zenith angle
- Parameters
sun ({float,numpy.array}, size=(m,), dtype=float, range=0...1) – the three different unit vectors in the direction of the sun.
- Returns
az ({float, numpy.ndarray}, unit=degrees) – azimuth angle of sun.
zn ({float, numpy.ndarray}, unit=degrees) – zenith angle of sun.
The angles related to the suns’ heading are as follows –
surface normal * sun ^ ^ / | | / ├-- zenith angle | / | / | /| |/ |/ | elevation angle └---- └--┴---
terrain tools¶
- dhdt.postprocessing.terrain_tools.d8_catchment(V_x, V_y, geoTransform, x, y, disperse=True, flat=True)¶
estimate the catchment behind a seed point
- Parameters
V_x (numpy.ndarray, size=(m,n)) – grids with horizontal and vertical directions
V_y (numpy.ndarray, size=(m,n)) – grids with horizontal and vertical directions
geoTransform –
x ({float, numpy.ndarray}) – map locations of the seeds
y ({float, numpy.ndarray}) – map locations of the seeds
- Returns
C – grid with the specific catchments
- Return type
numpy.ndarray, size=(m,n), dtype=boolean
See also
- dhdt.postprocessing.terrain_tools.d8_flow(Class, V_x, V_y, iter=20, analysis=True)¶
re-organize labelled dataset by steering towards convergent zones
- Parameters
Class (numpy.ndarray, size=(m,n)) – labelled array
V_x (numpy.ndarray, size=(m,n)) – grids with horizontal and vertical directions
V_y (numpy.ndarray, size=(m,n)) – grids with horizontal and vertical directions
iter (integer, {x ∈ ℕ | x ≥ 0}, default=20) – number of iterations to be executed
analysis (boolean) – keeps intermediate results, to see the evolution of the algorithm
- Returns
Class_new – newly evolved labelled array
- Return type
numpy.ndarray, size=(m,n)
See also
- dhdt.postprocessing.terrain_tools.get_d8_dir(V_x, V_y)¶
get directional data of direct neighbors, following [GM97].
- Parameters
V_x (numpy.ndarray, size=(m,n), float) – array of displacements in horizontal direction
V_y (numpy.ndarray, size=(m,n), float) – array of displacements in vertical direction
- Returns
d8_dir – compass directions starting from North going in clockwise direction
- Return type
numpy.ndarray, size=(m,n), range=0…7
Notes
It is important to know what type of coordinate systems is used, thus:
N 0 NW | NE 7 | 1 W----+----E 6----8----2 SW | SE 5 | 3 S 4References
- GM97
Garbrecht & Martz, “The assignment of drainage direction over flat surfaces in raster digital elevation models”, Journal of hydrology, vol.193 pp.204-213, 1997.
- dhdt.postprocessing.terrain_tools.get_d8_offsets(shape, order='C')¶
get the offsets of the direct neighbors for linear indexing
- Parameters
shape (list) – dimensions of the array of interest
order ({'C','F'}) – index order in either C-like or Fortran like order
- Returns
d8_offsets – relative offsets of neighbors
- Return type
numpy.ndarray, size=(8)
Notes
Two different orders are used:
0 1 2 3 0 4 8 4 5 6 7 1 5 9 8 9 ... 2 6 ... order 'C' order 'F'