Python conversion stuff

[mgcooper]

Starting a general list of ideas related to the python conversion.

Indents are currently three spaces. This is a holdover from a time when all my code was written on an eleven inch laptop. Plus three is way better than four. But we’ll need four for python. Use an autoformatter for this. Probably best to format the m files first, because formatting the py files may not be possible until they’re debugged.

Can try matlab to python conversion utilities such as smop and matlab2python. I tried both. They won’t work out of the box, but they might accelerate the process. A simple find and replace, keeping track of idiosyncrasies, might be faster. Start with one function, delete all ‘end’ statements, replace all ‘function’ statements with ‘def’ and so on. Patterns will emerge and soon you’ll have a list of conversions that can be applied to other functions.

The input met files are in a high level file format known as a timetable. I will convert them to flat binary files.

Another option is to first get it running in octave. This will reveal idiosyncratic issues like the met files, because octave is closer to python in terms of not being compatible with matlabic stuff like timetables.

[altanshapovalov]

Hi Matt,
Currently reinstalling MatLab to match your version (2021B) and I was curious if you'd recommend installing any additional toolboxes that are crucial for this model to run successfully

[mgcooper]

Hi Max. I will reply with more soon, but there are no required toolboxes. The interp1 issue is probably related to the roundn issue. Would I be correct that you replaced roundn with round, and then began receiving the interp1 error? If so, take the negative sign off of the second argument to round (formerly roundn, which uses negative numbers, whereas for round we need them to be positive). See if that solves the problems. Also make sure all paths in the repo including the 'util' folder are on your path.

[mgcooper]

Hi Max,
I replaced a few calls to custom functions (the ones in the util/ folder) with built-in functions. I also replaced the instances of roundn with round. If you have local commits or unstaged changes, you can do a pull with rebase/autostash and /or merge and deal with any conflicts. If you are unsure about this we can discuss tomorrow and walk through it. If you do not have any local commits or unstaged changes, you can pull my changes and the model should run.

There should be no toolbox dependencies or any issues with the version. But no harm in updating to R2021b. A few tools you can use to work through this stuff:

1. you can get a list of all required toolboxes and functions using something like this:
   [filelist,productlist] = matlab.codetools.requiredFilesAndProducts('icemodel.m')
   filelist will be a list of all the files required to run the code, productlist would list the toolboxes required, if any.

2. In any folder on your matlab path, you can click the little drop-down arrow in the upper right of your file bowser and navigate to 'Reports' then 'Dependency Report' and a gui will appear with information similar to the function call above.

Sometimes there are spurious results with both of these. I ran this stuff on my end and it should be good to go. Either way, we can sort it out tomorrow, but I wanted to let you know about the changes I pushed to python-feature which should take care of the errors you received.

[altanshapovalov]

Hi Matt,

Sorry for the long response, but to address your questions:

1. Yes, I replaced “roundn” with “round” in CVSPECTRAL.m and CVTHERMAL.m (both in the Functions folder). To answer your question more directly, then yes, I believe that the error with interp1 started showing up after I messed with the ‘roundn’ command. “Roundn” also seems to be used in a number of other functions but I have not touched those because no error on their behalf has been queued. Should I leave those alone if they won’t make a fuss or is it better to change it in every document for consistency?

2. for the paths, I always just right click on my “icemodel” folder and choose that all folders and subfolders are added to the path, so I believe that it’s all included.

3. The terminology is a bit new to me but after some googling, I think You’re suggesting that I ‘update’ my ‘util’ folder with whatever changes you ‘committed’ recently. I’m not too sure how to pull your updated content on GitHub Desktop. The only “pull origin” I see will just override the documents that I corrected, and I’m not sure if that will match whatever updates you created. In fact, the last updates I see to the files in your GitHub profile are 3 days old. Am I looking in the wrong place?
   

4. Lastly, I ran the tests that you recommended and didn’t get any toolbox recommendations. However, as I tried to run icemodel_run, I received an error saying that 'roundn' requires Mapping Toolbox. I hope that after I pull your changes, this problem will resolve itself. Otherwise, let me know if I should get that before the meeting so that I won’t waste our time.

See you tomorrow!

[altanshapovalov]

Hey Matt,

I opened up the code again and found a bunch of <<<< ==== and >>>>> lines in some .m files.

There is an explanation on this forum (here) and I think I understand the solution only halfway.

During our meeting you suggested that I create a branch of my own; should I create it going off the "main" branch of "python-feature"? Here is what GitHub Desktop is telling me for each

I have a fear of making your repository super messy so could you please direct me about 1) where I should create my branch, 2) how I can copy my slightly tweaked version of your code into there w/o those merge errors, and 3) how can i ensure that those merge errors don't show up again?
Whenever I improvise with organization and try to arrange things myself, I generally dig myself deeper into a hole so I'd like some guidance from you about this)

[altanshapovalov]

Hey Matt, quick UPD:

Johnny helped me create a branch of my own that I stemmed from python-feature (named it python-feature-max)
I'll make sure that the model can run again; I'll be tinkering with SEBFLUX.m: try to figure out how it works and which functions it uses

[mgcooper]

Hi Max,That sounds great. SEBFLUX is a logical place to start. You could write your own script that calls SEBFLUX to test it out. You could use the saved PHYSCONS.mat (see near the top of icemodel.m) or type them in yourself to get aquatinted with the physical constants that are used. You could also run icemodel and save the output, then load ice1.m and use the values as inputs to SEBFLUX, or you could make up dummy values, or whatever you want. Any questions feel free to ask. Have fun!

[altanshapovalov]

Hey Matt,

So I spent a lot of time staring at the code and got a slightly better understanding of it (at least in terms of structure).
My task for the week was to play with SEBFLUXm; this function has a number of outputs, but I thought that "balance" could be a good start for visualizing what this function yielded.

So, I extracted balance (along with Time, albedo, melt, freeze, drain, and surface temperature). For the first figure, I just wanted to compare Balance with Albedo. (I'm gonna number my questions to you for organization's sake)

I think a pattern can somewhat be seen..

1. is it even relevant looking at balance? what would it really tell us?
2. Is balance unitless?
3. Any idea as to why the albedo spikes so high? A 0.8 is an understandable constant for fresh snow, but then it gets even higher.. What's more curious to me is when it peaks: around april (for a while) and mid september. Is that just caused by snowstorm events? I guess my question is more about whether it makes sense to have an albedo of 0.96 mid-Sept.

For a second visual, I decided to compare melt, runoff, drain, and freeze all in one figure

4. Silly question, but again, what would be units for all of these? I'm speculating that they'd be the same.. (I'm guessing that if I was a good scholar, I would've read the relevant literature and knew the answer to this question, but perhaps it's easier to just confirm with you))
5. I'm really curious as to why all 4 variables just plateau after a certain point. No more data? Or does it truly become stable? Also, do they all start at zero because we just wanna track the change? Asking this because by the end of the year, the values don't plateau back on zero.

Lastly, decided to compare surface temperature, albedo, and melt (since in the previous figure, it was the highest of the 4)

What's pleasant to me in this figure is that some sort of correlation can be observed for all three subplots.
6. Was the comparison of the variables that I pulled in any way useful? Should I try to visualize some other ones? I totally ignored all the different Q's and other variables in ice1.

7. Will it prove useful to perform similar juxtaposition of data with a different station? As in, BEHAR vs the other one (KANM?)
8. If i was to play with some variables and replace them with my own, which should I mess with? Albedo? If so, then how would I change it, because there are 8,700+ data points for it. Perhaps multiply by 0.90 to see how the overall balance would be affected if the overall albedo was lowered by 10%? (but then would that even be realistic..)
9. If this analysis isn't that useful to spend time on, which other direction could I move in order to i) better understand the code and/or ii) prepare for the transfer from matlab to python?

[altanshapovalov]

Also, on a completely different note:
My brain can't help but come up with silly ideas and I had this especially pervasive one, so I decided to implement it and show it to you, even though I recognize that it's probably super unprofessional, especially so early in the project.
Buuuut, based on what you've written in Matlab, it's quite clear that you're a coding wiz so this is only fair

I hope you dig it)

[altanshapovalov]

Hey Matt,

Just a minor UPD (to show what I've been doing and where I'm at)
After many hours, I i) installed python 3.11, ii) installed VS Code (as advised by Johnny), iii) imported pip and PIL, iv) overlapped 2 images.
Here's the result of that:

Some problems I ran into:

• I can make the head by itself larger, but when I try to paste it on the other, I can only control its position but not size, so that head is a bit too small and i'm not sure how to fix it within Python (i.e., I didn't sit in front of the code for long enough to figure it out by myself);
• Overlaying one image on top of another was piece of cake. I had the matlab logo and your cutout head (fully opaque), and was able to overlap that, but getting your head transparent in Python is very difficult. I can make JUST the head image transparent, but overlaying that python-created-transparent head onto the logo just did not work. Something about converting to RGBA is not connecting in my brain. After a loooot of forum digging, I just thought "why not just make the image itself opaque and then overlay it" and that worked, but that's kinda cheating...
• I wanna tackle the letters next but I don't think I can get the letters to look as nice as the original... I'll probably crop the letters out or overlay them with a textbox, but it will clearly look foreign. Will that be okay or should I try to get as close as I can to that original-looking Mattlab?
• Lastly, I'm just curious if you have any ideas about making VS Code interface look more like MatLab (I really miss the "comment" selections button); perhaps this is a question more for Johnny

This task is fun but python is really iffy about its formatting; how it reports errors is bananas! So messy
I now see why people like MatLab so much)

[mgcooper]

WOW bravo! I am amazed you got your python and vscode workflow running this quickly! I don't think you should worry about those details about the image. The idea was just to get python installed and running in an editor and figuring out the basics of importing packages and stuff. You've managed to do that way faster than I expected. Now you can decide whether you want to keep developing stuff in python or matlab. I think a little of both is good at this stage, just enough matlab that you can efficiently run the model and understand the syntax to translate to python, which is a better tool for you to master in the long run.

As for vscode, I have a question. My guess is that you are running in interactive mode, meaning you are writing commands in a script and sending those commands to the interpreter one by one or a few at a time, and immediately seeing the results (the alternative would be that you are running a .py file from the terminal command line all at once, which I assume you are not). But do you know if you are running in an iPython interactive session? If so, you should have a second window where you see your output, similar to a jupyter notebook. This is a screenshot with a green box around the iPython interactive session in my vscode:

This is how you want things set up. You can highlight the code on the left and hit ctrl-enter to send it to the interactive session on the right. You can insert section breaks using # %% and then ctrl-enter (or cmd-enter) will send the entire section to the interpreter. The orange box is around the terminal. The red circle is the 'run' button that would send the file to the terminal. I would avoid that.

Either way, the python error / debugging stack trace is annoying compared with matlab. The main thing to understand is that the errors will be listed in a certain order (the stack) so you can trace the error. Imagine you enter an invalid parameter as input to a function call. That input gets sent to some function, but that function sends the input to another function, which sends it to a third function, and in this third function the parameter actually gets used and the error occurs. The error stack will show three error messages, listed in the order the functions were called. So you normally just look at the last (or first) message which shows the true source of the error. But sometimes it is helpful to see the entire stack especially when you are developing code so they show it all. You'll get the hang of it over time.

Couple other questions

• did you install python via homebrew? Or did you install anaconda?
• are you using an environment manager such as pyenv or conda?
• did you install the python vscode extension?

My guess is you installed anaconda and when you launch vscode you are using the anaconda python interpreter which has an iPython backend, in which case I should have just asked this up front and avoided a lot of the details above. But if not, it would be good to get that set up to make your life easier. Over time you will want to set your vscode environment settings to make it work exactly how you want it. You should be able to press cmd-shift-p and then select the "Preferences: Open User Settings (JSON)" whcih will open the settings.json file where you set all sorts of preferences.

In vscode, you should be able to "comment" selections by pressing cmd+/ . It should work if you just have the cursor on a line or if multiple lines are highlighted. I don't know if there is a button.

[altanshapovalov]

