search

alisonpchase / pigments-from-rrs

MIT License
1 stars 1 forks source link

Define the IOPs of water #9

Open cisaacstern opened 2 years ago

cisaacstern commented 2 years ago

@cisaacstern I'm diving back in, and looking at the next step in the data processing for our next incremental PR. A key step in the inversion of Rrs to get information (in our case, pigments) is defining the absorption and scattering spectra of water itself, so that we can account for and remove that signal from the measured light reflectance spectra.

This is a standard activity in optical oceanography computing, so I did a search and found some other GitHub repos where people have worked on this. In particular, there is this one: https://github.com/tadz-io/hydropt which is a full inversion model. In their case, they define a basis vector for all phytoplankton absorption, and where my method differs is that I separate out the absorption by different pigments. Anywho, they do have some code that defines the "IOPs" (Inherent Optical Properties, i.e. the absorption and scattering spectra) for "clear_nat_water". My matlab code uses the temperature and salinity of seawater to adjust the IOP spectra, so it may be different, but in case it's useful we could look at how they did the inversion of Rrs.

So I think our next PR should be to define the IOPs of water, which are known, but vary as a function of temperature and salinity, so maybe we can pass the temperature and salinity info in somewhere.

Originally posted by @alisonpchase in https://github.com/alisonpchase/pigments-from-rrs/issues/8#issuecomment-999001962

cisaacstern commented 2 years ago

@alisonpchase, made your earlier comment ☝️ into this new Issue for organizing our conversation.

If I understand correctly, the function get_water_iops from your earlier Matlab code represents the approach you would like to now replicate in Python. I've copied that function below (click its name in the list below to expand the code).

I hope its okay to copy sections of the Matlab code here (if not, let me know here or offline, and I can remove). This was the easiest way I could think to structure a conversation around next steps.

The hierarchical bullets below get_water_iops represent what I believe are all the other functions called from within this function, all of which we will need to replace with Pythonic versions in order to get our IOPs routine to work. If the sub-function is a link, that means its a builtin Matlab function, and the link goes to what I believe is the corresponding place in the Matlab documentation. If it's an expandable button, then it's a custom function, and the relevant code is copied within the expanded content.

My first question for you is: does the below organization seem accurate? Does the list below accurately reflect the full extent of what we will need to translate into Python, in order to implement your IOPs methods? If not, please let me know what I've overlooked or misunderstood.

