Maximising asif skills using near real time data

[SimonF92]

Im writing it based on their dT30 data from here, should be very easy to do and resistant to seasonal temperature changes:

https://data.meereisportal.de/download/buoys/2019T56_lastprofile.png

As you said, its why they included these parameters in the first place

//edit//
scratch that, ill do it on the dT120 instead, even more obvious

[uniquorn]

If it's not too much trouble, please include T1 in the final csv. It will make it easier to run temp vs thickness analysis

[SimonF92]

New code coming to you. I feel like this one might underestimate, but it is based on dT120 which has a very distinct and clear ice boundary.

The question is, IS it underestimating or it is that thickness?

Their data here,

https://data.meereisportal.de/gallery/index_new.php?lang=en_US&active-tab1=method&active-tab2=buoy&singlemap&buoyname=2019T56

//Edit//

 having read about it I have to say that I think its the real boundary

[SimonF92]

So we are getting closer, the ice:snow boundary is now very clear.

Can you see the bump labelled green? Thats something were discussing. That could either be a sharp melting event, or a bug.

As the air temperatures do not rise to facilitate this melt, is probably a bug.

The best way to tell is to slice those dates off the dataframe. If those bumps move with the slice, to a different date, then it is definitely a bug.

If the bumps disappear, then it is not a bug.


SUMMARY: ice thickness is blue-red lines

[SimonF92]

It wasnt a bug.

Any thoughts? Influx of warm water facilitating bottom melt?

[uniquorn]

Looking only at T64 where the chart shows 3 clear drops in estimated thickness. The temperature ani doesn't show any bottom melt events. There are, however, some higher air temperatures that flatten the curve that may occur at roughly the same time as the 3drops.
Not exactly a bug but perhaps an error in the method.
Some thermistors are recording persistently higher temperatures. Perhaps there was a problem during calibration
Was using Spyder but jupyter charts are good

Still need to see a time series of the stdev(10rows) charts to see it they are messing about during temperature rises

[SimonF92]

Its interesting you picked up on t64, i was also looking at it closely.

I used their proprietary plots to check ours. It seems there are major fluxes of heat (disruption) into the column at the same time we detect "melting"- its possible this is some kind of ocean current?

I dont think this can be overcome at our end at it is data- based.

https://data.meereisportal.de/gallery/index_new.php?lang=en_US&active-tab1=method&active-tab2=buoy&singlemap&buoyname=2019T64

Sorry for the low resolution

[uniquorn]

Those fluxes are in the atmosphere. I overlayed T64 thermistor1 onto daily temperature. Maybe our algo should detect air/snow first. There seems to be a defined line there, then snow/ice, then ice/water. Possibly only using mean heat temp change.

[SimonF92]

 Agreed- problems with the code and the temp anomalies during these "warm" events. Will figure something out.

Lovely overlay by the way, when this is tied up id love to get a look at some of your code/methods.

Its also on my agenda to do stdev timeseries, though that adds a new dimension to the data im a bit naive about (think of a csv file where each cell has depth) will try to get it working.

[uniquorn]

I was thinking of turning the 126 t64heat stdev images into an animation, just to get a feel for the movement.
Overlay was done in gimp. Drag both images in, measure the width of the charts in pixels, scale then set white to transparent.

[SimonF92]

I will have a bash, I think your skillset is better for animations than mine for sure, id love to see how you use gganimate in R.

I am working on this algo at the moment. Its going to have to be very well written, the changes are subtle but the boundary is there.

[SimonF92]

Ok so its definitely getting a bit more complicated but hopefully what im doing is logical.

The idea above would not work as its just too difficult to resolve. Instead, I am mapping the thermistors

RED: at the -5degree isotherm. (not perfectly mapped yet, currently is between -4 and -6 degrees)

GREEN : our perceived ice:ocean interface

BLUE :  a savitsky-golay fit on the red -5degree isotherm. As warm transient air temps penetrate the ice, they raise our RED isotherm (which is a data quirk and not real, the ice does not have bottom melt during this), which then disrupts our GREEN line. I will correct the GREEN against the delta (change) between RED and BLUE.

Will see how it looks. Im not going to correct until im sure iv got the best savgol method done

[uniquorn]

got these on the way to making the animation. Using some of your py code with this below               
        #try to save plot
        plt.plot(df1.loc['mean_t'])
        plt.ylabel('mean')
        plt.savefig("meantst5.png", dpi=400)
        #plt.plot(df1.loc['stdev'])
        #plt.ylabel('std')
        #plt.savefig("stdtst5.png", dpi=400)

At least 2 stdev stand out in what I think is the snow layer. I'm not sure what mean_t is yet but it looks the most useful.
click for full res

Check here for original macidR gganimate code and these nice 3d buoy presentations. https://rpubs.com/macid

[SimonF92]

got these on the way to making the animation. Using some of your py code with this below               
        #try to save plot
        plt.plot(df1.loc['mean_t'])
        plt.ylabel('mean')
        plt.savefig("meantst5.png", dpi=400)
        #plt.plot(df1.loc['stdev'])
        #plt.ylabel('std')
        #plt.savefig("stdtst5.png", dpi=400)

At least 2 stdev stand out in what I think is the snow layer. I'm not sure what mean_t is yet but it looks the most useful.
click for full res

Check here for original macidR gganimate code and these nice 3d buoy presentations. https://rpubs.com/macid

Nice plots and thanks for the link. Didnt know you could use the loc function in plt, thats cool. The algos snowy and thicky (partially) are based on the stdev plot.

Mean_t is the mean temperature at that measurement point (each per 10-single measurements) for each thermistor. Id be careful of using mean_t as we expect that to become very similar to the ice temperature in the coming weeks.

Still working on this smoothing, its driving me mad

[uniquorn]

I should have included
import matplotlib.pyplot as plt

Mean_t is the mean temperature at that measurement point (each per 10-single measurements) for each thermistor.

That's what I assumed, but why are the values from 0-5?
mean difference?

[SimonF92]

The values are between 0 and 5 as HEAT is not real temp, its some fancy thing they do to measure the heat capacity of that thermistor.

uniquorn what do you think of this?

We have an expression in our lab "cant see the forest for the trees". I have become so bogged down in trying to iron out the quirks in the data recently I lost sight of the big picture. Here is the end result? How would you feel about using these in your thickness maps?

[SimonF92]

I have actually been trying to get it working in python but it looks terrible compared to what you had

[uniquorn]

I get it. It's the mean temp rise over 120s. In which case it should still be valid as temps rise.

That's not so bad but could be misleading without dates. edit: many lat/lon's are bad on the temp_proc files, it's better to use the _ts file for location.
If those thicknesses compare closely enough with the basic method I'd be happy to use them though as mentioned previously I'm not too keen on smoothing.

added heat file ani for reference

[SimonF92]

I get it. It's the mean temp rise over 120s. In which case it should still be valid as temps rise.

That's not so bad but could be misleading without dates. edit: many lat/lon's are bad on the temp_proc files, it's better to use the _ts file for location.
If those thicknesses compare closely enough with the basic method I'd be happy to use them though as mentioned previously I'm not too keen on smoothing.

added heat file ani for reference

I completely agree, I am not a fan of smoothing under any circumstances.

In that case ill make a "basic method" script for use now and keep tweaking away at this over the coming weeks in the background.




///////

this is the basic method?

"""
Cell formula will be   =COUNTIF(D2:II2,"<-2.1")
We omit column IJ. I think it is some sort of end of line check.

Now subtract Tas=40
Cell formula will be   =COUNTIF(D2:II2,"<-2.1")-40

multiply by 2 because the thermistors are 2cm apart.
Cell formula will be   =(COUNTIF(D2:II2,"<-2.1")-40)*2

"""

[uniquorn]

Yes, except 40 only applies to T62. Need to identify the Therm number for the others.

[uniquorn]

Estimated thermistor numbers based on last entry in heat120 files (chart attached, cffr)
buoy   Th(snow)   Th(ice)
T56   24   35
T58   18   34
T62   40   51
T63   18   32
T64   24   47
T65   12   35
T66   42   49
T68   34   46
T70   32   42
T72   28   36

T69   25   35

[uniquorn]

Rough charts based on the data above. Some different curves.
Note that the heat120 method appears to shift the air-snow and snow-ice thermistor no.s upwards, possibly due to heat rising. I haven't corrected for this. Even so, that doesn't explain 0 or negative starting position unless data was sent before (or during) buoy installation. hey ho, it's an estimate.

[SimonF92]

Rough charts based on the data above. Some different curves.
Note that the heat120 method appears to shift the air-snow and snow-ice thermistor no.s upwards, possibly due to heat rising. I haven't corrected for this. Even so, that doesn't explain 0 or negative starting position unless data was sent before (or during) buoy installation. hey ho, it's an estimate.

Im actually working on a "basic method" script at the moment.

Within that script is the algo for detecting Th(snow) and Th(ice) thermistors .They look similar to your estimates, though not identical, ill send you it within 30mins.

New points in the code;

-basic method implemented
-lat and lon now taken from TS buoys, corrected for sampling period (to 1 per day)

