In [19]:
from charistools import hypsometry, modis_info
from netCDF4 import Dataset
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt

In [2]:
modis_info.area_per_500m_pixel_km2


Out[2]:
0.2146586737339221

In [3]:
help(hypsometry)


Help on module charistools.hypsometry in charistools:

NAME
    charistools.hypsometry

FILE
    /Users/brodzik/.conda/envs/pmesdr/lib/python2.7/site-packages/charistools/hypsometry.py

CLASSES
    Hypsometry
    
    class Hypsometry
     |  The Hypsometry class will read(write) CHARIS hypsometry data
     |  from(to) ASCII files.  The populated object will contain the
     |  comments from the beginning of the file, and a pandas
     |  DataFrame.
     |  
     |  2014-09-24 M. J. Brodzik brodzik@nsidc.org 303-492-8263
     |  National Snow & Ice Data Center, Boulder CO
     |  Copyright (C) 2014-2015 Regents of the University of Colorado at Boulder
     |  
     |  Methods defined here:
     |  
     |  __init__(self, filename=None, comments=[], data=Empty DataFrame
     |  Columns: []
     |  Index: [], verbose=False)
     |      from charistools import hypsometry
     |      
     |      hyp = hypsometry.Hypsometry(comments=[], data=pd.DataFrame(),
     |                                  filename=None, verbose=False)
     |      
     |      Initializer for a CHARIS hypsometry object
     |      
     |      filename : filename to open and read from
     |      comments : list of strings
     |      data : pandas DataFrame with elevations in columns and dates in rows
     |             index can be [None] for undated data, or
     |             should be a pandas.tseries.index.DatetimeIndex for dated
     |             contents
     |  
     |  append(self, elevation, modice, contour_m=100.0, min_contour_m=None, max_contour_m=None, verbose=False)
     |      hypsometry.append(elevation, modice, contour_m=100.,
     |                        min_contour_m=None,
     |                        max_contour_m=None,
     |                        verbose=False)
     |      
     |      Returns hypsometry object, with one row of MODICE area by the requested
     |      contour levels.
     |      
     |      Parameters
     |          elevation : ndarray
     |              raster with elevations in meters
     |          modice : ndarray, same size as elevation
     |              modice=2 for landice
     |          contour_m : value, default 100. m
     |              contour interval in meters
     |          min_contour_m : value, default None
     |              bottom of lowest contour interval, in meters
     |              if None, determine it from elevation data
     |          max_contour_m : value, default None
     |              top of highest contour interval, in meters
     |              if None, determine it from elevation data
     |          verbose : Boolean, default=False
     |      
     |      Example:
     |      from charistools import hypsometry
     |      hyps = hypsometry.Hypsometry(comments=['Hypsometry for MODICE from file.'])
     |      # Set elevation and modice raster ndarrays
     |      # Assumes elevation array is in meters
     |      # Assumes modice value for landice is 2
     |      hyps.append( elevation, modice )
     |      
     |      Returns hypsometry object, with one row MODICE area by the requested
     |      contour levels.
     |  
     |  data_by_doy(self, verbose=False)
     |      from charistools import hypsometry
     |      hyps = hypsometry.Hypsometry(filename='filename.txt')
     |      doy_series = hyps.data_by_doy()
     |      
     |      Returns a Series object, with the hypsometry data, summed by row (doy).
     |  
     |  read(self, filename, verbose=False)
     |      from charistools import hypsometry
     |      hyps = hypsometry.Hypsometry()
     |      hyps.read( filename, verbose=False )
     |      
     |      Reads elevation by date hypsometry data from filename into a pandas
     |      DataFrame.
     |      
     |      Assumes file format:
     |      1) 0 or more comment lines, beginning with '#'
     |      2) XX number of elevation bands
     |      3) lower bounds of each elevation band
     |      4) Data records by date
     |  
     |  write(self, filename, decimal_places=6, verbose=False)
     |      Writes elevation by date hypsometry data to filename,
     |      formatted for use with either CHARIS IDL or python
     |      modelling software.  See hypsometry.read for file format
     |      description.
     |      
     |      from charistools import hypsometry hyps =
     |      hypsometry.Hypsometry(comments=['first','second'])
     |      hyps.write( filename, verbose=False )
     |      
     |      If object comments do not begin with '#', then they will
     |      be prepended by this character in the output file.
     |  
     |  ----------------------------------------------------------------------
     |  Data and other attributes defined here:
     |  
     |  comments = []
     |  
     |  data = Empty DataFrame
     |  Columns: []
     |  Index: []