I want to tinker with your Matlab code a little more (just because I feel like there's still a lot for me to truly comprehend).
What would you say is the best approach: read/run functions line by line, OR make my own versions of functions and produce my own variables and figures? You gave me some notes/tasks earlier so I'll work with plotting Qh/Qsi/Qli etc
Also, I think I'm taking an intro-to-python course next quarter with Johnny so that should give me a solid boost with Python in the next term.

As for running VS Code: I am actually writing things in a .py file, but I do hit that "Run" button that you recommended straying from. I don't think I have those specific hotkeys set up for those functions but I'll dig around; the section advice was really helpful, I like it in matlab.

Addressing your questions:

1. did you install python via homebrew? Or did you install anaconda? – hmm.. neither? To get the 3.11 python, I just downloaded it from their website (https://www.python.org/downloads/). Then I just downloaded VS Code from their respective website. I've read about homebrew to use that as an alternative for installing pip (brew install pip), but I didn't actually use it. Frankly, clueless about both of those! Should I have them?
2. are you using an environment manager such as pyenv or conda? – I don't think I do. All I'm using is VS Code (and the default/recommended interpreter python 3.11.0 64 bit, but I think that's irrelevant). Should I go educate myself on adopting one of these?
3. did you install the python vscode extension? – I believe so... When I first opened VS Code, i think it asked me for things to install, and I just instinctively chose Python. So I think I have that covered. Is there a way to double check?
   Fyi, here is the recommended interpreter that I am using
   

On a very different note (and I should also address this to Johnny & Sarah): I'm progressively getting bombarded with the notion that writing skills are increeedibly important (for many reasons we all know), and I believe that! But I'm curious if you (and Johnny and Sarah) have any tips on me improving my writing skills, whether that's for scientific literature, grant proposals, or basic applications for scholarships. Would you recommend practice writing daily? Or perhaps reading certain materials? (scientific papers are not the easiest to read so I'm not sure if that is the kind of language I want to adopt)
For leisure, I like reading books in every other language but english, but I think that if I want to improve, there's no way around it)

[JohnnyRyan1]

Couple of quick comments, Max - you’re actually due to TA my Python course next quarter since there was no else qualified. I’ll make sure you’re well prepared.

In regards to package managers, we recently decided to migrate to venv so you can ignore the conda stuff for now.

I agree that more tinkering with the model code in matlab is a good idea. It will be a great learning experience converting it to Python but only worth starting that when you are a little more familiar with the model.

Lastly writing of course is one the more important skills to have in this world. You’ll get plenty of practice next quarter. I’d like you to submit a FINSST proposal which is due in Feb 11. We should discuss some ideas next time we all meet.

[altanshapovalov]

Hey Matt,

I'm a bit stuck (with downloading data) and want to get clarification (on variables)

My current task is to download the reanalysis data (specifically, M2T1NXFLX from here). Johnny gave me really good details on the directions (attached), but there are 2 points that I need help clearing up. The former is more directed for you, Matt.
merra2-notes.txt

Preamble: The penultimate step before downloading the reanalysis data is to pick out which variables I'd like to extract, and the list is extensive. I took that list and pasted them into the excel file that also has some of the main icemodel vars. The goal is to compare which variables from merra2 match those in icemodel.
vars_compare.xlsx

1. Matt, can you take a quick glance at the excel file and see if the green "matches" and yellow "maybe worth matching up" are accurate?
   Later, I realized that there are waaay more variables in icemodel that are loaded in, e.g., PHYSCONS. For many of those constants, I'm not sure where to look in order to get their definition. Do I need to account for them?
   The ultimate question is which of those variables in the reanalysis data match with icemodel variables // which ones should I tick?

2. I decided to skip ahead and just select 5 variables (green ones in excel sheet). I then wanted to download the data but found the only output to be netCDF (a large .txt file with lines of individual .nc files). Although there are directions on how to read the data, I am not sure how to approach this [UPD: resolved!]. The given directions – although for Mac – are meant to be executed in... Unix/shell (aka just the terminal?)? I tried running some of their lines and got an error at step 2 (type touch.netrc). Am i supposed to do it in Python? Can I find similar codes in matlab? Am i not looking hard enough in your code how you dealt with this? I saw some of your "ncread" commands but couldn't quite understand which variables/values you extracted, Matt.

[altanshapovalov]

UPD

I was able to download the nc files (thanks again, Johnny!); the directions were good, I just neglected a 'space'.
What I am currently slightly stuck with is exporting the data off .nc-files and putting into a neat timetable like "ice1" and "met".
I can use the ncread command and extract the numbers, but the output is a "1x1x24 Double"; I have 366 .nc-files, so each file is for 1 day in the year 2016, and file has 24 values for each hour. My problem is to extract those 24 numbers and put them in 24 rows of mr2 (my equivalent of Matt's timetables "ice1" and "met"), then call in the next .nc-file and fill the next 24 rows, etc.
I know that I need to use the power of wildcards and nested loops, but I can't get the details right just yet)

Below is my code (couldn't attach as an .m-file so it's in .txt format). I'm also attaching the folder with the data (zipfile; in case you want to give it a go). Fyi, inside "icemodel" (that contains folders like drive, functions, input, etc), I created my own folder called "drive_python", in which I put the folder "test1_5vars". This is why my first executable line of code is like that
test1_5vars.txt
test1_5vars.zip
This should be within my abilities to figure out, but I'm curious if you could push me in the right direction again so that I don't waste too much time on it.

And my question regarding which variables I should ultimately pick out still stands. I speculate that the 5 that I chose are insufficient.

[mgcooper]

Hi Max, Glad to see you are organizing this information. I have a similar spreadsheet attached here. It is not well organized. The most relevant sheets are ‘Greenland’ and ‘Forcing’. The Forcing sheet has the variables needed for forcing icemodel (not evaluating the output). There is a note that the incoming longwave is missing from the table. In the past I used LWGNT. I also made some notes on your spreadsheet. The ‘Greenland and Global’ sheet lists additional variables useful for evaluating the model. I think this is what you want. Note that Merra2 has “instantaneous” and “time averaged” values. You want time averaged. Download the Merra2 technical document and compare the filename definitions with the ones in my spreadsheet. That will clarify which variables are useful in general (which “collections” as they’re called). I attached a script to help download the data. It doesn’t actually download the data, it builds a shell script ‘wget_list.sh’ that does. It actually builds four scripts b/c there are four “collections”. The final step is to put those sh scripts wherever you want to save the data and then run them. I attached one of the sh scripts as an example, compare it to the text file you got and you should see how it works, but the scripts won’t work anymore b/c Merra2 filenames change over time when revisions are made. There is a better way you could do it though. Instead of writing the wget_list.sh scripts, you could just execute the lines that are written to that script in the loop using the matlab ‘system’ command or in python I think it’s called ‘subprocess’. So in summary you would just read in that text file and loop over each line, compose the wget command, and send it to ‘system’ or ‘subprocess’ (each line downloads one file, so your just writing a loop to download all the files). This is a good skill set to have so I suggest spending some time figuring it out. Good luck! From: max ***@***.***> Date: Saturday, November 5, 2022 at 9:53 PM To: mgcooper/icemodel ***@***.***> Cc: mgcooper ***@***.***>, Author ***@***.***> Subject: Re: [mgcooper/icemodel] Python conversion stuff (Discussion #2) Hey Matt, I'm a bit stuck (with downloading data) and want to get clarification (on variables) My current task is to download the reanalysis data (specifically, M2T1NXFLX from here<https://disc.gsfc.nasa.gov/>). Johnny gave me really good details on the directions (attached), but there are 2 points that I need help clearing up. The former is more directed for you, Matt. merra2-notes.txt<https://github.com/mgcooper/icemodel/files/9945172/merra2-notes.txt> Preamble: The penultimate step before downloading the reanalysis data is to pick out which variables I'd like to extract, and the list is extensive. I took that list and pasted them into the excel file that also has some of the main icemodel vars. The goal is to compare which variables from merra2 match those in icemodel. vars_compare.xlsx<https://github.com/mgcooper/icemodel/files/9945146/vars_compare.xlsx> 1. Matt, can you take a quick glance at the excel file and see if the green "matches" and yellow "maybe worth matching up" are accurate? Later, I realized that there are waaay more variables in icemodel that are loaded in, e.g., PHYSCONS. For many of those constants, I'm not sure where to look in order to get their definition. Do I need to account for them? The ultimate question is which of those variables in the reanalysis data match with icemodel variables // which ones should I tick? 2. I decided to skip ahead and just select 5 variables (green ones in excel sheet). I then wanted to download the data but found the only output to be netCDF (a large .txt file with lines of individual .nc files). Although there are directions<https://disc.gsfc.nasa.gov/data-access#mac_linux_wget> on how to read the data, I am not sure how to approach this. The given directions – although for Mac – are meant to be executed in... Unix/shell (aka just the terminal?)? I tried running some of their lines and got an error at step 2 (type touch.netrc). Am i supposed to do it in Python? Can I find similar codes in matlab? Am i not looking hard enough in your code how you dealt with this? I saw some of your "ncread" commands but couldn't quite understand which variables/values you extracted, Matt. — Reply to this email directly, view it on GitHub<#2 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACKRV6GE5GREWWOCFKCEF53WG42UJANCNFSM6AAAAAAQ2QPN64>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[altanshapovalov]

Hey Matt, quick question:
For your input weather station data, you have options for 15min and 1hr data. Based on some info I found, AWSs provide hourly data.
So, did you have to interpolate the hourly data to create a 15min subset or did you just download it as 15min data?

[mgcooper]

Interpolated. The 15 m is for numerical stability in the model … topic for another discussion. From: max ***@***.***> Date: Monday, December 5, 2022 at 11:27 AM To: mgcooper/icemodel ***@***.***> Cc: mgcooper ***@***.***>, Author ***@***.***> Subject: Re: [mgcooper/icemodel] Python conversion stuff (Discussion #2) Hey Matt, quick question: For your input weather station data, you have options for 15min and 1hr data. Based on some info<https://essd.copernicus.org/articles/13/3819/2021/> I found, AWSs provide hourly data. So, did you have to interpolate the hourly data to create a 15min subset or did you just download it as 15min data? — Reply to this email directly, view it on GitHub<#2 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACKRV6C5A4H76WGI5CQFG3DWLY6ZLANCNFSM6AAAAAAQ2QPN64>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[altanshapovalov]

Gotcha.
I now have a less quick question)

I wanted to play with the albedo and see how it would affect Runoff.
I took Merra2 albedo (M2T3NXGLC), interpolated it from 3hr to hourly (so that it could be the same size as ice1's 8784 rows), and wanted it to be the data that will replace met's/ice1's albedo (since ice1.albedo = met.albedo); fyi, I called that albedo m2a.albedo.

After a lot of code-perusing and head-scratching, I decided to replace met.albedo with m2a.albedo in the beginning of "icemodel.m", right after "METINIT(opts)" gets run. So, let me first you the default of IceModel running and showing Runoff changes. Here, the albedo is that of KAN M:

Now, here's me replacing met.albedo with m2a.albedo:

I also couldn't help but notice the other Runoffs in the figure, one of which is MERRA. After more digging, I found that MERRA is actually "MER", which loads "merra_behar_2016" (which is different from merra_kanm_2016).
Comparing the different albedos, I quickly figured out that m2a.albedo = merra_kanm_2016 albedo, and not the MER (aka merra_behar_2016) albedo.

My initial question/concern: since the MER and m2a albedos are so similar, then why do they produce such different runoffs? (based on the second figure)
Also, is it bad that I'm still calling things KAN M when we're actually using behar? How truly different are they? (let me know if you explained this somewhere already, i wouldn't want you to repeat yourself again)

Lastly, I think that my m2a is close to merra_kanm (which it obviously should be) but my data clearly has an extra spike or two. because of this, should I ditch my dataset and use what you have already?

[JohnnyRyan1]

I can have a go at answering these... I think it's likely that Behar is the average MERRA albedo for a 60km2 watershed upstream of KAN-M. Since most of these grid cells are higher elevation than KAN-M, they are slightly brighter, hence slightly less runoff than your single MERRA grid cell. There could be a number of reasons why m2a does not exactly equal merra_kanm. It could be older version of MERRA2, there could also be some smoothing applied to remove the spikes. I wouldn't worry too much about it. Overall this is good progress and I'm glad you're getting some runoff outputs.

[mgcooper]

I won’t be able to answer until Wednesday, but just a few ideas: * Johnny’s explanation sounds like the answer – you may be comparing a watershed-average to point-scale values * My forcing files are watershed-averages for the watershed files, and point-scale values for the met station files * Take some time to study the README.md and what happens in METINIT.m (also icemodel_opts where the met file name is built) to understand the forcing files because there are two types – “Met” files, and “Data” files * To understand that last part, pay attention to the file-naming convention in the README and the ‘swapvar’ concept and how it’s implemented in METINIT * Also try putting a debug stop in icemdoel_opts prior to the part where the filenames are built and/or in METINIT then step through one line at a time and see how the vars get swapped out, you can pass in your merra albedo as an extra variable to the function so you can have it in your workspace and you can plot it against the values in my files to see how they compare, but probably easier to do this outside of the model * Overall I think Johnny may be right but I need to read through your email more carefully This is a good exercise to see how the model works. I am a little confused by the instability in the KAN-M output. The runoff should be smooth like the others. But we’ll have to revisit that. Also, I think Johnny’s explanation is right as far as why the albedo timeseries are different, but the reason the runoff is so different must be related to something else since your m2a matches my merra kanm albedo. Again, you will have to study the way the swapvar concept works to understand what’s going on. Depending on your setup, you may be using KAN-M forcing data with MERRA albedo, or you may be running MERRA forcing. It is hard to know exactly what is going into the forcing file without seeing the full setup and/or stepping through the initialization. From: Johnny Ryan ***@***.***> Date: Monday, December 5, 2022 at 4:17 PM To: mgcooper/icemodel ***@***.***> Cc: mgcooper ***@***.***>, Author ***@***.***> Subject: Re: [mgcooper/icemodel] Python conversion stuff (Discussion #2) I can have a go at answering these... I think it's likely that Behar is the average MERRA albedo for a 60km2 watershed upstream of KAN-M. Since most of these grid cells are higher elevation than KAN-M, they are slightly brighter, hence slightly less runoff than your single MERRA grid cell. There could be a number of reasons why m2a does not exactly equal merra_kanm. It could be older version of MERRA2, there could also be some smoothing applied to remove the spikes. I wouldn't worry too much about it. Overall this is good progress and I'm glad you're getting some runoff outputs. — Reply to this email directly, view it on GitHub<#2 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACKRV6DMFZHN564XJTEKNTTWL2AYTANCNFSM6AAAAAAQ2QPN64>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[altanshapovalov]

Nested message incoming:
Okay, after tweaking some things in plotRunoff [mainly I changed the x axis that was Q.t (Q.t = time of discharge that was available only from july 6th to the 13th] to June 1st - Sept 1st to see the entire summer), I was able to yield this graph:

However, i would think that for the period of July 7-13th, the lines would look identical to the zoomed in version that icemodel provides by default:

Did I simply mislabel some lines or did I just plot something else entirely? (I suppose the latter is a question that should be more addressed at me since i'm the one who tweaked the code)
Okay, I think that one of those lines might be the ADCP that I thought I removed (and i removed it because it needed things like Q.err that I commented out).
In any case, can you spot out which line(s) got mislabeled or should I just go back to the code and pay more attention to the mess that I cooked up? Comparing Fig 2 to Fig 1 and looking at approximately July 13 (end of figure 2), The top line is correctly labeled RACMO and is blue. But then I think that in fig 1, red and black are switched; the red (MAR) should actually be KANM and vice versa. Lastly, I'm assuming that MERRA got labeled correctly.

[altanshapovalov]

Oh no, I think I figured out where I can check the lines AND isolate just KANM; somehow I figured this out while washing the dishes. I'll try to answer my own question tomorrow morning.

[altanshapovalov]

Hey @mgcooper , I have an odd problem. I forced the icemodel with merra2 albedo and i was able to run the model just fine, seeing the increase in runoff over the years. However, when I tried forcing it with modis albedo, the model ran but showed just a fat zero for the runoff change. I was able to determine that the error happens within icemodel.m, before the POSTPROC function. So, somewhere in the loop. Specifically, I noticed that the variables k_eff and T, by the end of the loop running, were all zero. I tried running the first iteration and comparing it simultanesouly with the successful merra2 albedo run, and the first iteration's answers were identical. So, it must be something wrong in the further iterations.
Do you have any idea why icemodel might be rejecting the modis albedo that I'm feeding it? Same size, no gaps/NaNs, just different albedo values. My speculation is that it doesn't like that it starts with "1" rather than some fraction like 0.849, but not sure if it cares for that at all. Perhaps you could tell by the graphs of albedo. Otherwise, I will also send
modis_simp_timetable.mat.zip
you my .mat file with the 'simplified' modis albedo (in case you want to try forcing it yourself; fyi, i just replaced met.albedo with this in the beginning of icemodel; this worked for merra2 albedo).
In the meantime, I want to put in some error message inside the loop that will yell at which iteration the k_eff turns to zero.
Sorry for the lack of formalities, this error is puzzling me deeply; hope you met the new year well)

[mgcooper]

Hi Max. If the modis albedo equals 1.0 that may be the problem. The albedo should never equal 1.0, it should always be between 0 and 1. If the albedo is strictly between 0 and 1 with no gaps/Nan, then I'd have to look in more detail to figure it out.

There used to be a series of checks in METINIT that checked for these types of things but I may have removed it because I have other functions that build the met files and do that checking. If you start building your own met files, we should discuss how to do that using the functions I wrote. Similar issues can arise if the wind speed is ever equal to zero or if the relative humidity is very very small.

As a side note, normally if T is all Nan, nearly every other variable defined on the thermal grid will be nan, but it may depend on where the solution is within the timestep. If it's just keff and T, that may be because keff gets updated first and the error occurs immediately after, prior to any other variables being updated.

Try setting any albedo values that equal or exceed 1.0 to 0.98 and see if that fixes it. Also, I cannot remember if we discussed this but I noticed last time that the modeled runoff wasn't very "smooth" and meant to remind you that icemodel needs 15 m forcing whereas skinmodel can run with hourly forcing. For testing purposes, I suggest running skin model. It only takes about 5 seconds to run one year. Have you run skin model yet?

[mgcooper]

Also - if you wanna jump on a Zoom sometime soon we can walk through running skinmodel b/c I am pretty sure the repo you cloned did not have a working version

[altanshapovalov]

Thanks for the tip – getting rid of the albedo values of 1 solved the problem! For the FINESST proposal that I'm writing right now, johnny remarked that it might be best if i stick to using icemodel for this since i'm trying to demonstrate icemodel's capabilities. So, since i don't need to worry about running skinmodel just this second, let's reschedule our zoom meeting for it sometime in the future when i'll need to get it going; I do appreciate the notion) However, I am planning to send you a rough draft of this proposal either thursday or friday this week. It's due on Feb 7 so hopefully by next week I could finalize it more or less.

[altanshapovalov]

On another note, @mgcooper any tips on extracting albedo data for a specific lat/lon from a 4D double? I tried looking at your scripts like scripts/preproc/met/mk_metfile_MAR_point.m and you have some kind of magical projfwd and projsipsn commands that turn your lat lons to x/y? Not usre if that will help but I remember you mentioning being cautious with reading nc files and how my loops are constructed but back then i just either 3d or 2d doubles.

Edit (for context): I have just one .nc file and it's for 2016. In there, the dimensions are 366 (time), 150 (x), 270 (y), and some sectors. The LAT and LON attributes are both 150x270 (x,y). AL2 (albedo) has the size of 150x270x1x366 (dimensions: x,y,SECTOR1_1,TIME). Should I... put AL2 into a variables (using ncread), then separate the x and y, and pick out the lat lon?

[mgcooper]

Hi Max, To read merra2 data files, you can use the function linked below. It depends on ‘ncparse’ which is linked below that. If you’re just reading albedo, you can ignore the comments in readMerra2.m regarding unit conversions. Getting the data for a particular location is slightly more complicated, especially using matlab. I recommend using a tool designed to work with netcdf files such as nco. Another option is cdo, but I am less familiar with it. I also suggest to download Panoply for quickly viewing netcdf data. There’s a bit of a learning curve with nco but well worth it. I am sure rasterio or xarray would also work if you want to use python. I always used to use the ‘raster’ package in R and loved it, but nowadays not so much. If you’re committed to using matlab, I can provide some tools but they’re not very user friendly. Basically, I read in the merra lat/lon variables and use the ‘inpolygon’ function to figure out which lat/lon values are inside a closed area, then I read in the entire data grid, and mask out the data using the result from inpolygon. Don’t be dumb like me though! Instead use nco, or xarray. https://github.com/mgcooper/matfunclib/blob/main/libhydro/readMerra2.m https://github.com/mgcooper/matfunclib/tree/main/libdata/readfiles/ncparse From: max ***@***.***> Date: Wednesday, January 25, 2023 at 5:43 PM To: mgcooper/icemodel ***@***.***> Cc: Matt Cooper ***@***.***>, Mention ***@***.***> Subject: Re: [mgcooper/icemodel] Python conversion stuff (Discussion #2) On another note, @mgcooper<https://github.com/mgcooper> any tips on extracting albedo data for a specific lat/lon from a 4D double? I tried looking at your scripts like scripts/preproc/met/mk_metfile_MAR_point.m and you have some kind of magical projfwd and projsipsn commands that turn your lat lons to x/y? Not usre if that will help but I remember you mentioning being cautious with reading nc files and how my loops are constructed but back then i just either 3d or 2d doubles. — Reply to this email directly, view it on GitHub<#2 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACKRV6ABTPKIBRDM34UQPRDWUHJDBANCNFSM6AAAAAAQ2QPN64>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[mgcooper]

PS here’s a good set of tips on nco. The linked example may do what you want. http://research.jisao.washington.edu/data_sets/nco/#example10 From: Matt Cooper ***@***.***> Date: Monday, January 30, 2023 at 7:02 PM To: mgcooper/icemodel ***@***.***>, mgcooper/icemodel ***@***.***> Cc: Mention ***@***.***> Subject: Re: [mgcooper/icemodel] Python conversion stuff (Discussion #2) Hi Max, To read merra2 data files, you can use the function linked below. It depends on ‘ncparse’ which is linked below that. If you’re just reading albedo, you can ignore the comments in readMerra2.m regarding unit conversions. Getting the data for a particular location is slightly more complicated, especially using matlab. I recommend using a tool designed to work with netcdf files such as nco. Another option is cdo, but I am less familiar with it. I also suggest to download Panoply for quickly viewing netcdf data. There’s a bit of a learning curve with nco but well worth it. I am sure rasterio or xarray would also work if you want to use python. I always used to use the ‘raster’ package in R and loved it, but nowadays not so much. If you’re committed to using matlab, I can provide some tools but they’re not very user friendly. Basically, I read in the merra lat/lon variables and use the ‘inpolygon’ function to figure out which lat/lon values are inside a closed area, then I read in the entire data grid, and mask out the data using the result from inpolygon. Don’t be dumb like me though! Instead use nco, or xarray. https://github.com/mgcooper/matfunclib/blob/main/libhydro/readMerra2.m https://github.com/mgcooper/matfunclib/tree/main/libdata/readfiles/ncparse From: max ***@***.***> Date: Wednesday, January 25, 2023 at 5:43 PM To: mgcooper/icemodel ***@***.***> Cc: Matt Cooper ***@***.***>, Mention ***@***.***> Subject: Re: [mgcooper/icemodel] Python conversion stuff (Discussion #2) On another note, @mgcooper<https://github.com/mgcooper> any tips on extracting albedo data for a specific lat/lon from a 4D double? I tried looking at your scripts like scripts/preproc/met/mk_metfile_MAR_point.m and you have some kind of magical projfwd and projsipsn commands that turn your lat lons to x/y? Not usre if that will help but I remember you mentioning being cautious with reading nc files and how my loops are constructed but back then i just either 3d or 2d doubles. — Reply to this email directly, view it on GitHub<#2 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACKRV6ABTPKIBRDM34UQPRDWUHJDBANCNFSM6AAAAAAQ2QPN64>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[mgcooper]

also see:
https://github.com/nco/nco

[altanshapovalov]

Hey Matt,

What is the difference between Ice1 and met Tsfc? I tried to look at the POSTPROC.m function to understand how the two differ, but i'm not sure which would be better for me to use. For context, I want to juxtapose summer temperature and runoffs. Here's a screenshot of POSTPROC where the Tsfc of ice1 and met get updated (the highlighted lines 39 and [order] 66) for convenience.

I also tried plotting them against one another and the difference seems to be minuscule. I'm guessing that I don't need Tair and that Tf is the freezing temperature

upd: JR said that Tair is actually what I need, which solves my conundrum, but I'm still curious how the calculation between Tsfc for ice1 and met differs

[mgcooper]

Hi Max, The ice1.Tsfc variable is the surface temperature simulated by icemodel. The met.Tsfc variable is the Tsfc variable from the met data source. For example, if a MAR metfile is used, then met.Tsfc is the surface temperature simulated by MAR. It isn’t used in the model, it is included for comparison with the icemodel simulated surface temperature in ice1. Matt From: max ***@***.***> Date: Tuesday, April 25, 2023 at 11:28 AM To: mgcooper/icemodel ***@***.***> Cc: Matt Cooper ***@***.***>, Mention ***@***.***> Subject: Re: [mgcooper/icemodel] Python conversion stuff (Discussion #2) Hey Matt, What is the difference between Ice1 and met Tsfc? I tried to look at the POSTPROC.m function to understand how the two differ, but i'm not sure which would be better for me to use. For context, I want to juxtapose summer temperature and runoffs. Here's a screenshot of POSTPROC where the Tsfc of ice1 and met get updated (the highlighted lines 39 and [order] 66) for convenience. [Screen Shot 2023-04-25 at 11 24 22]<https://user-images.githubusercontent.com/89253609/234367943-be4e5d60-5be7-4efa-ac98-c5fccbb28812.png> I also tried plotting them against one another and the difference seems to be minuscule. I'm guessing that I don't need Tair and that Tf is the freezing temperature [Screen Shot 2023-04-25 at 11 03 52]<https://user-images.githubusercontent.com/89253609/234368356-a6efe69c-de41-41de-80aa-d27bef0e3175.png> — Reply to this email directly, view it on GitHub<#2 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACKRV6GEAEXRKFDKOZMDON3XDAJVHANCNFSM6AAAAAAQ2QPN64>. You are receiving this because you were mentioned.Message ID: ***@***.***>