Support the Arctic Sea Ice Forum and Blog

Author Topic: Southern Hemisphere Sea Ice grid error?  (Read 9181 times)

ChadGreene

  • New ice
  • Posts: 8
    • View Profile
    • chadagreene.com
  • Liked: 0
  • Likes Given: 0
Southern Hemisphere Sea Ice grid error?
« on: September 15, 2015, 10:20:35 PM »
It looks like the NSIDC's georeferencing grids for sea ice data are offset by 400 km in the polar stereographic x direction.  In Matlab you can download and import a day's gridded sea ice concentration grid and corresponding lat,lon grids like this:

urlwrite('ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/pss25lons_v3.dat','pss25lons_v3.dat');
urlwrite('ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/pss25lats_v3.dat','pss25lats_v3.dat');
urlwrite('ftp://sidads.colorado.edu/pub/DATASETS/nsidc0081_nrt_nasateam_seaice/south/nt_20150914_f17_nrt_s.bin','test.bin');

fid = fopen('pss25lons_v3.dat');
lon = fread(fid,[316 332],'long','ieee-le')/100000;
fclose(fid);

fid = fopen('pss25lats_v3.dat');
lat = fread(fid,[316 332],'long','ieee-le')/100000;
fclose(fid);

fid = fopen('test.bin');
ci = fread(fid, [316 332]);
fclose(fid);


Using Antarctic Mapping Tools (http://www.mathworks.com/matlabcentral/fileexchange/47638), the grid looks okay until you put it in context with any other dataset such as a Bedmap2 Grounding line:

antmap
pcolorm(lat,lon,ci)
bedmap2 'grounding line'




To be clear, this is not an issue with the polar stereographic projection.  It is a problem with the lats and lons given in the pss25lats_v3.dat and pss25lons_v3.dat files.  This can be shown by plotting raw data in unprojected coordinates. 

In the sea ice concentration grid, values of 253 correspond to the coastline and 254 is the ice sheet.  Below I'm plotting each grid point associated with the ice sheet as black dots, and each grid point associated with the coast line as blue dots.  I'm also labeling McMurdo Station which is located at (166.67E,77.85S).  McMurdo is an identifiable location on any map because it sits on the coast and it was chosen by R.F. Scott because it is the southernmost point a ship can sail.

% Get indices of coastline:
coastline = ci==253;

% Get indices of land:
icesheet = ci==254;

figure
plot(lon(icesheet),lat(icesheet),'k.')
hold on
plot(lon(coastline),lat(coastline),'b.')

% Mark and label McMurdo Station:
plot(166.67,-77.85,'ro','linewidth',4)
text(166.67,-77.85,'MCM','vert','top',...
    'horiz','center','fontweight','bold','color','red')

xlabel 'longitude'
ylabel 'latitude'




McMurdo should sit directly on the coast.  To make it a bit more clear, we can double the grid for continuity and then plot the same as above.  Here's how I double the grid:

lat = [lat lat];
lon = [lon lon+360];
ci = [ci ci];




I've explored this a bit and it looks like the offset is 400 km in the polar stereographic easting direction.  It looks like the values are incorrect here as well: http://nsidc.org/data/polar_stereo/ps_grids.html.  I contacted NSIDC about the issue and they offered a proof that their grids are correct (http://nbviewer.ipython.org/gist/hwilcox/7e533c42abd38975a83a), but that proof does not address the issue that their grid incorrectly places McMurdo Station in the western hemisphere.

If I am missing something, please let me know.  Otherwise it may be best to manually georeference your sea ice data by

[x,y] = meshgrid(-3550000+12500:25000:4350000-12500,4350000-12500:-25000:-3950000+12500);
[lat,lon] = ps2ll(x,y,'TrueLat',-70,'EarthRadius',6378273,'Eccentricity',0.081816153);


and set

ci = flipud(rot90(ci));




Neven

  • Administrator
  • First-year ice
  • Posts: 9503
    • View Profile
    • Arctic Sea Ice Blog
  • Liked: 1336
  • Likes Given: 618
Re: Southern Hemisphere Sea Ice grid error?
« Reply #1 on: September 16, 2015, 02:15:11 PM »
Welcome Chad, and nice post. I've approved your membership, so you should be able to comment freely now.
The enemy is within
Don't confuse me with him

E. Smith

Wipneus

  • Citizen scientist
  • Young ice
  • Posts: 4220
    • View Profile
    • Arctische Pinguin
  • Liked: 1025
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #2 on: September 16, 2015, 07:07:13 PM »
...
I've explored this a bit and it looks like the offset is 400 km in the polar stereographic easting direction.  It looks like the values are incorrect here as well: http://nsidc.org/data/polar_stereo/ps_grids.html.  I contacted NSIDC about the issue and they offered a proof that their grids are correct (http://nbviewer.ipython.org/gist/hwilcox/7e533c42abd38975a83a), but that proof does not address the issue that their grid incorrectly places McMurdo Station in the western hemisphere.

...

If I understand correctly you say that the lon/lat values from the NSIDC files may perhaps be correct, but then the ice maps cannot be. If you where right, that would include maps from many other sources as well, not all (maybe not any) would have used these files.
About your claim that McMurdo in the map is in the western hemisphere. Attached is an ice map on which I placed an orange dot for McMurdo. I hope I got that about right.

The orange dot is on the right side of the image, that IS the eastern hemisphere.

400km is 16 25km cell widths. That is exactly the difference between the width and the height of the antarctic maps, can it be that some coordinate mix-up has occurred somewhere? 

ChadGreene

  • New ice
  • Posts: 8
    • View Profile
    • chadagreene.com
  • Liked: 0
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #3 on: September 16, 2015, 09:00:34 PM »
If I understand correctly you say that the lon/lat values from the NSIDC files may perhaps be correct...

I suspect the lat/lon values from the NSIDC files are *not* correct.  They are obtained by unprojecting a polar stereographic grid using X values that are exactly 400 kilometers offset.  Mapping each pixel in the sea ice grid to corresponding lats and lons given by pss25lats_v3.dat and pss25lons_v3.dat places the southernmost bit of ocean (McMurdo Sound) 95 kilometers into the western hemisphere. However, in the real world the southernmost bit of ocean is found 305 kilometers into the eastern hemisphere.   

Below I plot a day's sea ice concentration grid three ways in Matlab:

1. Directly map each pixel to pss25lats_v3 and pss25lons_v3.
2. Build the correct grid in polar stereographic coordinates and then transform to the correct lats and lons.
3. Recreate the incorrect pss25lats_v3 and pss25lons_v3 values by subtracting 400,000 meters from the correct X values, then perform inverse polar stereographic coordinate transformation.

% Download NSIDC's grid and a daily field of sea ice concentration:
urlwrite('ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/pss25lons_v3.dat','pss25lons_v3.dat');
urlwrite('ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/pss25lats_v3.dat','pss25lats_v3.dat');
urlwrite('ftp://sidads.colorado.edu/pub/DATASETS/nsidc0081_nrt_nasateam_seaice/south/nt_20150914_f17_nrt_s.bin','test.bin');

% Import lats:
fid = fopen('pss25lons_v3.dat');
pss25lons_v3 = fread(fid,[316 332],'long','ieee-le')/100000;
fclose(fid);

% Import lons:
fid = fopen('pss25lats_v3.dat');
pss25lats_v3 = fread(fid,[316 332],'long','ieee-le')/100000;
fclose(fid);

% Import a daily field of sea ice concentration data:
fid = fopen('test.bin');
ci = fread(fid, [316 332]);
fclose(fid);

% Plotting pss25lats in context with any other datasets reveals discrepancy: 
subplot(1,3,1)
antmap
pcolorm(pss25lats_v3,pss25lons_v3,ci)
bedmap2 'coast'

% We can build a correct grid in polar stereographic cartesian coordinates:
correctx = -3550000+12500:25000:4350000-12500;
y = 4350000-12500:-25000:-3950000+12500;
[Y,correctX] = meshgrid(y,correctx);

% Convert correct x,y to correct lat and lon:
[correctlat,correctlon] = ps2ll(correctX,Y,'TrueLat',-70,'EarthRadius',6378273,'Eccentricity',0.081816153);

% The correct grid lines up nicely with other datasets:
subplot(1,3,2)
antmap
pcolorm(correctlat,correctlon,ci)
bedmap2 'coast'

% NSIDC's values are off by 400,000 meters: 
nsidcX = correctX - 400000;

% NSIDC's pss25lats_v3.dat and pss25lons_v3.dat can be obtained like this: 
[nsidclat,nsidclon] = ps2ll(nsidcX,Y,'TrueLat',-70,'EarthRadius',6378273,'Eccentricity',0.081816153);

% The lat/lons derived from incorrect X values are the same as the grid
% provided in NSIDC's pss25lats_v3.dat and pss25lons_v3:   
subplot(1,3,3)
antmap
pcolorm(nsidclat,nsidclon,ci)
bedmap2 'coast'



« Last Edit: September 16, 2015, 09:05:52 PM by ChadGreene »

Wipneus

  • Citizen scientist
  • Young ice
  • Posts: 4220
    • View Profile
    • Arctische Pinguin
  • Liked: 1025
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #4 on: September 17, 2015, 08:41:36 AM »

I suspect the lat/lon values from the NSIDC files are *not* correct.  They are obtained by unprojecting a polar stereographic grid using X values that are exactly 400 kilometers offset. 

NSIDC already proved this not to be true. You did not object to this but instead complained about something else.

Quote
Mapping each pixel in the sea ice grid to corresponding lats and lons given by pss25lats_v3.dat and pss25lons_v3.dat places the southernmost bit of ocean (McMurdo Sound) 95 kilometers into the western hemisphere. However, in the real world the southernmost bit of ocean is found 305 kilometers into the eastern hemisphere.   

I am trying to add the other part of the proof you seem to require before looking at other possibilities.
Look at the image that I attached to my post above. McMurdo clearly lies in the right half of the image, that is the eastern hemisphere. You did not dispute this.
The pss25lons_v3.dat data if I read it into the correct array IS symmetrical, that is the values on the left are negative (lons are mapped -180 to +180), and on the right the same values but positive. So McMurdo, when plotted by these lons must ly in the correct hemisphere.
I conclude there is nothing in this file that supports your 16 cells shift.

I am not offering help in debugging a program in an unfamiliar language, but I am curious:
 
Quote
% Import a daily field of sea ice concentration data:
fid = fopen('test.bin');
ci = fread(fid, [316 332]);
fclose(fid);

How is matlab figuring out that the data in this file is preceded by a 300 byte binary header?
http://nsidc.org/data/docs/daac/nsidc0081_ssmi_nrt_seaice.gd.html


ChadGreene

  • New ice
  • Posts: 8
    • View Profile
    • chadagreene.com
  • Liked: 0
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #5 on: September 17, 2015, 07:06:19 PM »
Hi Wipneus,

Thanks for putting time and thought into this.  It is possible that I have made a mistake somewhere, so I appreciate your efforts to help me track down where I'm messing up.  If I am misunderstanding the sea ice grids I certainly want to know what I'm missing.

I do not dispute that the figure you posted looks correct, but my argument is that the sea ice concentration grids look good enough when they are plotted alone, but errors become evident when the georeferenced sea ice grids are plotted in context with other datasets.  How did you place the dot at McMurdo--by specifying geo coordinates?  How did you register the locations of each pixel of the sea ice concentration grid?  Can you add a graticule and maybe a coastline from a different dataset to provide context? 

What I keep seeing is that if I map each pixel in a sea ice concentration grid to corresponding lats and lons provided in pss25lats_v3.dat and pss25lons_v3.dat, things just don't line up.  I tried to show that as simply as possible, meaning without any coordinate transformations or fancy mapping programs, in Figure 3 of my original post.  Figure 3 of my original post shows dots marking the center locations of each pixel, and if you count them it's 16 pixels off.  When I plot the grids in projected coordinates the 16 pixel offset is manifest as a 400 km error in the polar stereographic easting direction.

Thanks for the note about fread.  That's another discrepancy I found between the documentation and the data files.  In Matlab,

ci = fread(fid, [316 332]);

is shorthand for

ci = fread(fid, [316 332],'uint8=>double');

and both of the above work well for importing sea ice data.  But specifying the 300 byte header skip by

ci = fread(fid, [316 332],'uint8=>double',300);

causes errors.  Unfortunately, it's a binary file so we can't just open it and read it directly.  That is a whole different complaint of mine--Especially for these small 316x332 grids there is no excuse for encoding the data in binary format.  And the data grids are not even in the same format as the corresponding lat, lon grids!  Why make the data so inaccessible?  Why perpetuate the notion that climate data are locked away in an ivory tower for only the eyes of elite experts to see?  Why make it so difficult to debug little georeferencing problems?  What poor design. 
« Last Edit: September 17, 2015, 08:42:22 PM by ChadGreene »

Wipneus

  • Citizen scientist
  • Young ice
  • Posts: 4220
    • View Profile
    • Arctische Pinguin
  • Liked: 1025
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #6 on: September 18, 2015, 10:47:52 AM »
Hi Wipneus,

Thanks for putting time and thought into this.  It is possible that I have made a mistake somewhere, so I appreciate your efforts to help me track down where I'm messing up.  If I am misunderstanding the sea ice grids I certainly want to know what I'm missing.

I am using the SIC data files myself, including  the cell area files, but not the lon/lat files. So the possibility of something wrong in those exist. After reading that NSIDC double checked it and found nothing wrong makes that a very small possibility IMO.
That leaves the SIC data files to consider, my first post dealt with that: very unlikely.

In all my judgement is that there is very likely something wrong in your reasoning, unfortunately most of that reasoning is expressed in matlab. That makes it rather difficult to follow and my knowledge of matlab is insufficient to point out  where the error lies.

Quote
I do not dispute that the figure you posted looks correct, but my argument is that the sea ice concentration grids look good enough when they are plotted alone, but errors become evident when the georeferenced sea ice grids are plotted in context with other datasets.  How did you place the dot at McMurdo--by specifying geo coordinates?  How did you register the locations of each pixel of the sea ice concentration grid?  Can you add a graticule and maybe a coastline from a different dataset to provide context? 

No I placed the orange dot manually after consulting a map on wikipedia.
As you agreed that the place is more or less correct, the argument goes as follows:
1) McM lies on the right half side of the map;
2) direct inspection of the longitude values in the pss25lons_v3 array shows that all values in the right half are from 0-180o

