main classes
general sentinel-2 product class
- class s2d2.sentinel2_product.Sentinel2Product(path: str | PathLike | None = None)[source]
Notes
For Sentinel-2 Level 1C the metadata is scattered over the files and folders. In order to make this organization more understandable, the following visualization is made:
* 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
- load_metadata() None[source]
Load meta-data from the product file (MTD_MSIL1C.xml)
Notes
The metadata structure of the Sentinel2 product is as follows:
* MTD_MSIL1C.xml └ n1:Level-1C_Tile_ID ├ n1:General_Info │ ├ Product_Info │ │ ├ PRODUCT_URI <- .SAFE folder name │ │ ├ PROCESSING_LEVEL │ │ └ PROCESSING_BASELINE │ ├ SPACECRAFT_NAME │ ├ SENSING_ORBIT_NUMBER │ └ SENSING_ORBIT_DIRECTION ├ 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
generic functions
advanced processing of flight paths functions
- s2d2.orbit_tools.acquisition_angles(p_x: ndarray, g_x: ndarray) tuple[ndarray, ndarray][source]
given satellite and ground coordinates, estimate observation angles
- s2d2.orbit_tools.calculate_correct_mapping(grid: Sentinel2Anglegrid, inclination: float = 98.5621, revolutions_per_day: float = 14.30824258387262, radius: float | None = None, mean_altitude: float | None = None)[source]
- Parameters:
grid (Sentinel2Anglegrid) –
- with the following entries:
- zenith{numpy.ndarray, numpy.masked.array}, size=(k,l,h)
observation angles of the different detectors/bands
- azimuth{numpy.ndarray, numpy.masked.array}, size=(k,l,h)
observation angles of the different detectors/bands
- bandnumpy.ndarray, size=(h,)
number of the band, corresponding to the third dimension of ‘zenith’
- detectornumpy.ndarray, size=(h,)
number of the detector, corresponding to the third dimension of ‘zenith’
- geotransformtuple, size=(6,)
geotransform of the grid of ‘zenit’ and ‘azimuth’
- crsosgeo.osr.SpatialReference() object
coordinate reference system (CRS)
inclination (float, unit=degrees) – angle of the orbital plane in relation to the equatorial plane
revolutions_per_day (float) – amount of revolutions a satellite platform performs around the Earth
mean_altitude (float, unit=meter) – estimate of mean satellite altitude above the Earth surface
- Returns:
l_time (numpy.ndarray, size=(p,), unit=seconds) – asd
lat, lon (float, unit=degrees) – ground location of the satelite at time 0
radius (float, unit=metre) – distance away from Earths’ center
inclination (unit=degrees) – tilt of the orbital plane in relation to the equator
period (float, unit=seconds) – time it takes to complete one revolution
time_para (numpy.ndarray, size=(p,b)) – polynomial fitting parameters for the different bands (b)
combos (numpy.ndarray, size=(h,2)) – combinations of band and detector pairs, corresponding to the third dimension of ‘grid.zenith’
- s2d2.orbit_tools.earth_axes() tuple[float, float][source]
get axis length of the WGS84 ellipsoid
- Returns:
major_axis (float, unit=meter) – the largest axis of the ellipsoid
minor_axis (float, unit=meter) – the shortest axis of the ellipsoid
See also
- s2d2.orbit_tools.earth_eccentricity() float[source]
get the eccentricity of the WGS84 ellipsoid
- Returns:
eccentricity – amount of deviation from circularity
- Return type:
See also
- s2d2.orbit_tools.ground_vec(lat_arr: ndarray, lon_arr: ndarray, eccentricity: float | None = None, major_axis: float | None = None) ndarray[source]
get ground coordinates in Cartesian system
- Parameters:
lat_arr (numpy.ndarray, size=(m,n), unit=degrees) – location of observation angles
lon_arr (numpy.ndarray, size=(m,n), unit=degrees) – location of observation angles
eccentricity (float) – eccentricity squared, if None default is WGS84 (6.69E-10)
major_axis (float, unit=meters) – length of the major axis, if None default is WGS84 (6.37E6)
- Returns:
g_x – ground location in Cartesian coordinates
- Return type:
numpy.ndarray, size=(m*n,3)
- s2d2.orbit_tools.line_of_sight(lat_arr: ndarray, lon_arr: ndarray, zn_arr: ndarray, az_arr: ndarray, eccentricity: float | None = None, major_axis: float | None = None) tuple[ndarray, ndarray][source]
- Parameters:
lat_arr (numpy.ndarray, size=(m,n), unit=degrees) – location of observation angles
lon_arr (numpy.ndarray, size=(m,n), unit=degrees) – location of observation angles
zn_arr (numpy.ndarray, size=(m,n), unit=degrees) – polar angles of line of sight to the satellite, also known as, delcination and right ascension.
az_arr (numpy.ndarray, size=(m,n), unit=degrees) – polar angles of line of sight to the satellite, also known as, delcination and right ascension.
eccentricity (float) – eccentricity squared, if None default is WGS84 (6.69E-10)
major_axis (float, unit=meters) – length of the major axis, if None default is WGS84 (6.37E6)
- Returns:
sat (numpy.ndarray, size=(m*n,3)) – observation vector from ground location towards the satellite
g_x (numpy.ndarray, size=(m*n,3)) – ground location in Cartesian coordinates
- s2d2.orbit_tools.observation_calculation(ltime, sat, g_x, radius, inclination, period, omega_0, lon_0)[source]
- Parameters:
ltime (numpy.array, size=(m,1)) –
sat (numpy.ndarray, size=(m,3)) – observation vector from ground location towards the satellite
g_x (numpy.ndarray, size=(m,3)) – ground location in Cartesian coordinates
radius (float, unit=meter) – radius towards the orbiting satellite
inclination (float, unit=degrees, range=-180...+180) – inclination of the orbital plane with the equator
period (float, unit=seconds) – time it takes to revolve one time around the Earth
omega_0 (float, unit=degrees) – angle towards ascending node, one of the orbit Euler angles
lon_0 (float, unit=degrees) – ephemeris longitude
- Returns:
v_x
- Return type:
numpy.array, size=(m,3)
- s2d2.orbit_tools.orbital_fitting(sat, g_x, inclination, lat=None, lon=None, radius=None, period=None, revolutions_per_day=None, mean_altitude=None, gps=None, convtol=0.001, orbtol=1.0, maxiter=20, printing=False)[source]
- Parameters:
sat (numpy.ndarray, size=(m*n,3)) – observation vector from ground location towards the satellite
g_x (numpy.ndarray, size=(m*n,3)) – ground location in Cartesian coordinates
lat (float, unit=degrees) – location of satellite within orbital track
lon (float, unit=degrees) – location of satellite within orbital track
radius (float, unit=meter) – radius towards the orbiting satellite
inclination (float, unit=degrees, range=-180...+180) – inclination of the orbital plane with the equatorial plane
period (float, unit=seconds) – time it takes to revolve one time around the Earth
- s2d2.orbit_tools.partial_obs(ltime, sat, g_x, lat, lon, radius, inclination, period)[source]
numerical differentiation, via pertubation of the observation vector
- s2d2.orbit_tools.standard_gravity() float[source]
provide value for the standard gravitational parameter of the Earth
- Returns:
mu – standard gravity
- Return type:
float, unit= m**3 * s**-2
- s2d2.orbit_tools.wgs84_param() tuple[float, float][source]
get paramters of the WGS84 ellipsoid
- Returns:
major_axis (float, unit=meter) – largest axis of the ellipsoid
flattening (float) – amount of compression of the ellipsoid
See also
Notes
The flattening (f) stem from the major- (a) en minor-axis (b) via:
\[f = (a - b)/a\]or via the eccentricity (e):
\[f = 1 - \sqrt{ 1 - e^2}\]
basic reading of geo-referenced imagery
- s2d2.mapping_input.read_geo_image(fname, boi=None, no_dat=None)[source]
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
read_geo_infobasic function to get meta data of geographic imagery
Examples
>>> import os >>> from s2d2.mapping_input import read_geo_image >>> fpath = os.path.join(os.getcwd(), "data.jp2" ) >>> I, spatialRef, geotransform, targetPrj = read_geo_image(fpath)
- s2d2.mapping_input.read_geo_info(fname)[source]
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
coordinate transformations
- s2d2.mapping_tools.ecef2llh(xyz)[source]
transform 3D cartesian Earth Centered Earth fixed coordinates, to spherical angles and height above the ellipsoid
- Parameters:
xyz (np.array, size=(m,3), unit=meter) – np.array with 3D coordinates, in WGS84. In the following form: [[x, y, z], [x, y, z], … ]
- Returns:
llh – np.array with angles and height. In the following form: [[lat, lon, height], [lat, lon, height], … ]
- Return type:
np.array, size=(m,2), unit=(deg,deg,meter)
- s2d2.mapping_tools.ecef2map(xyz, spatial_ref)[source]
transform 3D cartesian Earth Centered Earth fixed coordinates, to map coordinates (that is 2D) in a projection frame
- Parameters:
xyz (np.array, size=(m,3), float) – np.array with 3D coordinates, in WGS84. In the following form: [[x, y, z], [x, y, z], … ]
spatial_ref (osgeo.osr.SpatialReference) – target projection
- Returns:
xyz – np.array with planar coordinates, within a given projection frame
- Return type:
np.array, size=(m,2), float
- s2d2.mapping_tools.ll2map(ll, spatial_ref)[source]
transforms angles to map coordinates (that is 2D) in a projection frame
- Parameters:
llh (np.array, size=(m,2), unit=(deg,deg)) – np.array with spherical coordinates. In the following form: [[lat, lon], [lat, lon], … ]
spatial_ref (osgeo.osr.SpatialReference) – target projection system
- Returns:
xyz – np.array with 2D coordinates. In the following form: [[x, y], [x, y], … ]
- Return type:
np.array, size=(m,2), unit=meter
Examples
Get the Universal Transverse Mercator (UTM) cooridnates from spherical coordinates:
>>> import numpy as np >>> from osgeo import osr >>> proj = osr.SpatialReference() >>> proj.SetWellKnownGeogCS('WGS84') >>> lat, lon = 52.09006426183974, 5.173794246145571# Utrecht University >>> proj.SetUTM(32, lat>0)
>>> xy = ll2map(np.array([[lat, lon]]), proj) >>> xy array([[ 237904.03625329, 5777964.65056734, 0. ]])
- s2d2.mapping_tools.map2ll(xy, spatial_ref)[source]
transforms map coordinates (that is 2D) in a projection frame to angles
- Parameters:
xy (np.array, size=(m,2), unit=meter) – np.array with 2D coordinates. In the following form: [[x, y], [x, y], … ]
spatial_ref (osgeo.osr.SpatialReference) – target projection system
- Returns:
ll – np.array with spherical coordinates. In the following form: [[lat, lon], [lat, lon], … ]
- Return type:
np.array, size=(m,2), unit=(deg,deg)
transfer between mapping and image domain
- s2d2.image_coordinate_tools.create_local_crs()[source]
create spatial refence of local horizontal datum
- Returns:
crs – the coordinate reference system via GDAL SpatialReference description
- Return type:
{osr.SpatialReference, string}
See also
- s2d2.image_coordinate_tools.get_bbox(geotransform, rows=None, cols=None)[source]
given array meta data, calculate the bounding box
- Parameters:
geotransform (tuple, size=(6,)) – georeference transform of an image.
rows (integer, {x ∈ ℕ | x ≥ 0}) – amount of rows in an image.
cols (integer, {x ∈ ℕ | x ≥ 0}) – amount of collumns in an image.
- Returns:
bbox – bounding box, in the following order: min max X, min max Y
- Return type:
np.ndarray, size=(4,), dtype=float
See also
Notes
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | i | x --------+--------> --------+--------> | | | | image | j map | based v based |
- s2d2.image_coordinate_tools.map2pix(geotransform, x, y)[source]
transform map coordinates to local image coordinates
- Parameters:
geotransform (tuple, size=(6,1)) – georeference transform of an image.
x (np.array, size=(m), ndim={1,2,3}, dtype=float) – horizontal map coordinate.
y (np.array, size=(m), ndim={1,2,3}, dtype=float) – vertical map coordinate.
- Returns:
i (np.ndarray, ndim={1,2,3}, dtype=float) – row coordinate(s) in local image space
j (np.ndarray, ndim={1,2,3}, dtype=float) – column coordinate(s) in local image space
See also
Notes
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | i | x --------+--------> --------+--------> | | | | image | j map | based v based |
- s2d2.image_coordinate_tools.pix2map(geotransform, i, j)[source]
transform local image coordinates to map coordinates
- Parameters:
geotransform (tuple, size={(6,1), (8,1)}) – georeference transform of an image.
i (np.array, ndim={1,2,3}, dtype=float) – row coordinate(s) in local image space
j (np.array, ndim={1,2,3}, dtype=float) – column coordinate(s) in local image space
- Returns:
x (np.array, size=(m), ndim={1,2,3}, dtype=float) – horizontal map coordinate.
y (np.array, size=(m), ndim={1,2,3}, dtype=float) – vertical map coordinate.
See also
Notes
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | i | x --------+--------> --------+--------> | | | | image | j map | based v based |
- s2d2.image_coordinate_tools.pix_centers(geotransform, rows=None, cols=None, make_grid=True)[source]
provide the pixel coordinate from the axis, or the whole grid
- Parameters:
geotransform (tuple, size={(6,), (8,)}) – georeference transform of an image.
rows (integer, {x ∈ ℕ | x ≥ 0}, default=None) – amount of rows in an image.
cols (integer, {x ∈ ℕ | x ≥ 0}, default=None) – amount of collumns in an image.
make_grid (bool, optional) – Should a grid be made. The default is True.
- Returns:
X (np.ndarray, dtype=float) –
“make_grid” : size=(m,n)
otherwise : size=(1,n)
Y (np.ndarray, dtype=float) –
“make_grid” : size=(m,n)
otherwise : size=(m,1)
Notes
Two different coordinate system are used here:
indexing | indexing ^ y system 'ij'| system 'xy' | | | | j | x --------+--------> --------+--------> | | | | image | i map | based v based |
metric conversions
- s2d2.unit_conversion.datenum2datetime(t)[source]
some write a date as a number “YYYYMMDD” conver this to numpy.datetime
- Parameters:
t ({int,numpy.array}) – time, in format YYYYMMD
- Returns:
dt – times in numpy time format
- Return type:
{numpy.datetime64, numpy.array}, type=numpy.datetime64[D]
See also
- s2d2.unit_conversion.datetime2calender(dt)[source]
Convert array of datetime64 to a calendar year, month, day.
- Parameters:
dt (np.datetime64) – times of interest
- Returns:
cal – calendar array with last axis representing year, month, day
- Return type:
numpy array
- s2d2.unit_conversion.datetime2datenum(dt)[source]
convert numpy.datetime to a date as a number “YYYYMMDD”
- Parameters:
dt ({numpy.datetime64, numpy.array}, type=numpy.datetime64[D]) – times in numpy time format
- Returns:
t – time, in format YYYYMMD
- Return type:
{int, numpy.array}
See also
- s2d2.unit_conversion.datetime2doy(dt)[source]
Convert array of datetime64 to a day of year.
- Parameters:
dt (np.datetime64) – times of interest
- Returns:
year (integer, {x ∈ ℕ}) – calender year
doy (integer, {x ∈ ℕ}) – day of year
- s2d2.unit_conversion.deg2arg(θ)[source]
adjust angle to be in bounds of an argument angle
- Parameters:
θ (unit=degrees) –
- Returns:
θ
- Return type:
unit=degrees, range=-180…+180
See also
Notes
The angle is declared in the following coordinate frame:
^ North & y | - <--|--> + | +----> East & x
- s2d2.unit_conversion.deg2compass(θ)[source]
adjust angle to be in bounds of a positive argument angle,like a compass
- Parameters:
θ (unit=degrees) –
- Returns:
θ
- Return type:
unit=degrees, range=0…+360
See also
Notes
The angle is declared in the following coordinate frame:
^ North & y | - <--|--> + | +----> East & x
- s2d2.unit_conversion.deg2dms(ang)[source]
convert decimal degrees to degree minutes seconds format
- Parameters:
ang ({float,np.array}, unit=decimal degrees) – angle(s) of interest
- Returns:
deg ({integer,np.array}) – degrees
min ({integer,np.array}, range=0…60) – angular minutes
sec ({float,np.array}, range=0…60) – angular seconds
- s2d2.unit_conversion.dms2deg(deg, min, sec)[source]
convert degree minutes seconds format to decimal degrees
- Parameters:
deg ({integer,np.array}) – degrees
min ({integer,np.array}, range=0...60) – angular minutes
sec ({float,np.array}, range=0...60) – angular seconds
- Returns:
ang – angle(s) of interest
- Return type:
{float,np.array}, unit=decimal degrees
handling formats and specifics
sentinel-2 tiling structure
- s2d2.handler.mgrs.get_bbox_from_tile_code(tile_code, tile_path=None)[source]
Get the bounds of a certain MGRS tile
- Parameters:
tile_code (string, e.g.: '05VMG') – MGRS tile coding
tile_path (string) – Path to the geometric metadata
- Returns:
bounding box, in the following order: min max X, min max Y
- Return type:
numpy.ndarray, size=(1,4), dtype=float
- s2d2.handler.mgrs.get_geom_from_tile_code(tile_code, tile_path=None)[source]
Get the geometry of a certain MGRS tile
- Parameters:
tile_code (string) – MGRS tile coding, e.g.: ‘05VMG’
tile_path (string) – Path to the geometric metadata
- Returns:
Geometry of the MGRS tile, in lat/lon
- Return type:
shapely.geometry.polygon.Polygon
See also
get_tile_code_from_geom,get_bbox_from_tile_code
- s2d2.handler.mgrs.get_mgrs_tile(lat, lon)[source]
return a military grid reference system zone designation string. This zoning is used in the tiling of Sentinel-2.
- Parameters:
- Returns:
mgrs_code
- Return type:
string
Notes
- The MGRS 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
“CC” 100 km square identifier
- The following acronyms are used:
CRS : coordinate reference system
MGRS : US military grid reference system
WKT : well known text
UTM: universal transverse mercator projection
- s2d2.handler.mgrs.get_tile_codes_from_geom(geom, tile_path=None)[source]
Get the codes of the MGRS tiles intersecting a given geometry
- Parameters:
geom ({shapely.geometry, string, dict, GeoDataFrame, GeoSeries}) – geometry object with the given dict-like geojson geometry, GeoSeries, GeoDataFrame, shapely geometry or well known text, i.e.: ‘POLYGON ((x y, x y, x y))’
tile_path (string) – Path to the geometric metadata
- Returns:
MGRS tile codes
- Return type: