the Air Vent

Because the world needs another opinion

Sea Ice – Overlapping Satellite Record Analysis

Posted by Jeff Id on March 11, 2012

Despite criticisms from the super warmist advocacy guild (aka SWAG), this post is a continuation of our investigation into  sea ice from the satellite data.  I have never thoroughly vetted the knitting of satellite data from the 5 satellites which form the curves “n07″,”f08″, “f11″, “f13″ and  “f17″ are the designators.  Unfortunately, the task is quite daunting in that necessary data to do it thoroughly is not readily available on line.  Still, we can start with what we have, so to that end,  I have re-written the sea ice code again, this time the anomaly is taken over the whole range and leap years are accounted for.   It still isn’t terribly clean but it is improving.

ssq = function (x) {sum(x ^ 2)}

plotTrend = function (x, st = 1978, en = c(2012, 12), y.pos = NA, x.pos = NA, main.t = "Untitled") {

### Get trend
trend = lm(window(x, st, en) ~ I(time(window(x, st, en))))

### Initialize variables
N = length(window(x, st, en))
I = seq(1:N) / frequency(x)
rmI = ssq(I - mean(I))

### Calculate sum of squared errors
SSE = ssq(trend[[2]]) / (N - 2)

### Calculate OLS standard errors
seOLS = sqrt(SSE / rmI) * 10

### Calculate two-tailed 95% CI
ci95 = seOLS * 1.964

### Get lag-1 ACF
resids = residuals(trend)
r1 = acf(resids, lag.max=1, plot=FALSE)$acf[[2]]

### Calculate CIs with Quenouille (Santer) adjustment for autocorrelation
Q = (1 - r1) / (1 + r1)
ciQ = ci95 / sqrt(Q)

### Plot data
plot(window(x, st, en), main=main.t, ylab="Anomaly Km^2")
grid()

### Add trendline and x-axis
abline(h = 0, col = 2)
abline(trend, col = 4, lwd = 2)

### Annotate text
text(x = ifelse(is.na(x.pos), st, x.pos), y = ifelse(is.na(y.pos), max(x, na.rm = TRUE) - 0.2 * max(x, na.rm = TRUE), y.pos), paste("Slope: ", signif(trend$coef[2] * 10, 3), "+/-", signif(ciQ, 3), "Km^2/Decade\n(corrected for AR(1) serial correlation)\nLag-1 value:", signif(r1, 3)), cex = 0.8, col = 4, pos = 4)

###  Return
r = list(trend$coef[2] * 10, ciQ)
}

cal.Isig=function(dat=dat, st=0,en=2020)
{
###  Get trend
fm=lm(window(dat, st, en)~I(time(window(dat, st, en))))

###  Initialize variables
N=length(window(dat, st, en))
I=seq(1:N)/frequency(dat)
rmI=ssq(I-mean(I))

###  Calculate sum of squared errors
SSE=ssq(fm[[2]])/(N-2)

###  Calculate OLS standard errors
seOLS=sqrt(SSE/rmI)

###  Calculate two-tailed 95% CI
ci95=seOLS*1.964

###  Get lag-1 ACF
resids=residuals(lm(window(dat, st, en)~seq(1:N)))
r1=acf(resids, lag.max=1, plot=FALSE)$acf[[2]]

###  Calculate CIs with Quenouille (Santer) adjustment for autocorrelation
Q=(1-r1)/(1+r1)
ciQ=ci95/sqrt(Q)

return(c(fm$coef[2],ciQ))
}

getnorthicedata=function(file=file)
{
a=file(file,"rb")
header= readBin(a,n=300,what=integer(),size=1,endian="little",signed=FALSE)
data=readBin(a,n=304*448,what=integer(),size=1,endian="little",signed=FALSE)
close(a)
data
}

##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
#calculate northern sea ice area from gridded satelite data

Nfilenames=list.files(path="C:/agw/sea ice data/north/", pattern = NULL, all.files = TRUE, full.names = FALSE, recursive = TRUE)
startday=as.integer(unclass(as.POSIXct(strptime("1977-12-31", '%Y-%m-%d' )))/86400)
year= substr(Nfilenames,9,12)
month=substr(Nfilenames,13,14)
day=substr(Nfilenames,15,16)
dat=paste(year,month,day,sep="-")
oday=as.integer(unclass(as.POSIXct(strptime(dat, '%Y-%m-%d' )))/86400)-startday
yday=as.POSIXlt(dat)$yday  #day of the year. Used for anomaly calculations.
sat=substr(Nfilenames,18,20)
satlist=levels(factor(sat))

farea=rep(NA,length(Nfilenames))

