-
Notifications
You must be signed in to change notification settings - Fork 3
Sample dataset analysis
Below are details about the datasets to store in seis-lab-data catalog. These were initially provided by the client and then further analyzed.
The information provided by the client is implicitly defining some of the required properties that a survey-related record instance should have. It seems a possible dataset schema could include
class Dataset:
acronym: str
domain_type: DatasetDomain
category: DatasetCategory
stage: DatasetWorkflowStage
primary_files: list[Asset]
secondary_files: list[Asset]
class Asset:
media_type: str
storage_url: URLAccording to the client each dataset belongs to a certain domain:
- D1. Geophysical
- D2. Geotechnical
It will likely be beneficial if it is possible to add additional domain types, either at runtime or as an initialization parameter
According to the client each dataset has a certain category:
- C1. Bathymetry
- C2. Backscatter
- C3. Seismic
- C4. Magnetometer/gradiometer
- C5. Superficial sediment samples
- C6. Cores
- C7. CPT tests
It will likely be beneficial if it is possible to add additional domain types, either at runtime or as an initialization parameter
The client classifies each dataset according to the following stage:
- S1. Raw data
- S2. Quality control data
- S3. Processed data
- S4. Interpreted data
It will likely be beneficial if it is possible to add additional domain types, either at runtime or as an initialization parameter
- Acronym: MBES/BATHY (aka mbes-bathy)
- Method type: D1. Geophysical
- Data type: C1. Bathymetry
- Purpose of this data: Multibeam echo sounder - bathymetry data
- Acronym: MBES/BSK (aka mbes-bsk)
- Method type: D1. Geophysical
- Data type: C2. Backscatter
- Purpose of this data: Multibeam echo sounder - backscatter data
GDAL can read it natively since it should be a GEOTiff file
from osgeo import gdal, osr
# Open the file
ds = gdal.Open("data_example/02_MBES_Example/MBES_BCKS/FIG-24-M001_Backscatter_RGB_4m.tif")
# most metadata should be available off the shelf
# Bands
print(f"Number of bands: {ds.RasterCount}")
for i in range(1, ds.RasterCount + 1):
band = ds.GetRasterBand(i)
print(f"Band {i}: dtype={gdal.GetDataTypeName(band.DataType)}")
# GeoTransform
gt = ds.GetGeoTransform()
print("GeoTransform:", gt)
# gt gives: (originX, pixelWidth, 0, originY, 0, pixelHeight)
# Projection
proj_wkt = ds.GetProjection()
srs = osr.SpatialReference(wkt=proj_wkt)
print("Projection (WKT):", proj_wkt)
# Try to identify EPSG code
srs.AutoIdentifyEPSG()
print("Authority:", srs.GetAuthorityName(None))
print("EPSG Code:", srs.GetAuthorityCode(None))
# Metadata
print("Metadata:", ds.GetMetadata())- Acronym: SBP
- Method type: D1. Geophysical
- Data type: C3. Seismic
- Purpose of this data: Sub-bottom Profiling
This format is used by:
| dataset | stage | primary output | secondary output |
|---|---|---|---|
| mbes-bathy | S3. Processed data | x | - |
This seems to refer to a directory containing project files related to the CARIS HIPS and SIPS software, which is a product for processing hydrographic survey data.
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S1. Raw data | x | - |
| mbes-bathy | S2. Quality control data | x | - |
| mbes-bathy | S3. Processed data | - | x |
| mbes-bathy | S4. Interpreted data | x | x |
| mbes-bsk | S1. Raw data | x | - |
| mbes-bsk | S2. Quality control data | - | x |
| mbes-bsk | S4. Interpreted data | x | x |
Python's package osgeo provides a gdal module that can read .xyz data into da Dataset object. From a Dataset the Bands (layers or channels) can be read and region of a layer be extracted. Elements of a Band are the points of the dataset.
+--------------------------------------------------------------+
| Dataset |
| (whole raster file, e.g., example.xyz, example.tif) |
| |
| ┌───────────────────┐ ┌───────────────────┐ ┌───────────┐|
| | Band 1 | | Band 2 | | Band 3 ||
| | (full layer data) | | (full layer data) | | ... ||
| | e.g., Red channel | | e.g., Green chan. | | ... ||
| └───────────────────┘ └───────────────────┘ └───────────┘|
| ↑ |
| | **Each Band covers the *whole* extent** |
| | with same XSize × YSize as dataset |
| |
+--------------------------------------------------------------+
Gdal offers GDal.Grid() or even the command line gdal_grid tool to generated raster images by interpolating pixels values. Sparse points data are typically stored in .csv files
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S4. Interpreted data | x | - |
| mbes-bsk | S4. Interpreted data | x | - |
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S3. Processed data | x | - |
| mbes-bsk | S3. Processed data | x | - |
This seems to be a CARIS Spatial Archive file, a raster format for gridded bathymetry and elevation, as per https://docs.safe.com/fme/2017.1/html/FME_Desktop_Documentation/FME_ReadersWriters/csar/csar.htm
The only utility I can could come up with is stored at the ocvniciusnog repository on GitHub. This repo ships a series of python script that serve as a batch utility to run a carisbatch command that, to my understanding, is in the repo as a compiled Python file. As such this is a binary os/architecture dependent file
The requirements section is the README.md file explicitly refers to a Teledyne Caris Hips and Sips (with a valid license) product. So it's not clear to me how even how this stuff can be technically and legally run
Specific details about TCH&S can be found at the (hideously formatted) URL https://www.teledynecaris.com/docs/7.0/hips%20and%20sips//index.html#page/CARIS%2520HIPS%2520and%2520SIPS%2520Help%2FHIPS%2520NOTICES.01.1.html%23
The utility seems to be quite flexible and able to manipulate files in other formats, such as the .kmall format
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S3. Processed data | - | x |
| mbes-bathy | S4. Interpreted data | x | x |
| mbes-bsk | S4. Interpreted data | x | x |
A technical description can be found at lp360. that quotes ArcView™ 3.2 Help System. The explanation refers to a header file with extension .hdr which seems necessary in order to interpret correctly the data in .flt
The Binary Raster File format is a simple format that can be used to transfer raster data between various applications. It consists of two files, the IEEE floating-point file and a supporting ASCII header file. The header file must have the same name as the data file, but with a .hdr file extension. The header data includes the following keywords and values:
- ncols - number of columns in the data set.
- nrows - number of rows in the data set.
- xllcenter or xllcorner - x-coordinate of the center or lower-left corner of the lower-left cell.
- yllcenter or yllcorner - y-coordinate of the center or lower-left corner of the lower-left cell.
- cellsize - cell size for the data set.
- nodata_value - value in the file assigned to cells whose value is unknown. This keyword and value is optional. The nodata_value defaults to -9999.
- byteorder - the byte order of the binary cell values. You can choose between two keywords, msbfirst or lsbfirst. Msbfirst is used for cell values written with the most significant bit first. Lsbfirst is used for cell values written with the least significant bit first.
Header file example, ncols 480 nrows 450 xllcorner 378923 yllcorner 4072345 cellsize 30 nodata_value -32768 byteorder msbfirst The data file will be a matrix of 32 bit signed IEEE floating-point values. The file will have a line of binary numbers for each row in the data set. The first line of data is the top row of the data set, moving from left to right. The Binary Raster File format can be imported using the Spatial or 3D Analyst extensions for ArcView™ 3.x or ArcGIS™. References: ArcView™ 3.2 Help System
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S1. Raw data | - | x |
| mbes-bathy | S2. Quality control data | - | x |
| mbes-bathy | S3. Processed data | - | x |
| mbes-bathy | S4. Interpreted data | x | x |
| mbes-bsk | S1. Raw data | - | x |
| mbes-bsk | S2. Quality control data | - | x |
| mbes-bsk | S3. Processed data | - | x |
| mbes-bsk | S4. Interpreted data | x | x |
This format belongs to ESRI Arc Geodatabase Format family and to my understanding is actually a set of database files stored in a directory. The format should be readable with gdal's OpenFileGDB driver. This is in contrast with what stated in the fileformat.com which claims "File GDB is ESRI proprietary format and its specifications are not available publicly. Due to this reason, the file format details for FIle GDB could not be documented anywhere other than some sources who have done it partially by reverse engineering". The reverse engineering project is available from GitHub. But I guess if gdal can squeeze out of a gdb directory all we need to know we can get away with it
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S3. Processed data | x | x |
| mbes-bathy | S4. Interpreted data | x | - |
| mbes-bsk | S3. Quality control data | x | - |
| mbes-bsk | S3. Processed data | x | - |
| mbes-bsk | S4. Interpreted data | x | - |
GeoTiff is clearly a TIFF image format adopting an extended (and I guess standard) set of specialized IFDs tags and it's supported by GDAL for creation, copy and georeferencing operations.
Gdal standard operations should work on these files. For example this should get the spatial extension of a map
from osgeo import gdal
# Open the dataset
ds = gdal.Open("your_raster.tif")
# Get GeoTransform
# gt = (originX, pixelWidth, rotationX, originY, rotationY, pixelHeight)
gt = ds.GetGeoTransform()
# Raster size
width = ds.RasterXSize
height = ds.RasterYSize
# Calculate bounds
minx = gt[0]
maxy = gt[3]
maxx = minx + width * gt[1] + height * gt[2]
miny = maxy + width * gt[4] + height * gt[5]
extent = (minx, miny, maxx, maxy)
print(f"Spatial extent (minX, minY, maxX, maxY): {extent}")
# Close dataset
ds = None
but I guess this is a good option if the geometry returned is expressing the extent in coordinates that have been correcly localized
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S1. Raw data | x | - |
| mbes-bsk | S1. Raw data | x | - |
This seems to be data from a Kongsberg product (the Kongsberg multibeam echosounder systems - MBES). The format is binary and based on datagrams and is documented at https://www.kongsbergdiscovery.online/sis/kmall/html/index.html, where the C/C++ structures defining the datagrams are also available. A python module exists and it's available from Github (last update 2 years ago, last tagged release 5 years ago). The author claims the package is incomplete, but it seems to me it's able to pull some of the basic information from a .kmall file. Spatial and temporal ranges determination can be inferred for example from datagrams of the type '#SPO' and #MRZ
the #SPO datagram C/C++ struct definition is
struct EMdgmSPOdataBlock_def
{
uint32_t timeFromSensor_sec;
uint32_t timeFromSensor_nanosec;
float posFixQuality_m;
double correctedLat_deg;
double correctedLong_deg;
float speedOverGround_mPerSec;
float courseOverGround_deg;
float ellipsoidHeightReRefPoint_m;
int8_t posDataFromSensor[MAX_SPO_DATALENGTH];
};
the Python package produces reports that are instances of pandoc objects, but provides also support for direct datagram reading. Notice: the following example works with the last tagged release of the module. Documentation on GitHub consistently refers to a capitalized KMALL module, which probably means they changed the name of the module/class in the unreleased code in the repository
>>> import kmall
>>> km = kmall.kmall('data/0007_20190513_154724_ASVBEN.kmall')
>>> km.report_packet_types()
Count Size: Min Size Max Size
MessageType
b'#CPO' 87 13572 156 156
b'#IIP' 1 1054 1054 1054
b'#IOP' 1 1332 1332 1332
b'#MRZ' 472 32556126 68688 69120
b'#SCL' 87 6612 76 76
b'#SKM' 86 1159680 13368 13632
b'#SPO' 87 13572 156 156
b'#SVP' 1 292 292 292
b'#SVT' 2181 148308 68 68
>>> km.OpenFiletoRead()
>>> SPOIndex["ByteOffset"].iloc[0]
16862
>>> SPOIndex["ByteOffset"].iloc[1]
446474
>>> km.FID.seek(SPOIndex["ByteOffset"].iloc[0],0)
16862
>>> dg = km.read_EMdgmSPO()
>>> dg
{'header': {'numBytesDgm': 156, 'dgmType': b'#SPO', 'dgmVersion': 0, 'systemID': 40, 'echoSounderID': 2040, 'dgtime': 1557762444.784905, 'dgdatetime': datetime.datetime(2019, 5, 13, 15, 47, 24, 784905)},
'cmnPart': {'numBytesCmnPart': 8, 'sensorSystem': 0, 'sensorStatus': 577, 'padding': 16087},
'sensorData': {'timeFromSensor_sec': 1557762444, 'timeFromSensor_nanosec': 778000000, 'datetime': datetime.datetime(2019, 5, 13, 15, 47, 24, 778000), 'posFixQuality_m': 1.7999999523162842, 'correctedLat_deg': 45.465634228499994, 'correctedLong_deg': -83.61456595566666, 'speedOverGround_mPerSec': 2.569525957107544, 'courseOverGround_deg': 119.93270111083984, 'ellipsoidHeightReRefPoint_m': 138.61500549316406, 'posDataFromSensor': b'$INGGK,154724.778,051319,4527.93805371,N,08336.87395734,W,4,31,1.8,EHT138.615,M*7A'}}
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S1. Raw data | - | x |
| mbes-bathy | S2. Quality control data | - | x |
| mbes-bathy | S3. Processed data | - | x |
| mbes-bathy | S4. Interpreted data | - | x |
.shp files are a type of vector data file format used in geographic information systems (GIS) to store the location, shape, and attributes of geographic features like points, lines, and polygons. They are part of a collection of files that work together to represent spatial data, with the .shp file specifically containing the geometry data.
GDAL can read ESRI Shapefiles and GDAL Manual Page stresses out that
Normally the OGR Shapefile driver treats a whole directory of shapefiles as a dataset, and a single shapefile within that directory as a layer. In this case the directory name should be used as the dataset name. However, it is also possible to use one of the files (.shp, .shx or .dbf) in a shapefile set as the dataset name, and then it will be treated as a dataset with one layer.
ChatGPT produced an example easy to understand
from osgeo import ogr
# Open the Shapefile
shapefile = "your_file.shp"
ds = ogr.Open(shapefile)
if ds is None:
raise RuntimeError(f"Could not open {shapefile}")
# Get the first layer (a Shapefile usually has one layer)
layer = ds.GetLayer(0)
# Loop through features
for feature in layer:
# Geometry (spatial data)
geom = feature.GetGeometryRef()
print("Geometry type:", geom.GetGeometryName())
print("WKT:", geom.ExportToWkt()) # WKT representation
# Attribute fields
for i in range(feature.GetFieldCount()):
field_name = feature.GetDefnRef().GetFieldDefn(i).GetName()
field_value = feature.GetField(i)
print(f"{field_name}: {field_value}")
print("-" * 40)
I also found an (to me) interesting introduction
# Open the shapefile
ds = ogr.Open(folder + "EU_vector/gadm36.shp")
# An .shp file should be made by one single layer
# Get the first layer in the file
lyr = ds.GetLayer()
# And from a layer the spatial object reference
# Spatial reference
sr = lyr.GetSpatialRef()
# A spatial reference is not yet what we need, it
# has to be used to get
wkt = sr.ExportToWkt()
'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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
Also QGis is quoted as able to read .shp. Is it possibile to devise a command line utility program based on QGis libraris that can extract the metadata we need?
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| mbes-bathy | S1. Raw data | - | x |
| mbes-bathy | S2. Quality control data | - | x |
| mbes-bathy | S3. Processed data | - | x |
This format is used as:
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| sbp | S1. Raw data | x | - |
This seems to be a format used by EdgeTech - the specification can be found in this pdf document:
https://www.edgetech.com/wp-content/uploads/2023/04/0023492_Rev_N.pdf
The documentation refers to a jsfdefs.h C header that actually is not available from the EdgeTech website, but an apparently EdgeTech produced JSFDefs.h file is available from AVP-MOOS. Within this repository you'll find leitura_sonar.c, an example of reading data from a .jsf file. In principle this might be a base from which we could build a Python extension (comments in the file are written in Portuguese) or a C/C++ command line program that can print the metadata of interest
is a standard file format used in the geophysical industry for recording digital seismic data. In essence, this model serves as a blueprint for what a SEG-Y file should look like. The latest standard revision is available from lib.seg.org
This format is used as
| dataset | stage | primary file | secondary file |
|---|---|---|---|
| S1. Raw Data |
The segy file structure was defined long ago (some sources quote 1973 other 1975) and has an apparently simple structure. The text header obviously ships metadata crucial to a correct interpretation of the following sections but to my understanding there is no standard and any interpretation is application specific. My overall impression is that this is a format specification that helps to define other application/manufacturer specific data formats
SEG-Y Revision 0:
┌──────────┐ ┌────────┐ ┌──────────────┐ ┌──────────────┐
│ Txt File │ │ Binary │ │ Trace 1 │ │ Trace N │
│ Header │───►│ Header │───►│ Header 240B │─ ... ─►│ Header 240B │
│ 3200B | | 400B | | Data | | Data |
└──────────┘ └────────┘ └──────────────┘ └──────────────┘
The most complete information I could find is located at readthedocs.io and refers to a Python segy module that works as an helper to define SEG-Y data formats
In summary the key components are
-
Quoting from the aforementioned source "A Text File Header that stores the information required to parse the textual file header of the SEG-Y file". It sounds quite a circular specification. I guess the information in this section is instrumental to properly understanding the textual information in the following (Trace?) sections
-
Application specific binary data
-
Extended Text Header (not shown in the diagram above)
-
Traces: traces are individual seismic traces
Segyio is an actively maintained repository for reading/writing SEG-Y datasets. The package documentation is as well located on readthedocs.io.
Segpy is a piece of software created by Sixty
North, but it requires additional license if we want to redistribute it as part
of our project. It can be installed calling pip install segpy
The P1/11 (P111) Geophysical position data exchange format. A QGIS plugin exists to read such file format, I don't know how to harness this stuff for the IPMA project though.