twVolc / PyCamPermanent

Permanent PiCam (SO2) installation project software
GNU General Public License v3.0
2 stars 2 forks source link

Save ICA masses in emission rate file #233

Open twVolc opened 3 weeks ago

twVolc commented 3 weeks ago

There's times where flow_glob fails to find a plume speed (gets inf as it things there is no lag between lines). It's then impossible to pull useful data from the sequence (technically you could look at the flow_nadeau emission rate, then divide by it's plume speed to get the ICA mass, but this is convoluted and won't work when the plume speed is calculated to be 0, which it sometimes is).

I think it would therefore be useful for the line's ICA mass to be saved in the emission rate file, so that a user could manually generate a plume speed and then get to an emission rate themselves by scaling the ICA mass.

ubdbra001 commented 1 week ago

This is from the fluxcalc.py in pyplis:

def det_emission_rate(cds, velo, pix_dists, cds_err=None, velo_err=None,
                      pix_dists_err=None, mmol=MOL_MASS_SO2):
    """Determine emission rate.

    :param cds: column density in units cm-2 (float or ndarray)
    :param velo: effective plume velocity in units of m/s (float or ndarray)
        Effective means, that it is with respect to the direction of the normal
        vector of the plume cross section used (e.g. by performing a scalar
        product of 2D velocity vectors with normal vector of the PCS)
    :param pix_dists: pixel to pixel distances in units of m (float or ndarray)

    """
    if cds_err is None:
        logger.info("Uncertainty in column densities unspecified, assuming 20 % of "
              "mean CD")
        cds_err = mean(cds) * 0.20
    if velo_err is None:
        logger.info("Uncertainty in plume velocity unspecified, assuming 20 % of "
              "mean velocity")
        velo_err = mean(velo) * 0.20

    if pix_dists_err is None:
        logger.info("Uncertainty in pixel distance unspecified, assuming 10 % of "
              "mean pixel distance")
        pix_dists_err = mean(pix_dists) * 0.10

    C = 100**2 * mmol / N_A
    phi = sum(cds * velo * pix_dists) * C
    dphi1 = sum(velo * pix_dists * cds_err)**2
    dphi2 = sum(cds * pix_dists * velo_err)**2
    dphi3 = sum(cds * velo * pix_dists_err)**2
    phi_err = C * sqrt(dphi1 + dphi2 + dphi3)
    return phi, phi_err

Would the ICA mass be: cds * pix_dists * C?

twVolc commented 1 week ago

Would the ICA mass be: cds * pix_dists * C?

Yes it would be the sum of that vector. So sum(cds * pix_dists * C)

ubdbra001 commented 1 week ago

Okay great, that should be easy to access and add as we're passing all those values (except C) to pyplis from (e.g.) here: https://github.com/twVolc/PyCamPermanent/blob/67c8b0f76b004e48397f3b51cef5f5b4b759c528/pycam/so2_camera_processor.py#L3311-L3312

Looks like there are multiple places self.det_emission_rate_kgs() is run so we'd maybe have to make changes around them all?

twVolc commented 1 week ago

The ICA will be the same for each flow-type, so I think for each image sequence we should only have to do this once, i.e. once for each line in the above for-loop rather than having to do this separately for each flow. So hopefully we'd only have to add this calculation to one place? Then perhaps it could be added to the res dictionary with a key of 'ICA' or something like that? Then hopefully this can be saved in the emission rates file

ubdbra001 commented 1 week ago

Should it appear in all the output files or just the flow_glob one? Also, I'm assuming it should be different for each ICA line?

twVolc commented 1 week ago

Good question. It's only important that it's saved somewhere, but I guess duplicating it wouldn't be bad. I suppose if you didn't select flow_glob you'd miss out on ICA being saved if it were only saved in flow_glob, so perhaps we should duplicate it in all files.