Augment existing HYPSO NetCDF files with GCP as metadata

[sivertba]

---

name: Feature request about: Suggest an idea for this project title: 'Augment existing H1 NetCDF files with GCP as metadata' labels: 'enhancement' assignees: ''

---

Is your feature request related to a problem? Please describe. The existing NetCDF files for HYPSO-1 lack GCP (Ground Control Points) data, which is essential for accurate georeferencing for the HYPSO-1 system. instead of everyone recreating it every time we should aim to add it to the existing NetCDF structure, as additional metadata.

Describe the solution you'd like a Python script, module or functionality within this library to augment existing NetCDF files for HYPSO-1 with GCP as they are generated by QGIS, and include it as metadata. The metadata should include CRS type as well as the GCP data, some of which is included in the output of QGIS.

Describe alternatives you've considered Manually adding GCP data to the NetCDF files each is time-consuming and prone to errors, but is what is currently being done.

[sivertba]

@DennisNTNU edit as you see fit :)

[sivertba]

There should be one attribute with GCP, and a function in the package that generate lat lon from this field. Each EPSG4326 shall be used.

[DennisNTNU]

An example ground control point file from QGIS looks as follows:

#CRS: GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["unknown"],AREA["World"],BBOX[-90,-180,90,180]],ID["EPSG",4326]]
mapX,mapY,sourceX,sourceY,enable,dX,dY,residual
-123.78811622954093252,49.60202492727613333,71.2374193548386,-244.703870967742,1,-0.68054821121393161,0.05132248897356817,0.68248067054014361
-123.53404392959831171,49.71177170954961611,24.7480645161288,-103.077096774194,1,0.02281588512232702,-0.19742732986620126,0.19874132733785663
-123.97589197069532929,49.69682406327249424,43.5596774193547,-378.852258064516,1,2.20129444895331972,-0.25568247795295918,2.21609358568740644
-124.09554647713179065,49.76901651708252672,20.7004838709675,-457.876451612903,1,-1.36606890483406573,0.43126349280026943,1.4325265976507493
-124.16941820568965227,49.44625230574087027,136.885322580645,-449.472903225807,1,-0.06528727290262282,-0.00403313903467506,0.06541172840962518
-124.51066025119462211,49.366955591260961,174.58564516129,-632.577741935484,1,-0.8775470859813268,-0.12411940793720078,0.88628128466136502
-123.91745626326216723,49.17331080719448977,225.21252688172,-248.14752688172,1,-1.70085575183418314,0.77908533266048607,1.8707977560693525
-123.70903772370824925,49.01965574018832683,274.040483870967,-85.5247311827957,1,-1.76565158881521711,-0.39372114423332505,1.80901682482560111
-124.18898278649413669,48.85669744064466613,346.434354838709,-359.873602150538,1,-2.58262875342182951,-1.09175986219139531,2.8039099262803302
-124.45851875416626342,48.92015217175199382,330.513870967741,-533.919569892473,1,0.50640769989584555,0.33166744590073449,0.60535283363020398
-123.75747457566673404,48.61533190408727023,416.322580645161,-45.8198924731183,1,2.72611147658517439,0.88062431913709815,2.86481813981709443
-124.40403651124780993,48.55973779267353763,455.64193548387,-442.714086021505,1,0.74450057919645474,0.79583200785606323,1.08978424339502622
-123.8728899360034319,48.50488848758529059,460.370537634408,-101.612258064516,1,2.35844712541700119,-0.71952259166948807,2.46576268187157321
-123.91408337745231449,48.06868337145770198,613.844516129031,-40.5516129032252,1,-3.89277605159020368,-0.07963489032147208,3.89359051565399339
-124.73750986074801972,48.16938180652307011,605.312473118279,-573.341720430106,1,-0.86175772277302087,0.14764172081720517,0.87431370256036145
-124.48824807007541438,47.67597871244522878,769.819856630824,-338.55634408602,1,-1.31857022341216634,0.1876081003639456,1.33184992900528343
-124.33597056974001305,47.40486270780001377,859.04652329749,-195.327670250895,1,0.2353700479253007,0.33099645527454413,0.40614986503096567

Its a CSV file, with the first line defining the target coordinate reference system (CRS). The second line are the names of the fields of the CSV file. Each following line represents one ground control point (GCP).

The target CRS shall be EPSG:4326, such that the target CRS coordinates are (longitude,latutude).

The source hypso image coordinates are (l,S-s). l refers to line index, corresponding to the frames/lines dimension. S-s corresponds to the 'samples', i.e. swath dimension. Capital S represents the 'samples' dimension size, i.e. 684 for a nominal cube, and s represents the samples index. Example: given a points in an hypso image, if QGIS assigns coordiantes (212.37, -272.44) to the point, then the second coordinate is actually 684-272.44 = 411.56. This coordinate transformation needs to be taken into account, and converted when reading the GCPs and storing them in the netCDF file.

The remaining four fields are not relevant to be added to the netCDF file, except that points that do not have 1, but 0 in the 'enabled' field should be skipped.

Code should be added to the hypso-package that takes a path to such a GCP file and a corresponging hypso image netCDF fle as input, and adds the CRS and GCP info to the netCDF file.

By fitting a mapping function, e.g. a 2D second degree polynomial, to the GCPs, the images can be georeferenced with a high degree of accuracy. Code for doing that should also be added to the hypso-package.

Here is an example python code snippet that generate lat,lon values for every spatial pixel in an hypso image, given the GCPs:

import numpy as np
import pandas as pd
import sklearn.preprocessing as skp #import PolynomialFeatures
import sklearn.linear_model as sklm #import LinearRegression
import matplotlib.pyplot as plt

def generate_lat_lon_from_GCPs(gcps, cube_dimensions):

    input_points = gcps[:,[2,3]]
    input_points[:,1] = -gcps[:,3]
    input_data_lat = gcps[:,1]
    input_data_lon = gcps[:,0]

    poly = skp.PolynomialFeatures(degree=2)
    input_features = poly.fit_transform(input_points)

    model_lat = sklm.LinearRegression(fit_intercept=False)
    model_lat.fit(input_features, input_data_lat)

    model_lon = sklm.LinearRegression(fit_intercept=False)
    model_lon.fit(input_features, input_data_lon)

    data_lines = cube_dimensions[0]
    data_samples = cube_dimensions[0]

    points = np.zeros([data_lines,data_samples,2])
    for i in range(data_lines):
        for j in range(data_samples):
            points[i,j,0] = j+0.5
            points[i,j,1] = i+0.5
    points.shape = (data_lines*data_samples,2)

    lats = model_lat.predict(poly.transform(points))
    lons = model_lon.predict(poly.transform(points))

    lats.shape = (data_lines,data_samples)
    lons.shape = (data_lines,data_samples)
    points.shape = (data_lines,data_samples,2)

    return lats,lons

gcps_pd = pd.read_csv('gcp.points', skiprows=1)
gcps = gcps_pd.values[:,0:4]

lats, lons = generate_lat_lon_from_GCPs(gcps, [956,684])