Im also working on an additional program that lets you manually enter  Th(snow) and Th(ice) thermistors, when you wish to plot. In some ways I think the human eye might actually be better at this than an algorithm


//////

FYI all the code is pretty horrible at the moment, im reading a book called

 "Clean Code, A Handbook of Agile Software Craftmanship"

Once were happy, ill be tightening everything up very strictly

[SimonF92]

Iv added some new outputs to the console for you also;

[uniquorn]

Thinking about those rough charts, it would be easier to compare them if they all started on the same date. I think the shallow gradient was during the warm spell in October, some of the Tbuoy's data don't start till late oct/nov. Of course, if the ice is thin, it thickens quicker.

[SimonF92]

Thinking about those rough charts, it would be easier to compare them if they all started on the same date. I think the shallow gradient was during the warm spell in October, some of the Tbuoy's data don't start till late oct/nov. Of course, if the ice is thin, it thickens quicker.

Good point? I can slice it at a date if you think that would work?

December 1st?


//////

scratch that, you can change it to any date you wish within the code

[uniquorn]

ok. I added thickness*100 to the charts and got this
and I think I'm mixing Th numbers with thickness

[SimonF92]

ok. I added thickness*100 to the charts and got this
and I think I'm mixing Th numbers with thickness

I dont think you are, the easiest way to change to thickness from m to cm is to change:

 #ice thickness is therefore the number of these thermistors
        thickness=(len(therms)*2)/100

at line 254 on github to:
 
 #ice thickness is therefore the number of these thermistors
        thickness=(len(therms)*2)

and nothing else

[SimonF92]

My output on the action above

[uniquorn]

Comparison of ice thermistors from the two methods.
Buoy   Th(ice)   Th(ice)SF
T56   35   38
T58   34   31
T62   51   47
T63   32   29
T64   47   45
T65   35   32
T66   49   46
T68   46   44
T70   42   40
T72   36   33

I'll run an ani with the SF ice therms.

[SimonF92]

Here is a very quick way of checking my thermistor method, just stick it in its own cell and it will do its thing


////////////////////////////////////////////////////////////////////


import pandas as pd

class buoys:
   
    def __init__(self,url):
       
        self.url=url
        self.df=pd.read_csv(self.url)
       
    def thice(self):
       
       
       
        df_HEAT=self.df
        ##########
        #possible integer to change below
        ##########
        df_HEAT=df_HEAT.tail(10)
        df_HEAT=df_HEAT.drop(['time', 'latitude (deg)','longitude (deg)'], axis=1)       
        df_HEAT.loc['stdev'] = df_HEAT.std()
        df_HEAT.loc['mean_t'] = df_HEAT.mean()       
        df_HEAT=df_HEAT.tail(2)
        df_HEAT=df_HEAT.T
        df_HEAT.reset_index(level=0, inplace=True)
        df_HEAT["Temp_Change"]=df_HEAT.mean_t.diff()
        ##########
        #possible integer to change below
        ##########
        df_ocean=df_HEAT.tail(50)
        core_temperature=df_ocean.mean_t.mean()
        df_HEAT=df_HEAT.head(-5)
        df_HEAT["Deltacore"]=df_HEAT.mean_t-core_temperature
        ##########
        #possible integer to change below
        ##########
        thermistors_not_in_ice=df_HEAT.where(df_HEAT.Deltacore>1)
        thermistors_not_in_ice=thermistors_not_in_ice.dropna()
        Th_ice= thermistors_not_in_ice.tail(1)
        Th_ice=list(Th_ice.index.values)
        Th_ice=int(Th_ice[0])
       
        self.Th_ice=Th_ice
        return(self.Th_ice)
       
       
T56=buoys("https://data.meereisportal.de/download/buoys/2019T56_300234065176750_HEAT120_proc.csv")
T58=buoys("https://data.meereisportal.de/download/buoys/2019T58_300234065171790_HEAT120_proc.csv")
T62=buoys("https://data.meereisportal.de/download/buoys/2019T62_300234068706290_HEAT120_proc.csv")
T63=buoys("https://data.meereisportal.de/download/buoys/2019T63_300234068709320_HEAT120_proc.csv")
T64=buoys("https://data.meereisportal.de/download/buoys/2019T64_300234068701300_HEAT120_proc.csv")
T65=buoys("https://data.meereisportal.de/download/buoys/2019T65_300234068705730_HEAT120_proc.csv")
T66=buoys("https://data.meereisportal.de/download/buoys/2019T66_300234068706330_HEAT120_proc.csv")
T68=buoys("https://data.meereisportal.de/download/buoys/2019T68_300234068708330_HEAT120_proc.csv")
T70=buoys("https://data.meereisportal.de/download/buoys/2019T70_300234068705280_HEAT120_proc.csv")
T72=buoys("https://data.meereisportal.de/download/buoys/2019T72_300234068700290_HEAT120_proc.csv")

print(T56.thice())
print(T58.thice())
print(T62.thice())
print(T63.thice())
print(T64.thice())
print(T65.thice())
print(T66.thice())
print(T68.thice())
print(T70.thice())
print(T72.thice())
   

[SimonF92]

Ps, technically my thermistors should be +1 on its value because it is theoretically the first thermistor ABOVE the ice

Not sure about the 1st one



//////

to fix that, add the line

Th_ice=Th_ice+1

underneath

"Th_ice=int(Th_ice[0])"

[uniquorn]

Will look at that in a bit. First draft of the ani. My initial impression is that it's a shame to lose the path detail by using only 1frame/day. How about we use all the ts data but replicate the thicknesses.
Note that none of the buoys are over 2m thick using this method.  ctr
I have a feeling they are all going to experience sudden melt at this drift rate.

adjusted a few things ctffr
and added est snow thickness, though I can't do 2 colours at a time.

[uniquorn]

The 2 remaining cryosphereinnovation simb buoys are in broad agreement with our numbers. Snow depth perhaps more stable.
piomas indicating our results should be a touch thicker

[SimonF92]

Looks great, fantastic job on the animation, I agree about the frames, I will retain all ts data "fill in the blanks" for thickness

[uniquorn]

Another thought, maybe thickness change would be better for path colour retaining thickness as label. I can do that in R though.

[SimonF92]

Another thought, maybe thickness change would be better for path colour retaining thickness as label. I can do that in R though.

you can also change the size of the points based on a quantitative variable such as thickness


//////


youll struggle to do thickness change in R now ive stacked the rows n times (where n is the number of TS measurements per day)

ill add it into the code so it operates before this line

[SimonF92]

I wrote a script for plotting individual buoys in a "nice" way. Based on the figure you shared above

[blumenkraft]

So we have an almost linear growth of ice thickness and changing temperatures don't dramatically change the growth? That's stunning. 

[SimonF92]

So it would seem. Id say im probably less informed than yourself and uniquorn about this, but its maintained across the buoys

Sorry about the legend clipping here, python is being a real pain



FYI that is an animated gif

[uniquorn]

Easier to see slow down in thickening, probably due to temp, with a 'taller' scale
btw SF, Nice emulation ctr

[SimonF92]

We deliberately cut the start day to Dec 1st, but when you take out that line the changes in growth become more obvious.

We put that line in because not all the buoys start transmitting on the same day so the data-structure gets messy



///////

please ignore the dates, they are now bugged- code is like a house of cards

[uniquorn]

20 frames of estimated change in thickness. Still tinkering. ctffr

[SimonF92]

20 frames of estimated change in thickness. Still tinkering. ctffr

Looking good.

Would you consider not using viridis as the colourmap for this plot? There are some diverging colourmaps that might be better for easily spotting whether its loss or gain. Im not sure if they translate directly into R but there is a list here;

https://matplotlib.org/3.1.1/gallery/color/colormap_reference.html


Also, for the temps (which is a great feature) I can give you a once-daily temp dataframe so they dont have flats at each day- something I guess is because there are repeating measurements due to the TS fix.

[uniquorn]

R not happy with discrete colourmaps with the current coding. Please suggest your favourite continuous.
Was lazy with the chart but wanted to fill up the blank space.

[SimonF92]

"icechange <- colourmap(c("red", "blue", "green"), breaks=c(float(max_loss),0,
                     (float(max_gain))"

May allow you to make a continuous colourmap of your own choosing. Im pretty bad with R though.

https://www.rdocumentation.org/packages/spatstat/versions/1.63-3/topics/colourmap

[uniquorn]

Looking at thickness change colours and incorporating some bathy (scale appears to be wrong). ctffr

[SimonF92]

Looking at thickness change colours and incorporating some bathy (scale appears to be wrong). ctffr

starting to look really good

are the buoys moving in frames of 10 measurements-per-movement at the moment?

[uniquorn]

After asking for more rows I went back to 1frame/day. Text repel labels jump about too much at that resolution. If you are happy with the numbers I'll post sthing similar to this on the mosaic thread as a reference. Subsequent posts can then be higher res, shorter time frame, more distance between the labels.
For me, the drift path is important so thanks for the extra coding.

[SimonF92]

No problem, I think it looks great- well done.

Id consider having a "thickness" figure, exactly the same as your "change in thickness" except where the colourmap change is just thickness (you've done it before in viridis). I think the two will complement each other and will provide different benefit to different people in the community