Conclusion must be that McM lies in the eastern hemisphere. Values could still be incorrect but not like what you found.

Now you might have a look at part 2 of my reasoning. Do you see the same thing? , if not we can talk about why not. A look at this figure from NSIDC documentation may helpful:


Quote
Thanks for the note about fread.  That's another discrepancy I found between the documentation and the data files.  In Matlab,

ci = fread(fid, [316 332]);

is shorthand for

ci = fread(fid, [316 332],'uint8=>double');

and both of the above work well for importing sea ice data.  But specifying the 300 byte header skip by

ci = fread(fid, [316 332],'uint8=>double',300);

causes errors.  Unfortunately, it's a binary file so we can't just open it and read it directly.  That is a whole different complaint of mine--Especially for these small 316x332 grids there is no excuse for encoding the data in binary format.  And the data grids are not even in the same format as the corresponding lat, lon grids!  Why make the data so inaccessible?  Why perpetuate the notion that climate data are locked away in an ivory tower for only the eyes of elite experts to see?  Why make it so difficult to debug little georeferencing problems?  What poor design.

The 300 bytes is correct  I can tell you from experience.
Matlab can do an fseek() call, otherwise a dummy read of 300 bytes should do the trick.

About 300 as argument four in fread, the matlab docs say:

Quote
A = fread(fileID, sizeA, precision, skip) skips skip bytes after reading each value

"after", perhaps it skips 300 bytes after each element, so 316x332 times?

Wipneus

  • Citizen scientist
  • Young ice
  • Posts: 4220
    • View Profile
    • Arctische Pinguin
  • Liked: 1025
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #7 on: September 18, 2015, 12:46:22 PM »
OK, I could not resist. Here is the output of your code (third figure) from the top post, with an added:

fseek(fid,300,0)

before reading in the SIC data.


Wipneus

  • Citizen scientist
  • Young ice
  • Posts: 4220
    • View Profile
    • Arctische Pinguin
  • Liked: 1025
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #8 on: September 18, 2015, 01:03:31 PM »
So by treating the 300 byte header as data, and reading 316 byte lines of data, byte 17 of each line is now read at the beginning of the lines. That effectively shifts the bulk of the image 400 km to the left and 25 km up.

ChadGreene

  • New ice
  • Posts: 8
    • View Profile
    • chadagreene.com
  • Liked: 0
  • Likes Given: 0
Re: Southern Hemisphere Sea Ice grid error?
« Reply #9 on: September 18, 2015, 06:27:06 PM »
Wipneus, that's it! The

fseek(fid,300,0)

bit gets everything back in order.  After including the fseek line I also verified that

[x,y] = ll2ps(lat,lon,'TrueLat',-70,'EarthRadius',6378273,'Eccentricity',0.081816152829245);

does indeed return the x and y values described on the http://nsidc.org/data/polar_stereo/ps_grids.html page.  Your solution works perfectly--thanks for setting me straight.

A-Team

  • Young ice
  • Posts: 2977
    • View Profile
  • Liked: 944
  • Likes Given: 35
Re: Southern Hemisphere Sea Ice grid error?
« Reply #10 on: September 23, 2015, 01:17:48 AM »
Quote
Unfortunately, it's a binary file so we can't just open it and read it directly.  That is a whole different complaint of mine--Especially for these small 316x332 grids there is no excuse for encoding the data in binary format.  And the data grids are not even in the same format as the corresponding lat, lon grids!  Why make the data so inaccessible?  Why perpetuate the notion that climate data are locked away in an ivory tower for only the eyes of elite experts to see?  Why make it so difficult to debug little georeferencing problems?  What poor design!


Exactly. Why?

99% of what they do could simply be dropped into small intuitive grayscale layers (ie standard raster GIS that the rest of the geo world uses) from which people could instantly port into their favorite format.

What are they trying to show with this file anyway, just ice/no ice/land? That's not even an 8-bit grayscale, just a 2-bit black and white with a 2-bit land mask  layer. At very very coarse resolution.

Then there's the equivalent csv (comma, separated, value) format universally importable into spreadsheets, even when gigantic. Anyone and everyone can open csv and port it elsewhere as needed.

I suppose their answer would be that their format is so wonderfully general that all types of data now and forever can be put into these massive complex all-encompassing formal data protocols so robots could potentially seamlessly crawl the internet to grab and assemble all manner of climate data and run it to update models. Which is just one coarse 316x332 grid after another.

To me, it is a huge mistake. I've seen it all before -- opaque master formatting protocols, $2500 software with proprietary command line work environment, grandiose data mining schemes -- and way too many people waste scads of time trying to get out of the stupid format into something commonsensical but mostly are shut out from participation altogether.

I am going to say the problem probably resides in the DC of NSIDC. At some point, probably prior to cheap cloud storage, they applied for a huge grant to be a data center. The proposal obviously could not just talk about providing a big hard drive, there had to be some elaborate chatter about advanced data formats.

Somehow forcing round pegs, square pegs, little pegs into a data protocol was supposed to revolutionize something. I'm not sure what -- I've tried their site a few times and just been put off by the sprawl of xml. Every time I've unpacked one of their netCDF files, the content vs padding ratio has been appalling.

It's way overdone -- and the idea of coercing everyone into the same ecosystem (matlab and its obtuse files) was a poor one. Worse, very few people are depositing data there (or anywhere).
« Last Edit: September 23, 2015, 09:49:13 AM by A-Team »