get_water_iops ```matlab function [a_sw,bb_sw] = get_water_iops(wave,T,S) %function to obtain seawater absorption and backscattering spectra %pure water absorption from Mason et al 2016 for 250-550, Pope and Frye for %550-730 nm, and Smith and Baker for 730-800 nm %salt water backscattering from Zhang et al 2009 %corrected for in situ temperature and salinity conditions Sullivan et al. 2006 wl_water1 = [250 260 270 280 290 300 302 304 306 308 310 312 314 316 318 320 322 324 326 328 330 ... 332 334 336 338 340 342 344 346 348 350 352 354 356 358 360 362 364 366 368 370 372 ... 374 376 378 380 382 384 386 388 390 392 394 396 398 400 402 404 406 408 410 412 414 ... 416 418 420 422 424 426 428 430 432 434 436 438 440 442 444 446 448 450 452 454 456 ... 458 460 462 464 466 468 470 472 474 476 478 480 482 484 486 488 490 492 494 496 498 ... 500 502 504 506 508 510 512 514 516 518 520 522 524 526 528 530 532 534 536 538 540 ... 542 544 546 548 550 ]; wl_water2 = [552.5 555 557.5 560 562.5 565 567.5 570 572.5 575 577.5 580 582.5 585 587.5 590 592.5 595 597.5 ... 600 602.5 605 607.5 610 612.5 615 617.5 620 622.5 625 627.5 630 632.5 635 637.5 640 642.5 ... 645 647.5 650 652.5 655 657.5 660 662.5 665 667.5 670 672.5 675 677.5 680 682.5 685 687.5 ... 690 692.5 695 697.5 700 702.5 705 707.5 710 712.5 715 717.5 720 722.5 725 727.5 730 732.5 ... 735.0 737.5 740.0 742.5 745.0 747.5 750.0 752.5 755.0 757.5 760.0 762.5 765.0 767.5 770.0 772.5 ... 775.0 777.5 780.0 782.5 785.0 787.5 790.0 792.5 795.0 797.5 800.0]; wl_water = [wl_water1 wl_water2]; aw1 = [58.7100 51.5000 43.5700 22.3000 9.3900 4.6700 4.3300 3.5600 3.1300 2.7500 2.3600 2.0500 ... 1.8500 1.7600 1.6300 1.4700 1.3300 1.2400 1.1800 1.1200 1.0700 1.0100 0.9900 0.9500 ... 0.9100 0.8500 0.8200 0.8100 0.8200 0.8400 0.8900 0.9400 0.9700 0.9800 0.9900 1.0600 ... 1.1500 1.2000 1.2100 1.2200 1.2400 1.2700 1.2900 1.3300 1.3700 1.4300 1.4700 1.5100 ... 1.5500 1.6200 1.7000 1.7500 1.8500 1.9600 2.0800 2.2200 2.3700 2.4800 2.5700 2.5900 ... 2.6600 2.7100 2.8000 2.8800 3.0000 3.1200 3.2200 3.3100 3.4400 3.5800 3.7600 3.9500 ... 4.1700 4.4200 4.8000 5.2200 5.7400 6.2600 6.9100 7.5100 8.0800 8.4200 8.6300 8.7700 ... 8.9300 9.0900 9.3300 9.5500 9.7900 9.9900 10.3000 10.6500 11.0000 11.3800 11.7700 12.1400 ... 12.5400 12.9400 13.3600 13.9100 14.6000 15.4500 16.4800 17.7400 19.2600 20.7300 22.4200 24.2400 ... 26.6800 29.7100 33.0000 35.6900 37.3800 38.2100 38.7800 39.1700 39.6200 40.1700 40.8800 41.6200 ... 42.4200 43.3000 44.3600 45.4100 46.4500 47.5400 48.8200 50.4000 52.2400 54.2500 56.2900]; aw2 = [0.0593 0.0596 0.0606 0.0619 0.064... 0.0642 0.0672 0.0695 0.0733 0.0772 0.0836 0.0896 0.0989 0.11 0.122 0.1351 0.1516 0.1672 0.1925 0.2224 0.247 0.2577... 0.2629 0.2644 0.2665 0.2678 0.2707 0.2755 0.281 0.2834 0.2904 0.2916 0.2995 0.3012 0.3077 0.3108 0.322 0.325 0.335... 0.34 0.358 0.371 0.393 0.41 0.424 0.429 0.436 0.439 0.448 0.448 0.461 0.465 0.478 0.486 0.502 0.516 0.538 0.559... 0.592 0.624 0.663 0.704 0.756 0.827 0.914 1.007 1.119 1.231 1.356 1.489 1.678 1.7845 1.9333 2.0822 2.2311 2.3800... 2.4025 2.4250 2.4475 2.4700 2.4900 2.5100 2.5300 2.5500 2.5400 2.5300 2.5200 2.5100 2.4725 2.4350 2.3975 2.3600... 2.3100 2.2600 2.2100 2.1600 2.1375 2.1150 2.0925 2.0700]; a_water = [aw1*1e-3 aw2]; % 1e-3 comes from Mason et al. 2016, Table 2 a_pw = interp1(wl_water,a_water,wave,'linear','extrap'); %temp and salinity correction for water absorption (need to know at what T it was measured): if S==0 || isnan(S)==1,S=35;end if T==0 || isnan(T)==1,T=22;end T_pope=22.0; % use salt water scattering from Zhang et al 2009 [betasw124, bb_sw, beta90sw, theta] = betasw124_ZHH2009(wave, S, T); % use the temp and salinity corrections from Sullivan et al. 2006 [psiT,psiS] = tempsal_corr(wave); %temperature and salinity corrections: a_sw = (a_pw+psiT*(T-T_pope)+psiS*S); end ```
cisaacstern commented 2 years ago

