geomagpy / magpy

MagPy (or GeomagPy) is a Python package for analysing and displaying geomagnetic data.
BSD 3-Clause "New" or "Revised" License
49 stars 27 forks source link

wrong conversions in IMF values when D in hundredths of minute #122

Closed stephanbracke closed 1 year ago

stephanbracke commented 2 years ago

When using IMF files with a declination in the format, It will not be correctly interpreted and cause missing interpretation of values being xyz ibstead of hdz. The cause of this is that in the code format_imf there is no interpretation of the file values with the corresponding headers everithing is always considered as being nanoTesla values also the DataComponents headers is always XYZF which causes problems later on you can see this conversion in the code for readIMF

https://github.com/geomagpy/magpy/blob/e7761aff7f2c33e46298ef7cbbce987b96e009c7/magpy/lib/format_imf.py#L1750-L1775

The same is true for the write imf. The way I solved is to extract the method that adds the headers into a method that is inspired by the code you find in the iaga2002 format

#components = XYZF or HDZF etc

def add_col_info_headers(header, components):

    header["col-x"] = components[0].upper()
    header["col-y"] = components[1].upper()
    header["col-z"] = components[2].upper()
    header["unit-col-x"] = 'nT'
    header["unit-col-y"] = 'nT'
    header["unit-col-z"] = 'nT'
    header["unit-col-f"] = 'nT'
    if components.endswith('g'):
        header["unit-col-df"] = 'nT'
        header["col-df"] = 'G'
        header["col-f"] = 'F'
    else:
        header["col-f"] = 'F'
    # print ("VAR", varstr)
    if components in ['dhzf', 'dhzg']:
        header["unit-col-y"] = 'deg'
        header['DataComponents'] = 'HDZF'
    elif components in ['ehzf', 'ehzg']:
        # consider the different order in the file
        header["col-x"] = 'H'
        header["col-y"] = 'E'
        header["col-z"] = 'Z'
        header['DataComponents'] = 'HEZF'
    elif components in ['dhif', 'dhig']:
        header["col-x"] = 'I'
        header["col-y"] = 'D'
        header["col-z"] = 'F'
        header["unit-col-x"] = 'deg'
        header["unit-col-y"] = 'deg'
        header['DataComponents'] = 'IDFF'
    elif components in ['hdzf', 'hdzg']:
        header["unit-col-y"] = 'deg'
        header['DataComponents'] = 'HDZF'
    else:
        header['DataComponents'] = 'XYZF'

I also added two methods proprietary at the imfv format that encodes and decodes the values according to IMFV standards based on the units

def decode(type, value)-> float:
    """
    type : variable is deg or nT
    value the value to decode
    """

    if type == 'deg':
        return value/6000
    else:
        return value/10
def encode(type, value)-> float:
    """
    type : variable is deg or nT
    value the value to encorde
    """
    if type == 'deg':
        return value*6000
    else:
        return value*10

These methods make it now simple to adapt the readIMF and writeIMF methods

# in READIMF
....
                headers['DataComponents'] = block[4]
                headers['DataSensorAzimuth'] = float(block[8])/10/60
                headers['DataSamplingRate'] = '60 sec'
                headers['DataType'] = block[5]
                format = block[4].lower()
                add_col_info_headers(headers, format)
....
               if not int(data[0+index]) > 999990:
                      array[1].append(decode(headers["unit-col-x"], float(data[0 + index])))
                  else:
                      array[1].append(np.nan)
                  if not int(data[1+index]) > 999990:
                      array[2].append(decode(headers["unit-col-y"], float(data[1 + index])))
                  else:
                      array[2].append(np.nan)
                  if not int(data[2+index]) > 999990:
                      array[3].append(decode(headers["unit-col-z"], float(data[2 + index])))
                  else:
                      array[3].append(np.nan)
                  if not int(data[3+index]) > 999990:
                      array[4].append(decode(headers["unit-col-f"],float(data[3+index])))
........

# in writeIMF

elemtype = datastream.header['DataComponents']
       if not isnan(elemx):
            x = encode(header['unit-col-x'],elemx)
        else:
            x = 999999
        if not isnan(elemy):
            y =encode(header['unit-col-y'],elemy)
        else:
            y = 999999
        if not isnan(elemz):
            z = encode(header['unit-col-z'],elemz)
        else:
            z = 999999
        if not isnan(elemf):
            f = encode(header['unit-col-f'],elemf)

The method add_col_info_headers could be externalised and used in all conversions making sure that all headers are always set with the correct names and formats needed for the GUI

leonro commented 1 year ago

corrected with 20dda573a713d838668f89f5172415552b0ceac6