# -*- coding: utf-8 -*- """ .. :module:: backend :platform: Unix :synopsis: This represents the main engine for extracting from the raw datasets. .. moduleauthor:: Matthew Dillon """ import zipfile import gdal import numpy import itertools import os import errno from osgeo import osr import sqlite3 import shutil # Classes class SNAPDataSet: """ Tools to work with a SNAP Dataset. :param filename: A ZIP dataset from SNAP """ def __init__(self, filename): self.filename = filename self.zip_data = self.load_dataset() self.name_list = sorted(self.zip_data.namelist()) # At this point, all SNAP datasets are zipped in a directory self.zip_dir = self.name_list[0] # Prefix in the form xyz_mean_C_abc_model_ (month_year.tif following) self.prefix = self.name_list[1][len(self.zip_dir):-11] # Assume some info about dataset from the filename components = filename.replace('.', '_').split('_')[:-1] if 'historical' in components: # HISTORICAL DATA self.model = 'CRU' self.scenario = components[components.index('CRU') + 1] self.resolution = components[components.index('CRU') - 1] else: # PROJECTION DATA for comp in components: if comp == 'AK': startmarker = components.index(comp) + 2 if comp.startswith('sres'): self.scenario = comp.replace('sres', '').upper() endmarker = components.index(comp) self.model = "_".join(components[startmarker:endmarker]) self.resolution = components[startmarker - 1] def load_dataset(self): """ Import ZIP dataset. :returns: A reference to a zipfile """ return zipfile.ZipFile(self.filename,"r") def dump_raw_temperatures(self, communities, extracted_temps, out): """ Given a set of extracted temperatures, generate csv output of data. :param communities: Python list of community names :param extracted_temps: Numpy array with extracted temps :param out: path to output directory """ min_year = numpy.min(extracted_temps['year']) max_year = numpy.max(extracted_temps['year']) time_years = max_year - min_year + 1 i = 0 for community in communities: community = community.decode('utf-8') community = community.replace(" ", "_") outdir = ''.join([out, '/', community]) mkdir_p(outdir) outfile = ''.join([outdir, '/', community, '_', self.model, '_', self.scenario, '_', str(min_year), '_', str(max_year),'.txt']) header = ' '.join([community.replace("_", " "), ',', str(min_year), '-', str(max_year), '\nAverage Monthly' \ ' Air Temperature (deg C)\nYear, ' \ 'Jan, Feb, Mar, Apr, May, Jun, Jul, ' \ 'Aug, Sep, Oct, Nov, Dec']) temp_data = extracted_temps[i, :]['temperature'].reshape(time_years, 12) file_data = numpy.zeros((time_years, 13)) file_data[:, 1:] = temp_data file_data[:, 0] = numpy.arange(min_year, max_year+1) numpy.savetxt(outfile, file_data, fmt=('%d', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f', '%7.1f'), delimiter=',', header=header) i += 1 class GeoRefData(SNAPDataSet): """ Use GDAL to work with the SNAP datasets. :param SNAPDataSet: A SNAPDataSet object """ def __init__(self, filename): SNAPDataSet.__init__(self, filename) test_tiff = self.read_geotiff_as_gdal(1, int(self.filename[-13:-9])) self.cols = test_tiff.RasterXSize self.rows = test_tiff.RasterYSize self.bands = test_tiff.RasterCount geotransform = test_tiff.GetGeoTransform() self.origin_x = geotransform[0] self.origin_y = geotransform[3] self.pixel_width = geotransform[1] self.pixel_height = geotransform[5] # Close the file test_tiff = None def read_geotiff_as_gdal(self, month, year): """ Read GeoTIFF Data in from ZIP dataset. :param month: desired month (1- or 2-digit integer) :param year: desired year (4-digit integer) :returns: A GDAL data object """ # A bit clunky, but here we assemble a SNAP-style geotiff filename tiff = ''.join(['/vsizip/', self.filename, '/', self.zip_dir, self.prefix, str(month).zfill(2), '_',str(year), '.tif']) gdal_data = gdal.Open(tiff) return gdal_data def read_geotiff_as_array(self, month, year): """ Read GeoTIFF Data in from ZIP dataset. :param month: desired month (1- or 2-digit integer) :param year: desired year (4-digit integer) :returns: A Numpy array """ gdal_data = self.read_geotiff_as_gdal(month, year) temp_band = gdal_data.GetRasterBand(1) temp_data = temp_band.ReadAsArray(0, 0, self.cols, self.rows) temp_band = None gdal_data = None return temp_data def ne_to_indices(self, northing, easting): """ Convert Northings and Eastings (NAD 83 Alaska Albers Equal Area Conic) to X,Y array indices. :param northing: position northing (in meters) :param easting: position easting (in meters) :returns: array indices that correspond to location """ x_ind = (easting - self.origin_x)/self.pixel_width y_ind = (northing - self.origin_y)/self.pixel_height x_ind = x_ind.astype(numpy.int, copy=False) y_ind = y_ind.astype(numpy.int, copy=False) return (x_ind, y_ind) def indices_to_ne(self, x_ind, y_ind): """ Convert index values to Northings and Eastings (NAD 83 Alaska Albers Equal Area Conic). :param x_ind: array x-index :param y_ind: array y-index :returns: position northings and eastings (in meters) corresponding to location """ northing = self.origin_y + (y_ind * self.pixel_height) easting = self.origin_x + (x_ind * self.pixel_width) return (northing, easting) def extract_points(self, northing, easting, start_year, end_year): """ Extract points from range of years between start and end at the specified points (Jan->Dec). Point locations should be numpy arrays. :param northing: position northing (in meters) :param easting: position easting (in meters) :param start_year: 4-digit year for start of analysis period :param end_year: 4-digit year for end of analysis period, same as start_year if only analyzing one year :returns: numpy array of extracted temperatures """ x_offsets, y_offsets = self.ne_to_indices(northing, easting) years = list(range(start_year, end_year + 1)) months = list(range(1, 13)) # Record structure: (Year, Month, Temperature) # Each row represents a community, each column is a monthly temp. extracted_temps = numpy.zeros((len(x_offsets), 12*len(years)), dtype={'names': ['year', 'month', 'temperature'], 'formats':['i4', 'i4', 'f4']}) i = 0 for year, month in itertools.product(years, months): #for year in years: # for month in months: temp_data = self.read_geotiff_as_array(month, year) # gdal rotates for some reason, so y,x extracted_temps[:, i]['temperature'] = temp_data[y_offsets, x_offsets] extracted_temps[:, i]['year'] = year extracted_temps[:, i]['month'] = month i += 1 return extracted_temps # Functions def mkdir_p(path): """ Function to emulate mkdir -p functionality. Pulled from: http://stackoverflow.com/q/600268/313548 :param path: path to create new directory at :returns: creates a path at the desired location, if one does not already exist """ try: os.makedirs(path) except OSError as exc: if exc.errno == errno.EEXIST and os.path.isdir(path): pass else: raise def wgs84_to_ne(latitude, longitude): """ Convert WGS84 lat/long to Northings and Eastings (NAD 83 Alaska Albers Equal Area Conic) :param latitude: WGS84 latitude (in decimal degrees) :param longitude: WGS84 longitude (in decimal degrees) :returns: transformed coordinates to Alaska Albers """ wgspoint = osr.SpatialReference() wgspoint.ImportFromEPSG(4326) nepoint = osr.SpatialReference() nepoint.ImportFromEPSG(3338) transform = osr.CoordinateTransformation(wgspoint, nepoint) return transform.TransformPoint(longitude, latitude) def ne_to_wgs(northing, easting): """ Convert Northings and Eastings (NAD 83 Alaska Albers Equal Area Conic) to WGS84 lat/long . :param northing: AK Albers in meters :param easting: AK Albers in meters :returns: transformed coordinates in WGS84 lat long """ wgspoint = osr.SpatialReference() wgspoint.ImportFromEPSG(4326) nepoint = osr.SpatialReference() nepoint.ImportFromEPSG(3338) transform = osr.CoordinateTransformation(nepoint, wgspoint) return transform.TransformPoint(easting, northing) if __name__ == '__main__': print("nothing to see here...")