Let's make sure the bulleted list above is indeed an accurate reflection of the work ahead. If so, I'll preview my suggestion that the next thing we work on are function definitions and docstrings for each of these custom functions. An example would be,

# Matlab function definition line:
#    `function [a_sw,bb_sw] = get_water_iops(wave,T,S)`

# Python definition and docstring

def get_water_iops(wave, T, S):
    """Description of what this function does.

    Parameters
    ----------
    wave: (input type for wave?)
        Description of wave.
    T: (input type for T?)
        Description of T.
    S: (input type for S?)
        Description of S.

    Returns
    -------
    a_sw: (return type for a_sw?)
        Description of a_sw.
    bb_sw: (return type for bb_sw?)
        Description of bb_sw.
    """
    pass   # this function doesn't do anything yet, so we're just going to `pass` for now

These definitions could go into a new module, e.g. pigments_from_rrs/iops.py.

cisaacstern commented 2 years ago

@alisonpchase, one more design note for now. The repo you found provides a good example for how to manage reference vectors in a Pythonic way. That is, store the data itself as CSVs, and then assign it to a capitalized constant with pd.read_csv, as in, e.g.: https://github.com/tadz-io/hydropt/blob/master/hydropt/bio_optics.py

This is more idiomatic for Python than defining the vectors within the function bodies, as you've done in the Matlab functions. It's also more flexible and maintainable: this way we'll be able to swap out those reference vectors without touching the functions themselves.

In terms of specific implementation, though, it looks like importlib.resources is a more contemporary way of managing this than pkg_resources (the approach used in tadz-io/hydropt/).

alisonpchase commented 2 years ago

@cisaacstern delayed big thanks for this framework to get going with the water iops functions. I will work on the function definitions in a new module, maybe pigments_from_rrs/water_iops? Should I create a PR for this? Or just work on the code and commit it to repo? I think my understanding is that the PR isn't necessary since it's my own repo, but can be a nice way for us to keep track of something we are collaborating on?

It should be fine to have the copied matlab code posted as long as we are referencing the papers. Currently the references are mentioned in the comments of the higher level functions - we can make sure we have complete references in our docstrings maybe, as we get the functions organized.

cisaacstern commented 2 years ago

maybe pigments_from_rrs/water_iops? Should I create a PR for this? Or just work on the code and commit it to repo? I think my understanding is that the PR isn't necessary since it's my own repo, but can be a nice way for us to keep track of something we are collaborating on?

Yep, pigments_from_rrs/water_iops sounds like a perfect spot for it.

Good question re: PRs, to which the answer is: definitely make a PR for this. 😄 While it is strictly true that you can commit changes directly to this repo, you are 100% correct that in our collaborative context, the PR makes it infinitely easier to keep track of what the new material is. Even if you were working completely on your own, using PRs is still a good way to compartmentalize changes to your code.

Good luck with this! And definitely let me know how I can help as you go along (no question too small!).

alisonpchase commented 2 years ago

Thanks!! sounds good. So to check my workflow: I checkout a new branch, work on my .py module of functions, commit it, then open a PR on it to have a decided space for us to collaborate on it?

cisaacstern commented 2 years ago

Yep! Let's say your new branch is called iops.

Once you've committed some code to your local branch, you can

git push origin iops

To push that branch to a new branch of the same name on GitHub.

Then GitHub should prompt you to open a PR from the new branch. Once the PR is open, all commits pushed to the associated branch will be tracked in the PR.