DATA
    print_function = _Feature((2, 6, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0)...



In [7]:
hyps = hypsometry.Hypsometry(comments=['some comments'])

In [9]:
hyps.comments


Out[9]:
['some comments']

In [10]:
hyps.data


Out[10]:

In [11]:
%cd /Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc
%ls


/Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc
MODICE.v0.4.h22v04.1strike.min05yr.mask.nc
MODICE.v0.4.h22v04.2strike.min05yr.mask.nc
MODICE.v0.4.h22v04.3strike.min05yr.mask.nc
MODICE.v0.4.h22v05.1strike.min05yr.mask.nc
MODICE.v0.4.h22v05.2strike.min05yr.mask.nc
MODICE.v0.4.h22v05.3strike.min05yr.mask.nc
MODICE.v0.4.h23v04.1strike.min05yr.mask.nc
MODICE.v0.4.h23v04.2strike.min05yr.mask.nc
MODICE.v0.4.h23v04.3strike.min05yr.mask.nc
MODICE.v0.4.h23v05.1strike.min05yr.mask.nc
MODICE.v0.4.h23v05.2strike.min05yr.mask.nc
MODICE.v0.4.h23v05.3strike.min05yr.mask.nc
MODICE.v0.4.h23v06.1strike.min05yr.mask.nc
MODICE.v0.4.h23v06.2strike.min05yr.mask.nc
MODICE.v0.4.h23v06.3strike.min05yr.mask.nc
MODICE.v0.4.h24v04.1strike.min05yr.mask.nc
MODICE.v0.4.h24v04.2strike.min05yr.mask.nc
MODICE.v0.4.h24v04.3strike.min05yr.mask.nc
MODICE.v0.4.h24v05.1strike.min05yr.mask.nc
MODICE.v0.4.h24v05.2strike.min05yr.mask.nc
MODICE.v0.4.h24v05.3strike.min05yr.mask.nc
MODICE.v0.4.h24v06.1strike.min05yr.mask.nc
MODICE.v0.4.h24v06.2strike.min05yr.mask.nc
MODICE.v0.4.h24v06.3strike.min05yr.mask.nc
MODICE.v0.4.h25v05.1strike.min05yr.mask.nc
MODICE.v0.4.h25v05.2strike.min05yr.mask.nc
MODICE.v0.4.h25v05.3strike.min05yr.mask.nc
MODICE.v0.4.h25v06.1strike.min05yr.mask.nc
MODICE.v0.4.h25v06.2strike.min05yr.mask.nc
MODICE.v0.4.h25v06.3strike.min05yr.mask.nc
MODICE.v0.4.h26v05.1strike.min05yr.mask.nc
MODICE.v0.4.h26v05.2strike.min05yr.mask.nc
MODICE.v0.4.h26v05.3strike.min05yr.mask.nc
MODICE.v0.4.h26v06.1strike.min05yr.mask.nc
MODICE.v0.4.h26v06.2strike.min05yr.mask.nc
MODICE.v0.4.h26v06.3strike.min05yr.mask.nc

In [13]:
file = Dataset('MODICE.v0.4.h23v05.1strike.min05yr.mask.nc', mode='r', format='NETCDF4')

In [14]:
file


Out[14]:
<type 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format UNDEFINED):
    Title: MODICE mask for a minimum number of years.
    Institution: National Snow & Ice Data Center, Boulder, CO USA
    Source: MODICEv04
    History: Created on Tue Jun 16 14:28:19 2015 by combine_modice.pro
    Comment: Mask locations with 2 indicate MODICE for >= min_years
    MODIS_tile_id: h23v05
    Start_year: 2000
    End_year: 2014
    Min_years: 5
    MODICE_FILES: 2000/h23v05/MODICE.v0.4.h23v05.2000.1strike.landice.mask.tif, 2001/h23v05/MODICE.v0.4.h23v05.2001.1strike.landice.mask.tif, 2002/h23v05/MODICE.v0.4.h23v05.2002.1strike.landice.mask.tif, 2003/h23v05/MODICE.v0.4.h23v05.2003.1strike.landice.mask.tif, 2004/h23v05/MODICE.v0.4.h23v05.2004.1strike.landice.mask.tif, 2005/h23v05/MODICE.v0.4.h23v05.2005.1strike.landice.mask.tif, 2006/h23v05/MODICE.v0.4.h23v05.2006.1strike.landice.mask.tif, 2007/h23v05/MODICE.v0.4.h23v05.2007.1strike.landice.mask.tif, 2008/h23v05/MODICE.v0.4.h23v05.2008.1strike.landice.mask.tif, 2009/h23v05/MODICE.v0.4.h23v05.2009.1strike.landice.mask.tif, 2010/h23v05/MODICE.v0.4.h23v05.2010.1strike.landice.mask.tif, 2011/h23v05/MODICE.v0.4.h23v05.2011.1strike.landice.mask.tif, 2012/h23v05/MODICE.v0.4.h23v05.2012.1strike.landice.mask.tif, 2013/h23v05/MODICE.v0.4.h23v05.2013.1strike.landice.mask.tif, 2014/h23v05/MODICE.v0.4.h23v05.2014.1strike.landice.mask.tif
    dimensions(sizes): Columns(2400), Rows(2400)
    variables(dimensions): uint8 modice_min_year_mask(Rows,Columns), float32 latitude(Rows,Columns), float32 longitude(Rows,Columns)
    groups: 