#set up hole mask here
starthmask=!(getnorthicedata(paste("C:/agw/sea ice data/north/",Nfilenames[10],sep=""))==251)
endhmask=!(getnorthicedata(paste("C:/agw/sea ice data/north/",Nfilenames[9000],sep=""))==251)
diffmask=!(starthmask==endhmask)
holemask=starthmask #insures no step

lat=readBin("C:/agw/sea ice data/sea ice coordinates/psn25lats_v3.dat", integer(), endian = "little",n=304*448)/100000
lon=readBin("C:/agw/sea ice data/sea ice coordinates/psn25lons_v3.dat", integer(), endian = "little",n=304*448)/100000
latmask=lat>66.5  #arctic circle

#################################
#load area values into 2d matrix 366days each year x 35 days
for(i in 1:(length(Nfilenames)-1))
{
#open filename
fn=paste("C:/agw/sea ice data/north/",Nfilenames[i],sep="")
data=getnorthicedata(fn)
print(dat[i])
#Limit data/land mask to all info 25 - 10%
datamask=data25

#store sum in trend value
farea[i]=sum(data[(datamask&holemask)])/250*625
}

st=oday[1]
en=oday[length(oday)]

Ifarea = approx(oday,farea,st:en)#interpolated series for anom
plot(Ifarea,type="l")
yd=as.POSIXlt((st:en)*86400, origin="1977-12-31")$yday+1  #day of the year. Used for anomaly calculations.
yr=as.POSIXlt((st:en)*86400, origin="1977-12-31")$year-77  #day of the year. Used for anomaly calculations.

yrs=max(yr)
anomarry=array(NA,dim=c(yrs,366))
for(i in 1:length(yr))
{
anomarry[yr[i],yd[i]]=Ifarea$y[i]
}
an=colMeans(anomarry,na.rm=TRUE)
plot(an)

## North sea ice area, no infilling all satellites
NIA=farea
tm=1978+oday/365
## North sea ice area anomaly, no infilling all satellites
NIAA=NIA-an[yday+1]
plot(tm,NIAA,type="l")
## Assemble all satellites into timeseries
allNSAT=ts(array(NA,dim=c(max(oday),length(satlist))),start=1978,deltat=1/365.25)

The code now allows anomaly calculation using infilled data and plotting of data which has no infilling whatsoever in the main series. This is important because in the early 80′s and before, sea ice data was taken every other day. The trend in sea ice means that missing years have a different total offset than years with data. If not corrected, the result is to add a sawtooth pattern to the anomaly calculation making things less significant – and we don’t want SWAG on our butts.

What I have done is to take the data with its true time and arrange it into a linear series. That series is infilled with linear interpolation using the “approx” function in R. This infilled data is assembled into a matrix and a daily average is calculated for use in anomaly offset only. The original non-infilled data is then corrected for seasonal anomaly and is shown in the following plot.

Figure 1

The purpose of this was twofold: First, I want to get to a numerically clean solution for anomaly calculation. Second, I would like to investigate the performance of individual satellites.  In a single series the satellite trend is shown again with appropriate color coding.

Figure 2

By chance the first satellite I plotted was F13 by itself (dark blue).  It had a very unusual look so I spent the next two hours reading about diurnal corrections, more on that in a future post.   I then took a look at how cleanly the two satellites matched during the conversion point.

Figure 3

The between the two independent sensors is quite good but visually F17 is slightly above F13 at the beginning 2007 yet closer at the end 2008.   To investigate, we can plot the difference to see if there is any trend between the satellites.

Figure 4

Now that is interesting.   The new F17 satellite has a biased trend of about 20km^2/year that nearly reaches two sigma significance!  Sometimes we need to override the numbers with our minds though.  That is an UGLY breakpoint in the middle of the series, not impossible, not likely natural though and it is not likely a continuing error.  Now that change doesn’t represent much ice considering the peak 14 million km^2 of Arctic sea ice BUT that is 200,000km^2/decade which happens to represent about 40% of the total measured sea ice loss trend in the Arctic.  There is an unexplained net offset of 35000Km^2 in the data though which, if  corrected, would increase the net downtrend.

Figure 5 - Total Northern Hemisphere Sea Ice Loss. Please excuse the C/Decade in the graph - it should be km^2/Decade.

Now that is interesting.  It sure gives you an idea of just how small of a difference between these instruments is important.   These guys can’t make any mistakes or the whole trend can be messed up.


