'''
This script ingests a 24km unpacked IMS ascii data file and converts it to a packed form, ready for ingest to ArcGIS. 

Use of this data should be acknowledged and cited as follows:

National Ice Center. 2008, updated daily. IMS Daily Northern Hemisphere Snow and Ice Analysis at 1 km, 4 km, and 24 km Resolutions, Version 1. 
[Indicate subset used]. Boulder, Colorado USA. NSIDC: National Snow and Ice Data Center. doi: http://dx.doi.org/10.7265/N52R3PMC. [Date Accessed].

Last edited June 19, 2017 G. Deemer

===============================================
Technical Contact
===============================================

NSIDC User Services
National Snow and Ice Data Center
CIRES, 449 UCB
University of Colorado
Boulder, CO 80309-0449  USA
phone: +1 303.492.6199
fax: +1 303.492.2468
e-mail: nsidc@nsidc.org

*************************************************

'''
import numpy as np

FileName = "ims1998280_24km_v1.1"
with open(FileName + ".asc", 'r') as f:
    # Skip 1365 bytes within the header.
    f.seek(1365)
    # Take all lines after the header, strip the end-line character and join all in one string.
    Lines = " ".join(line.strip('\n') for line in f)

# Remove all white space from Lines with Lines.split.
# If characters (s) in string are digits, convert to integers, then convert to a NumPy Array.
# Lastly, reshape to 1024 x 1024 24 km resolution.
DataArray = np.reshape(np.asarray([int(s) for s in Lines.split() if s.isdigit()]), [1024,1024])

# From the documentation (Table 3): 
# For unpacked data, integer value of 164 is sea ice, while 165 is snow-covered land.
# Convert 164 to 3 (sea ice) and 165 to 4 (snow covered land) to align with packed data.
DataArray = np.where(DataArray == 164, 3, DataArray)
DataArray = np.where(DataArray == 165, 4, DataArray)

# Create the output header to be saved to the new file.
HeaderText = "NCOLS 1024\nNROWS 1024\nXLLCORNER -12126597.0\nYLLCORNER -12126840.0\nCELLSIZE 23684.997"
np.savetxt(FileName + "_2packed.asc",  DataArray, header = HeaderText, delimiter = " ", fmt = '%s', comments = '')

print("Done")