In [15]:
modice = file.variables["modice_min_year_mask"][:]

In [17]:
np.shape(modice)


Out[17]:
(2400, 2400)

In [20]:
%cd /Users/brodzik/projects/CHARIS/elevation_data/SRTMGL3
%ls


/Users/brodzik/projects/CHARIS/elevation_data/SRTMGL3
#00notes.txt#
00notes.txt
00notes_retired.txt
CHARIS_DEM_Larissa.tfw
CHARIS_DEM_Larissa.tif
Koksu_Kunes.tfw
Koksu_Kunes.tif
Koksu_Kunes.zip
SRTMGL3.v0.1.h22v04.tif
SRTMGL3.v0.1.h22v05.tif
SRTMGL3.v0.1.h22v06.tif
SRTMGL3.v0.1.h23v04.tif
SRTMGL3.v0.1.h23v05.tif
SRTMGL3.v0.1.h23v06.tif
SRTMGL3.v0.1.h24v04.tif
SRTMGL3.v0.1.h24v05.tif
SRTMGL3.v0.1.h24v06.tif
SRTMGL3.v0.1.h25v04.tif
SRTMGL3.v0.1.h25v05.tif
SRTMGL3.v0.1.h25v06.tif
SRTMGL3.v0.1.h26v04.tif
SRTMGL3.v0.1.h26v05.tif
SRTMGL3.v0.1.h26v06.tif
SRTMGL3.v0.1.h27v04.tif
SRTMGL3.v0.1.h27v05.tif
SRTMGL3.v0.1.h27v06.tif
SRTMGL3_GCS_edited.tif
SRTMGL3_GCS_edited.tif.aux.xml
SRTMGL3_GCS_edited.zip
SRTMGL3_GCS_edited_langtang_area.tif
SRTMGL3_GCS_edited_langtang_area_utm45.tif
SRTMGL3_GCS_edited_uib_area.tif
SRTMGL3_sin.tfw
SRTMGL3_sin.tif
SRTMGL3_vakhsh.tif
du
srtmgl3_clip_by_basin.sh*
srtmgl3_clip_by_basin.sh~

In [22]:
elevation_file = 'SRTMGL3.v0.1.h23v05.tif'
dataset = gdal.Open(elevation_file, gdal.GA_ReadOnly)

In [24]:
dataset.RasterCount


Out[24]:
1

In [25]:
band = dataset.GetRasterBand(1)
elevation = band.ReadAsArray()

In [26]:
np.shape(elevation)


Out[26]:
(2400, 2400)

In [27]:
np.amin(elevation)


Out[27]:
-32768

In [28]:
np.amax(elevation)


Out[28]:
7710

In [29]:
hyps


Out[29]:
<charistools.hypsometry.Hypsometry instance at 0x10752cb00>

In [30]:
hyps.comments


Out[30]:
['some comments']

In [31]:
hyps.append(elevation,modice,min_contour_m=1400,verbose=True)


charistools.hypsometry : begin.
Out[31]:
()

In [32]:
hyps.data


Out[32]:
1400.0 1500.0 1600.0 1700.0 1800.0 1900.0 2000.0 2100.0 2200.0 2300.0 ... 6900.0 7000.0 7100.0 7200.0 7300.0 7400.0 7500.0 7600.0 7700.0 7800.0
Date
None 0 0 0 0 0 0 0 0 0 0 ... 20.607233 14.167472 11.806227 8.586347 4.507832 3.005221 2.575904 0.214659 0.214659 0

1 rows × 65 columns


In [33]:
hypsh23v05 = hyps

In [35]:
elevation_file_h24 = 'SRTMGL3.v0.1.h24v05.tif'
dataset24 = gdal.Open(elevation_file_h24, gdal.GA_ReadOnly)
band24 = dataset24.GetRasterBand(1)
elevation24 = band24.ReadAsArray()

In [36]:
np.amin(elevation24)


Out[36]:
111

In [37]:
np.amax(elevation24)


Out[37]:
8451

In [38]:
%cd /Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc
%ls


/Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc
MODICE.v0.4.h22v04.1strike.min05yr.mask.nc
MODICE.v0.4.h22v04.2strike.min05yr.mask.nc
MODICE.v0.4.h22v04.3strike.min05yr.mask.nc
MODICE.v0.4.h22v05.1strike.min05yr.mask.nc
MODICE.v0.4.h22v05.2strike.min05yr.mask.nc
MODICE.v0.4.h22v05.3strike.min05yr.mask.nc
MODICE.v0.4.h23v04.1strike.min05yr.mask.nc
MODICE.v0.4.h23v04.2strike.min05yr.mask.nc
MODICE.v0.4.h23v04.3strike.min05yr.mask.nc
MODICE.v0.4.h23v05.1strike.min05yr.mask.nc
MODICE.v0.4.h23v05.2strike.min05yr.mask.nc
MODICE.v0.4.h23v05.3strike.min05yr.mask.nc
MODICE.v0.4.h23v06.1strike.min05yr.mask.nc
MODICE.v0.4.h23v06.2strike.min05yr.mask.nc
MODICE.v0.4.h23v06.3strike.min05yr.mask.nc
MODICE.v0.4.h24v04.1strike.min05yr.mask.nc
MODICE.v0.4.h24v04.2strike.min05yr.mask.nc
MODICE.v0.4.h24v04.3strike.min05yr.mask.nc
MODICE.v0.4.h24v05.1strike.min05yr.mask.nc
MODICE.v0.4.h24v05.2strike.min05yr.mask.nc
MODICE.v0.4.h24v05.3strike.min05yr.mask.nc
MODICE.v0.4.h24v06.1strike.min05yr.mask.nc
MODICE.v0.4.h24v06.2strike.min05yr.mask.nc
MODICE.v0.4.h24v06.3strike.min05yr.mask.nc
MODICE.v0.4.h25v05.1strike.min05yr.mask.nc
MODICE.v0.4.h25v05.2strike.min05yr.mask.nc
MODICE.v0.4.h25v05.3strike.min05yr.mask.nc
MODICE.v0.4.h25v06.1strike.min05yr.mask.nc
MODICE.v0.4.h25v06.2strike.min05yr.mask.nc
MODICE.v0.4.h25v06.3strike.min05yr.mask.nc
MODICE.v0.4.h26v05.1strike.min05yr.mask.nc
MODICE.v0.4.h26v05.2strike.min05yr.mask.nc
MODICE.v0.4.h26v05.3strike.min05yr.mask.nc
MODICE.v0.4.h26v06.1strike.min05yr.mask.nc
MODICE.v0.4.h26v06.2strike.min05yr.mask.nc
MODICE.v0.4.h26v06.3strike.min05yr.mask.nc