19 Responses to “Sea Ice – Overlapping Satellite Record Analysis”

  1. stan said

    Jeff, experts never make mistakes.

  2. Jeff condon said

    Stan,

    I think you must be mixing up Tamino with experts?!
    ======

    For others, there really isn’t any error here though that I can find yet. It is rather amazing that you can see two trends as closely matched as figure 4 yet the difference is 40% of the sea ice “doom” trend.

  3. Layman Lurker said

    Jeff, do you know why there is overlaping data posted between f13 and f17 but not the other satellites?

  4. Jeff condon said

    I believe other sats died before launch of the next. There is one other overlap point 11,13 and there are two poles so that should improve our fidelity. The problem is that variation in time of day may cancel at the poles leaving us with the impression that all is well, when it might not be.

  5. curious said

    Hi Jeff – I’m confused! In the line

    “There is an unexplained net offset of 35000Km^2 in the data though which, if corrected, would increase the net downtrend.”

    I think you are refering to the step at approx 2007.4? So if this were “corrected” and the RHS of F17-F13 was pushed up 35000km^2, wouldn’t the net downtrend decrease ie be less of a negative slope? Or are you saying F13 should be pulled down? Or have I missed the point!?

  6. Layman Lurker said

    I believe other sats died before launch of the next.

    Jeff, if you look at Wentz (2010) table 1 and 2, figure 1 and 2, Wentz has graphed the raw data differences between satellites during the overlapping periods.

  7. Another Ian said

    Jeff,

    I just had to go into compatibility mode to view comments here – like WUWT trouble

  8. Kevin O'Neill said

    That is an UGLY breakpoint in the middle of the series, not impossible, not likely natural though and it is not likely a continuing error.

    The NASA Intersensor Calibration between F13 SSMI 1 and F17 SSMIS plots show discontinuities between days 19-33, 64-73, 189-201 (approximations from viewing their graphs). I assume this is the missing data (they say 35 days total). The ‘breakpoint’ falls into one of the periods of missing data.

  9. Matthew W said

    Dose this factor in in or not related?

    http://www.real-science.com/arctic-fraud-worse

  10. Kevin O'Neill said

    Matthew – you shouldn’t believe everything you read. Remember the S.S.Manhattan – the icebreaking oil tanker? She sailed the Northwest Passage in 1969. She had to turn south through Prince of Wales Strait after encountering ice 35 meters thick. Since 2007 that same area is usually open water.

  11. DocMartyn said

    Jeff, R.E. Fig 3. Do the different satellites have a different variation in the signals?

  12. Layman Lurker said

    #8 Kevin O’Neill

    Kevin, I presume that you are referring to Cavalieri et al in your comment. Actually, if you look at Jeff’s graph of the difference between F17 and F13, the break seems to occur at about day 168 (no missing data). At day 189 (begin missing data) it drops but then recovers at day 204 (at the end of missing data). You see similar dip and rebounds in Jeff’s graph in the 19 – 33 and 64 – 73 periods.

  13. Kevin O'Neill said

    LL, yes i was looking at Cavalieri et al. I was also looking at the break at 190. I didn’t pay much attention to 168 because I was looking at the wrong graph in Cavalieri – (a) extent instead of (b) area. The extent graph does show a similar magnitude decrease at 170.

  14. Jeff condon said

    The NSIDC record returns 100% of all data for both satellites in 2007-2008.

    I’m wondering now what they did. No time tonight – same old mantra.

  15. n-g said

    The time axes for Figs. 3 and 4 seem to be offset relative to each other by about 0.04 years. In Fig. 3, the N17 curve first goes above the N13 curve at about 2007.08. In Fig. 4, the positive difference first appears at about 2007.04. Similar offsets apply for the other two major positive difference periods.

  16. Jeff Condon said

    n-g,

    I don’t know how you can see that but there appears to be a few pixel difference in the plots. It is the same data so I’m lost.

  17. Doug proctor said

    But any mistakes increase the trend, right?

    Benefit of the doubt consistently goes to the sea-loss side. Unconscious or by belief?

  18. Jeff Condon said

    Doug,

    I’m online now and don’t understand.

    If you are saying that most mistakes are to the benefit of AGW, we have found the opposite to be dramatically true. If you are saying that in this case it may go the other way, I’m not too sure. There is a real trend difference that appears due to a step but may not be. Just because we have identified a step in this single year, doesn’t mean that they don’t exist all over the record. I find it interesting that such a small and otherwise undetectable step, can create a trend equal to 40% of ice loss.

  19. [...] I have quietly been spending some more hours of free time on sea ice and have a ton of data to present.  Previously, we found that there was a nearly statistically significant difference in trend between satellite NOAA 17 and 13 that would create a spurious trend increase.  In addition, we found a fairly strong offset in the data which would tend to go in the opposite direction.  Post here. [...]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
Follow

Get every new post delivered to your Inbox.

Join 134 other followers

%d bloggers like this: