<jupyter_start><jupyter_text>Extract a section from the NSIDC sea ice clim download from here: https://nsidc.org/data/NSIDC-0192/versions/3<jupyter_code>import numpy as np
import matplotlib.pylab as plt
# https://gis.stackexchange.com/questions/97828/reading-zipped-esri-bil-files-with-python
def ReadBilFile(bil):
    import gdal
    gdal.GetDriverByName('EHdr').Register()
    img = gdal.Open(bil)
    band = img.GetRasterBand(1)
    data = band.ReadAsArray()
    return data<jupyter_output><empty_output><jupyter_text>rename file to "mean.sep.1979-2020.s.bil" , create header file below (https://nsidc.org/support/22256950-How-do-I-import-the-sea-ice-trends-data-into-ArcGIS-)<jupyter_code>%%bash

cat > mean.sep.1979-2020.s.hdr <<END
nrows 332

ncols 316

nbands 1

nbits 8

layout bil

byteorder I

ulxmap -3937500

ulymap 4337500

xdim 25000

ydim 25000
END
ls mean.sep*bil mean.sep*hdr
a = ReadBilFile('mean.sep.1979-2020.s.bil')
# keep the sea ice fraction only
a_masked=np.ma.masked_greater(a, 100., copy=True) 
# quick look
plt.imshow(a_masked, vmin=0., vmax=100.)
plt.colorbar()<jupyter_output><empty_output><jupyter_text>find the indices for a section in the bitmap with NSIDCs python tools<jupyter_code>!git clone https://github.com/nsidc/polar_stereo.git<jupyter_output>Klone nach 'polar_stereo' ...
remote: Enumerating objects: 139, done.[K
remote: Total 139 (delta 0), reused 0 (delta 0), pack-reused 139[K
Empfange Objekte: 100% (139/139), 29.47 MiB | 7.48 MiB/s, Fertig.
Löse Unterschiede auf: 100% (74/74), Fertig.<jupyter_text>cd polar_stereo/source/  nsidc_polar_lonlat should do the job, use nsidc_polar_ij as a check<jupyter_code>from polar_convert import polar_xy_to_lonlat
import numpy as np


def nsidc_polar_ij(i, j, grid, hemisphere):
    """Transform from NSIDC Polar Stereographic I, J coordinates
    to longitude and latitude coordinates
    
    Args:
        i (int): an integer or integer array giving the x grid coordinate(s)
        j (int): an integer or integer array giving the y grid coordinate(s)
        grid (float): 6.25, 12.5 or 25; the grid cell dimensions in km
        hemisphere (1 or -1): Northern or Southern hemisphere
    
    Returns:
        If i and j are scalars then the result is a
        two-element list containing [longitude, latitude].
        If i and j are numpy arrays then the result will be a two-element
        list where the first element is a numpy array containing
        the longitudes and the second element is a numpy array containing
        the latitudes.

    Examples:
        print(nsidc_polar_ij(608, 896, 12.5, 1))
            [350.01450147320855, 34.40871032516291]
    """

    true_scale_lat = 70
    re = 6378.273
    e = 0.081816153

    if grid != 6.25 and grid != 12.5 and grid != 25:
        raise ValueError("Legal grid values are 6.25, 12.5, or 25")
    
    if hemisphere != 1 and hemisphere != -1:
        raise ValueError("Legal hemisphere values are 1 or -1")

    if hemisphere == 1:
        delta = 45
        imax = 1216
        jmax = 1792
        xmin = -3850 + grid/2
        ymin = -5350 + grid/2
    else:
        delta = 0
        imax = 1264
        jmax = 1328
        xmin = -3950 + grid/2
        ymin = -3950 + grid/2

    if grid == 12.5:
        imax = imax//2
        jmax = jmax//2
    elif grid == 25:
        imax = imax//4
        jmax = jmax//4

    if np.any(np.less(i, 1)) or np.any(np.greater(i, imax)):
        raise ValueError("'i' value is out of range: [1, " + str(imax) + "]")
    if np.any(np.less(j, 1)) or np.any(np.greater(j, jmax)):
        raise ValueError("'j' value is out of range: [1, " + str(jmax) + "]")

    # Convert I, J pairs to x and y distances from origin.
    x = ((i - 1)*grid) + xmin
    y = ((jmax - j)*grid) + ymin
    lonlat = polar_xy_to_lonlat(x, y, true_scale_lat, re, e, hemisphere)
    lon = lonlat[0] - delta
    lon = lon + np.less(lon, 0)*360
    return [lon, lonlat[1]]
import numpy as np
from polar_convert import polar_lonlat_to_xy


def nsidc_polar_lonlat(longitude, latitude, grid, hemisphere):
    """Transform from geodetic longitude and latitude coordinates
    to NSIDC Polar Stereographic I, J coordinates
    
    Args:
        longitude (float): longitude or longitude array in degrees
        latitude (float): latitude or latitude array in degrees (positive)
        grid (float): 6.25, 12.5 or 25; the grid cell dimensions in km
        hemisphere (1 or -1): Northern or Southern hemisphere
    
    Returns:
        If longitude and latitude are scalars then the result is a
        two-element list containing [I, J].
        If longitude and latitude are numpy arrays then the result will
        be a two-element list where the first element is a numpy array for
        the I coordinates and the second element is a numpy array for
        the J coordinates.

    Examples:
        print(nsidc_polar_lonlat(350.0, 34.41, 12.5, 1))
            [608, 896]
    """

    true_scale_lat = 70
    re = 6378.273
    e = 0.081816153

    if grid != 6.25 and grid != 12.5 and grid != 25:
        raise ValueError("Legal grid value are 6.25, 12.5, or 25")
    
    if hemisphere >= 0:
        delta = 45
        imax = 1216
        jmax = 1792
        xmin = -3850 + grid/2
        ymin = -5350 + grid/2
    else:
        delta = 0
        imax = 1264
        jmax = 1328
        xmin = -3950 + grid/2
        ymin = -3950 + grid/2

    if grid == 12.5:
        imax = imax//2
        jmax = jmax//2
    elif grid == 25:
        imax = imax//4
        jmax = jmax//4

    xy = polar_lonlat_to_xy(longitude + delta, np.abs(latitude),
                            true_scale_lat, re, e, hemisphere)
    i = (np.round((xy[0] - xmin)/grid)).astype(int) + 1
    j = (np.round((xy[1] - ymin)/grid)).astype(int) + 1
    # Flip grid orientation in the 'y' direction
    j = jmax - j + 1
    return [i, j]<jupyter_output><empty_output><jupyter_text>Define a section here:<jupyter_code>NSIDC_lats=np.arange(50.,90.,1.)  # 60°S to 90°S, in 0.5 steps
NSIDC_lons=-30.*np.ones_like(NSIDC_lats) # section along 30°W
hemisphere=-1 # Southern Hemisphere
grid=25 # res is 25km
# find the right pixels
NSIDC_section_idx=nsidc_polar_lonlat(longitude=NSIDC_lons, latitude=NSIDC_lats, grid=grid, hemisphere=hemisphere)
NSIDC_section_idx
# transform back, looks more or less ok
nsidc_polar_ij(i=NSIDC_section_idx[0], j=NSIDC_section_idx[1], grid=grid, hemisphere=hemisphere)<jupyter_output><empty_output><jupyter_text>check where the section lies<jupyter_code>NSIDC_section_idx=np.array(NSIDC_section_idx)
NSIDC_section_idx=NSIDC_section_idx[[1,0],:] # switch i and j
NSIDC_section_idx
# keep the sea ice fraction only
a_masked=np.ma.masked_greater(a, 100., copy=True) 
a_masked[tuple(NSIDC_section_idx)]=200
# this is the correct section at 30W
plt.imshow(a_masked, vmin=0., vmax=100.)
plt.colorbar()<jupyter_output><empty_output><jupyter_text>plot the section<jupyter_code># keep the sea ice fraction only
a_masked=np.ma.masked_greater(a, 100., copy=True) 
section30w=a_masked[tuple(NSIDC_section_idx)]
plt.plot(-NSIDC_lats, section30w)
def determine_extent_lat(zz_area, lats):
    
    maxlat=np.nanmax(lats[zz_area[:]>0.15])
    
    xp=[maxlat+1.,maxlat] # maxlat+1. is equal or below 0.15
    fp=[zz_area[lats==maxlat+1].data[0], zz_area[lats==maxlat].data[0]]

    return np.interp(0.15, fp, xp)
determine_extent_lat(section30w/100., -NSIDC_lats)
from matplotlib import cm
import matplotlib.colors as mcol
import numpy.matlib
from matplotlib.backends.backend_pdf import PdfPages
fig, ax = plt.subplots(figsize=(14,1))

lvIce = np.linspace(0.35,1.0,17) #12
cmIce = cm.bone(lvIce)
newCmap = mcol.ListedColormap(cmIce)

newCmap.set_under('0.2')

# set sizes
numfontsize=33
texfontsize=37

which='NSIDC'
NSIDCextent=determine_extent_lat(section30w/100., -NSIDC_lats)

X=-NSIDC_lats
Z=section30w

Y=np.arange(0,5)
Yrep=np.matlib.repmat(Y, np.shape(NSIDC_lats)[0], 1).T

im=plt.pcolor(np.array([X,X,X,X,X]),Yrep, np.array([Z,Z,Z,Z,Z]), cmap=newCmap, vmin=15., vmax=100., shading='nearest')

plt.vlines(NSIDCextent, 0., 4., color='cyan', linewidths=4, linestyle='dashed')
plt.axis([-78., 0., 0., 4.]);

# tick sizes
plt.xticks(fontsize=numfontsize, rotation=0);
plt.yticks(fontsize=numfontsize, rotation=0);

plt.xlabel('latitude',fontsize=texfontsize);
ax.axes.yaxis.set_visible(False)
ax.set_facecolor('0.2')

# annotation
plt.annotate(which, xy=(0.84, 0.22), xycoords='axes fraction', backgroundcolor = 'none', color='w', fontsize=33)

# Now adding the colorbar
cbaxes = fig.add_axes([.95, -.3, 0.027, 3.7]) 
cbar = plt.colorbar(im, orientation = 'vertical', cax = cbaxes, extend='min')  

cbar_label='sea ice [%]'
cbar.set_label(cbar_label, fontsize=texfontsize)
cbar.ax.tick_params(labelsize=numfontsize) 

#plt.tight_layout()

saveFIG=True
if saveFIG==True:
   pp = PdfPages('./Figure4a_ice.pdf')
   pp.savefig(fig,bbox_inches = 'tight')
   pp.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12150369210114119, 'fraction_numerical': 0.05728350861490266, 'mean_word_length': 3.3330101793504605, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '300708', 'license': 'MIT', 'n_tokens_mistral': 3735, 'n_tokens_neox': 3305, 'n_words': 852, 'path': 'notebooks/Figure4_NSIDCice.ipynb', 'repo_name': 'trackow/AntarcticSeaIce_NatCommun', 'size': '8938'}