In [39]:
file24 = Dataset('MODICE.v0.4.h24v05.1strike.min05yr.mask.nc', mode='r', format='NETCDF4')

In [40]:
modice24 = file24.variables["modice_min_year_mask"][:]

In [41]:
print(np.amin(modice24),np.amax(modice24))


(0, 2)

In [42]:
hypsh24v05 = hypsometry.Hypsometry(comments=['h24v05'])

In [43]:
hypsh24v05.append(elevation24,modice24,min_contour_m=1400,verbose=True)


charistools.hypsometry : begin.
Out[43]:
()

In [44]:
hypsh23v05.data


Out[44]:
1400.0 1500.0 1600.0 1700.0 1800.0 1900.0 2000.0 2100.0 2200.0 2300.0 ... 6900.0 7000.0 7100.0 7200.0 7300.0 7400.0 7500.0 7600.0 7700.0 7800.0
Date
None 0 0 0 0 0 0 0 0 0 0 ... 20.607233 14.167472 11.806227 8.586347 4.507832 3.005221 2.575904 0.214659 0.214659 0

1 rows × 65 columns


In [45]:
hypsh24v05.data


Out[45]:
1400.0 1500.0 1600.0 1700.0 1800.0 1900.0 2000.0 2100.0 2200.0 2300.0 ... 7600.0 7700.0 7800.0 7900.0 8000.0 8100.0 8200.0 8300.0 8400.0 8500.0
Date
None 0 0 0 0 0 0 0 0 0 0 ... 4.078515 2.575904 1.073293 0.643976 0.214659 0.214659 0.214659 0.214659 0.214659 0

1 rows × 72 columns


In [46]:
all = hypsh23v05.data + hypsh24v05.data

In [47]:
all


Out[47]:
1400.0 1500.0 1600.0 1700.0 1800.0 1900.0 2000.0 2100.0 2200.0 2300.0 ... 7600.0 7700.0 7800.0 7900.0 8000.0 8100.0 8200.0 8300.0 8400.0 8500.0
Date
None 0 0 0 0 0 0 0 0 0 0 ... 4.293173 2.790563 1.073293 NaN NaN NaN NaN NaN NaN NaN

1 rows × 72 columns


In [50]:
orig = hypsometry.Hypsometry(filename='/Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/IN_Hunza_at_Danyour.0100m.modicev04_1strike_area_by_elev.txt')

In [51]:
orig.data


Out[51]:
1400. 1500. 1600. 1700. 1800. 1900. 2000. 2100. 2200. 2300. ... 6800. 6900. 7000. 7100. 7200. 7300. 7400. 7500. 7600. 7700.
Date
NaN 0 0 0 0 0 0 0 0 0 0 ... 18.68 16.31 10.95 11.16 8.37 4.51 3.22 3.22 1.29 0.21

1 rows × 64 columns


In [ ]: