Commit a7859ff2 authored by Hoek, Steven's avatar Hoek, Steven
Browse files

Initial commit

parents
from .asciigrid import AsciiGrid;
from .floatingpointraster import FloatingPointRaster;
from .gridenvelope2d import GridEnvelope2D;
from .netcdf4envelope2d import Netcdf4Envelope2D;
from .hdf5raster import Hdf5Raster;
from .netcdf4raster import Netcdf4Raster;
\ No newline at end of file
import const;
import os.path;
import array;
from raster import Raster
from .gridenvelope2d import GridEnvelope2D;
class AsciiGrid(GridEnvelope2D, Raster):
"A raster represented by an ASCII file, with extension 'asc'"
# Constants
const.FILEXT = "asc";
const.MAXDIGITSPERCELL = 11;
# Data attributes - assign some dummy values for the mean time
name = "dummy." + const.FILEXT;
folder = os.getcwd();
"""
ncols = 1;
nrows = 1;
xll = 0.0;
yll = 0.0;
"""
cellsize = 100.0;
nodatavalue = -9999.0;
datatype = const.FLOAT;
dataformat='f'
# Private attributes
datafile = None;
currow = 0;
__mode = 'r';
__digitspercell = 7;
def __init__(self, filepath, *datatype):
# Initialise
DATAFILEXT = const.FILEXT
HEADEREXT = ""
# Retrieve the name from the filepath and assign - incl. extension
self.name = os.path.basename(filepath);
# Also derive the folder
self.folder = os.path.dirname(filepath);
# Finally set the datatype
if len(datatype) > 0:
if (datatype[0] == const.INTEGER):
self.datatype = const.INTEGER;
self.dataformat = 'i'
else:
self.datatype = const.FLOAT;
def open(self, mode, ncols=1, nrows=1, xll=0.0, yll=0.0, cellsize=100.0, nodatavalue=-9999.0):
# Initialise
super(AsciiGrid, self).open(mode);
# If file does not exist and mode[0] = 'w', create it!
if (mode[0] == 'w') and (not os.path.exists(self.folder + os.path.sep + self.name)):
self.datafile = file(self.folder + os.path.sep + self.name, 'w');
self.__mode = mode;
GridEnvelope2D.__init__(self, ncols, nrows, xll, yll, cellsize, cellsize);
self.cellsize = cellsize;
self.nodatavalue = nodatavalue;
self.writeheader();
return True;
else:
# Open the file
if os.path.exists(self.folder + os.path.sep + self.name):
self.datafile = open(self.folder + os.path.sep + self.name, mode[0]);
if (mode[0] == 'w'):
# Assign the data attributes
self.ncols = ncols;
self.nrows = nrows;
self.xll = xll;
self.yll = yll;
self.cellsize = cellsize;
self.nodatavalue = nodatavalue;
self.writeheader();
else:
# File is open - retrieve the data attributes from the header of the file
self.readheader();
# Also find out how many digits per cell were used - assume it's constant
pos = self.datafile.tell();
line = self.datafile.readline();
self.__digitspercell = ((1 + len(line)) / self.ncols) - 1;
self.datafile.seek(pos); # return to first line with data
GridEnvelope2D.__init__(self, self.ncols, self.nrows, self.xll, self.yll, self.cellsize, self.cellsize);
return True;
else: return False;
def readheader(self):
# Assume that the file is open; read header of the file and assign all attributes
if (self.datafile != None):
if (not self.datafile.closed):
hl = self.datafile.readline();
self.ncols = int(hl.replace('ncols', ''));
hl = self.datafile.readline();
self.nrows = int(hl.replace('nrows', ''));
hl = self.datafile.readline();
self.xll = float(hl.replace('xllcorner', ''));
hl = self.datafile.readline();
self.yll = float(hl.replace('yllcorner', ''));
hl = self.datafile.readline();
self.cellsize = float(hl.replace('cellsize', ''));
hl = self.datafile.readline();
if (self.datatype == const.INTEGER):
self.nodatavalue = int(hl.replace('NODATA_value', ''));
else:
self.nodatavalue = float(hl.replace('NODATA_value', ''));
else:
msg = "File " + self.name + " not found in folder " + self.folder;
raise IOError(msg);
def next(self, parseLine=True):
# Read the next row if possible, otherwise generate StopIteration
# Assume that the header lines have been read and are correct wrt. ncols and nrows
result = None;
try:
if (self.datafile != None):
if (not self.datafile.closed):
self.currow += 1;
if (self.currow > self.nrows):
raise StopIteration("Attempt to move beyond last row.");
# Allocate a new array with ncols of the right type
if (self.datatype == const.INTEGER):
result = array.array('l', self.ncols * [self.nodatavalue]);
else:
result = array.array('f', self.ncols * [self.nodatavalue]);
# Now fill the array - first translate whitespace into space
rawline = self.datafile.readline();
if parseLine:
i = 0;
for x in rawline.split():
if (i < self.ncols):
if (self.datatype == const.INTEGER):
result[i] = int(x);
else:
result[i] = float(x);
i = i + 1;
return result;
else: raise StopIteration("Attempt to read raster data from a closed file.");
else: raise StopIteration("Attempt to read raster data from an unassigned file.")
except Exception, e:
print "Error: " + str(e);
raise StopIteration;
"""
def hasSameExtent(self, obj):
if not isinstance(obj, AsciiGrid):
return False;
if self.ncols != obj.ncols:
return False;
if self.nrows != obj.nrows:
return False;
if abs(self.cellsize - obj.cellsize) > const.epsilon:
return False;
if abs(self.xll - obj.xll) > const.epsilon:
return False;
if abs(self.yll - obj.yll) > const.epsilon:
return False;
else:
return True;
"""
@staticmethod
def getFileExt():
return const.FILEXT;
def writeheader(self):
# Assume that the file is open; write header of the file with all attributes
if (self.datafile != None):
if (not self.datafile.closed):
try:
self.datafile.write("ncols " + str(self.ncols).rjust(const.MAXDIGITSPERCELL + 1) + "\n");
self.datafile.write("nrows " + str(self.nrows).rjust(const.MAXDIGITSPERCELL + 1) + "\n");
self.datafile.write("xllcorner " + str(self.xll).rjust(const.MAXDIGITSPERCELL + 1) + "\n");
self.datafile.write("yllcorner " + str(self.yll).rjust(const.MAXDIGITSPERCELL + 1) + "\n");
self.datafile.write("cellsize " + str(self.cellsize).rjust(const.MAXDIGITSPERCELL + 1) + "\n");
self.datafile.write("NODATA_value " + str(self.nodatavalue).rjust(const.MAXDIGITSPERCELL + 1) + "\n");
except Exception, e:
print e;
msg = "Header lines could not be written to file " + self.name + " in folder " + self.folder;
raise IOError(msg);
def writenext(self, sequence_with_data):
# Write the next line if possible, otherwise generate StopIteration
# We assume that exactly 1 row is included.
try:
if (self.datatype == const.INTEGER):
for k in range(0, self.ncols):
s = str(sequence_with_data[k]).rjust(const.MAXDIGITSPERCELL + 1);
self.datafile.write(s);
else:
totalwidth = const.MAXDIGITSPERCELL - 1
fmtstr = "{:" + str(totalwidth) + ".2f}"
for k in range(0, self.ncols):
s = fmtstr.format(sequence_with_data[k]).rjust(const.MAXDIGITSPERCELL + 1);
self.datafile.write(s);
return self.datafile.write("\n");
except Exception, e:
print e;
raise StopIteration
def flush(self):
self.datafile.flush();
def reset(self):
self.datafile.seek(0);
if (self.__mode[0] == 'r'):
self.readheader();
super(AsciiGrid, self).reset()
\ No newline at end of file
# Copyright (c) 2004-2015 Alterra, Wageningen-UR
# Steven Hoek (steven.hoek@wur.nl), June 2015
from raster import Raster
from gridenvelope2d import GridEnvelope2D
import os.path
from math import fabs
# Abstract class
class BandRaster(GridEnvelope2D, Raster):
# Constants
INTEL = 'I'; # Least significant byte first
UNSIGNEDINT = 'UNSIGNEDINT'
# Data attributes
nbands = 3 #default
nbits = 8 #default
nodatavalue = 256; #default
dataformat = "h" # Default data format (2-byte signed short int)
def __init__(self, filepath, *dataformat):
Raster.__init__(self, filepath)
# Retrieve the name from the filepath and assign - incl. extension
self.name = os.path.basename(filepath);
# Also derive the folder
self.folder = os.path.dirname(filepath);
# Finally set the dataformat
# TODO: get this from the header file?
if len(dataformat) > 0:
self.dataformat = dataformat[0]
def open(self, mode, ncols=1, nrows=1, nbands=1, xll=0.0, yll=0.0, cellsize=1.0, nodatavalue=-9999.0):
super(BandRaster, self).open(mode);
# If file does not exist and mode[0] = 'w', create it!
if (mode[0] == 'w'):
# Set data attributes and write header file anyhow
self.nbands = nbands;
self.cellsize = cellsize
self.nodatavalue = nodatavalue
self.envelope = GridEnvelope2D.__init__(self, ncols, nrows, xll, yll, cellsize, cellsize)
self.writeheader()
if (not os.path.exists(self.folder + os.path.sep + self.name)):
self.datafile = file(self.folder + os.path.sep + self.name, 'w') # File does not exist
else:
self.datafile = open(self.folder + os.path.sep + self.name, mode[0] + 'b') # File exists
return True
else:
# Open the file
if os.path.exists(os.path.join(self.folder, self.name)):
# Open the file and retrieve the data attributes from the header file
self.datafile = open(self.folder + os.path.sep + self.name, mode[0] + 'b');
self.readheader();
self.xll = self.xul;
if self.ycoords_sort == 'DESC':
self.yll = self.yul - self.nrows * self.dy;
else:
self.yll = self.yul + self.nrows * self.dy;
self.envelope = GridEnvelope2D.__init__(self, self.ncols, self.nrows, self.xll, self.yll, self.dx, self.dy)
return True
else: return False
def readheader(self):
# See if we can find the header file to use. Then read it and assign all attributes
hdrFilename = ""
if os.path.isfile(os.path.join(self.folder, self.name + "." + self.HEADEREXT)): #@UndefinedVariable
hdrFilename = os.path.join(self.folder, self.name + "." + self.HEADEREXT) #@UndefinedVariable
else:
name_wo_ext = os.path.splitext(os.path.join(self.folder, self.name))[0]
if os.path.isfile(name_wo_ext + "." + self.HEADEREXT ): #@UndefinedVariable
hdrFilename = os.path.join(name_wo_ext + "." + self.HEADEREXT) #@UndefinedVariable
if hdrFilename == "":
raise ValueError("Not sure about name of header file: " + hdrFilename + "?");
hdrFilename = os.path.normpath(hdrFilename)
if not os.path.exists(hdrFilename):
raise ValueError("Header file not found: " + hdrFilename);
hf = open(hdrFilename, 'r');
if (hf != None):
hl = hf.readline();
self.byteorder = str(hl.replace('BYTEORDER', '').strip());
if self.byteorder != self.INTEL:
raise ValueError("Unsupported byte order")
hl = hf.readline();
layout = str(hl.replace('LAYOUT', '').strip());
extuc = self.DATAFILEXT.upper()
if layout != extuc:
raise ValueError("Incorrect layout in header - apparently not a " + extuc + " file");
hl = hf.readline();
self.nrows = int(hl.replace('NROWS', ''));
hl = hf.readline();
self.ncols = int(hl.replace('NCOLS', ''));
hl = hf.readline();
self.nbands = int(hl.replace('NBANDS', ''));
hl = hf.readline();
self.nbits = int(hl.replace('NBITS', ''));
hl = hf.readline();
bandrowbytes = int(hl.replace('BANDROWBYTES', ''));
if self.nbits * self.ncols / 8 != bandrowbytes:
raise ValueError("Incorrect number of bytes per band row in header");
hl = hf.readline();
totalrowbytes = int(hl.replace('TOTALROWBYTES', ''));
if self.nbands * bandrowbytes != totalrowbytes:
raise ValueError("Incorrect total bytes per row in header");
hl = hf.readline();
if hl.find('PIXELTYPE') != -1:
# Assume ESRI style header
self.pixeltype = str(hl.replace('PIXELTYPE', '').strip());
hl = hf.readline();
ulxmap = float(hl.replace('ULXMAP', ''));
hl = hf.readline();
ulymap = float(hl.replace('ULYMAP', ''));
hl = hf.readline();
self.dx = float(hl.replace('XDIM', ''));
hl = hf.readline();
self.dy = float(hl.replace('YDIM', ''));
# ulxmap - The x-axis map coordinate of the center of the upper left pixel
# ulymap - The y-axis map coordinate of the center of the upper left pixel
self.xul = ulxmap - 0.5*self.dx
self.yul = ulymap + 0.5*self.dy
hl = hf.readline();
self.nodatavalue = int(hl.replace('NODATA', ''));
hf.close()
else:
# Assume that the rest of the information has to be read from a world file
name_wo_ext = os.path.splitext(os.path.join(self.folder, self.name))[0]
if os.path.isfile(name_wo_ext + "." + self.WORLDEXT):
# TODO: Adapt the following so that it accounts also for rotated mapsheets
sign = lambda x: (1, -1)[x<0];
hdrFilename = os.path.normpath(os.path.join(name_wo_ext + "." + self.WORLDEXT))
hf = open(hdrFilename, 'r')
hl = hf.readline();
self.dx = float(hl.strip());
hl = hf.readline();
self.roty = float(hl.strip());
hl = hf.readline();
self.rotx = float(hl.strip());
eps = 0.0001;
if abs(self.rotx)>eps or abs(self.roty)>eps:
raise NotImplementedError("Cannot handle rotated mapsheets yet!")
hl = hf.readline();
self.dy = fabs(float(hl.strip()));
if sign(float(hl.strip())) == 1.0: self.ycoords_sort = 'ASC';
hl = hf.readline();
self.xul = float(hl.strip()) - 0.5 * self.dx;
hl = hf.readline();
self.yul = float(hl.strip()) + 0.5 * self.dy;
hf.close();
self.cellsize = (self.dx + self.dy) / 2
else:
msg = "Unable to open header file " + hdrFilename + " in folder " + self.folder;
raise IOError(msg);
def writeheader(self):
# TODO: test!
# Write header file with all attributes
hdrFilename = os.path.join(self.folder, self.name + "." + self.HEADEREXT) #@UndefinedVariable
hdrFilename = os.path.normpath(hdrFilename)
try:
# Open the file if it exists, otherwise create it
if os.path.exists(hdrFilename):
hf = open(hdrFilename, 'w');
else:
hf = file(hdrFilename, 'w');
# Now write all the attributes
hf.write("BYTEORDER " + str(self.byteorder) + "\n")
hf.write("LAYOUT " + self.DATAFILEXT.upper() + "\n")
hf.write("NROWS " + str(self.nrows) + "\n")
hf.write("NCOLS " + str(self.ncols) + "\n")
hf.write("NBANDS " + str(self.nbands) + "\n")
hf.write("NBITS " + str(self.nbits) + "\n")
bandrowbytes = self.nbits * self.ncols / 8
hf.write("BANDROWBYTES " + str(bandrowbytes) + "\n")
hf.write("TOTALROWBYTES " + str(self.nbands * bandrowbytes) + "\n")
if self.dataformat.lower() == 'i':
if self.dataformat.isupper(): hf.write("PIXELTYPE UNSIGNEDINT\n")
else: hf.write("PIXELTYPE SIGNEDINT\n")
else: hf.write("PIXELTYPE FLOAT\n")
# ulxmap - The x-axis map coordinate of the center of the upper left pixel
ulxmap = self.xll + 0.5*self.dx
hf.write("ULXMAP " + str(ulxmap) + "\n");
# ulymap - The y-axis map coordinate of the center of the upper left pixel
ulymap = self.yll + self.nrows*self.dy - 0.5*self.dy # Check this!!!
hf.write("ULYMAP " + str(ulymap) + "\n");
hf.write("XDIM " + str(self.dx) + "\n");
hf.write("YDIM " + str(self.dy) + "\n");
hf.write("NODATA " + str(self.nodatavalue) + "\n");
except Exception, e:
msg = "Header file " + hdrFilename + " could not be written in folder " + self.folder;
raise IOError(msg + "(" + str(e) + ")");
def _is_sequence(self, arg):
return (not hasattr(arg, "strip") and hasattr(arg, "__getitem__") or hasattr(arg, "__iter__"))
def next(self, parseLine=True):
pass
def writenext(self, sequence_with_data):
# input is sequence type - e.g. list, array.array or numpy.array
pass
@staticmethod
def getDataFileExt(self):
return self.DATAFILEXT;
@staticmethod
def getHeaderFileExt():
return self.HEADEREXT
def close(self):
if self.datafile:
if not self.datafile.closed:
self.datafile.close();
def reset(self):
self.currow = 0;
self.datafile.seek(0);
def __exit__(self, exception_type, exception_value, traceback):
self.close()
def __enter__(self):
return self
\ No newline at end of file
from raster import Raster
import os
from libtiff.libtiff_ctypes import TIFFFieldInfo, TIFFDataType, FIELD_CUSTOM, add_tags
class BaseTiffRaster(Raster):
# Constants
DATAFILEXT = "tif";
HEADEREXT = "tfw"; # WORLD_EXT?
__BYTESPERCELL = 4;
# Data attributes - assign some dummy values for the mean time
name = "dummy.tif";
folder = os.getcwd();
nodatavalue = -9999.0;
byteorder = 'II'; # Little endian
roty = 0.0;
rotx = 0.0;
pageNumber = [0, 1]
modelPixelScale = [0.0, 0.0, 0.0]
modelTiepoint = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
ModelTransformation = 16 * [0.0]
GeoKeyDirectory = 32 * [0]
GeoDoubleParams = 1.0e-317
GeoAsciiParams = "GCS_WGS_1984|"
GdalMetadata = '<GDALMetadata>\n <Item name="RepresentationType" sample="0">ATHEMATIC</Item>\n</GDALMetadata>'
def __init__(self, filepath, *datatype):
pass
def open(self, mode, ncols=1, nrows=1, xll=0.0, yll=0.0, cellsize=1.0, nodatavalue=-9999.0):
super(BaseTiffRaster, self).open(mode)
extra_tags = [
TIFFFieldInfo(297, 2, 2, TIFFDataType.TIFF_SHORT, FIELD_CUSTOM, True, False, "PageNumber"),
TIFFFieldInfo(33550, 3, 3, TIFFDataType.TIFF_DOUBLE, FIELD_CUSTOM, True, False, "ModelPixelScaleTag"),
TIFFFieldInfo(33922, 6, 6, TIFFDataType.TIFF_DOUBLE, FIELD_CUSTOM, True, False, "ModelTiepointTag"),
TIFFFieldInfo(34264, 16, 16, TIFFDataType.TIFF_DOUBLE, FIELD_CUSTOM, True, False, "ModelTransformationTag"),
TIFFFieldInfo(34735, 32, 32, TIFFDataType.TIFF_SHORT, FIELD_CUSTOM, True, False, "GeoKeyDirectoryTag"),
TIFFFieldInfo(34736, -1, -1, TIFFDataType.TIFF_DOUBLE, FIELD_CUSTOM, True, False, "GeoDoubleParamsTag"),
TIFFFieldInfo(34737, -1, -1, TIFFDataType.TIFF_ASCII, FIELD_CUSTOM, True, False, "GeoAsciiParamsTag"),
TIFFFieldInfo(42112, -1, -1, TIFFDataType.TIFF_ASCII, FIELD_CUSTOM, True, False, "GDAL_METADATA"),
TIFFFieldInfo(42113, -1, -1, TIFFDataType.TIFF_ASCII, FIELD_CUSTOM, True, False, "GDAL_NODATA")
]
add_tags(extra_tags)
def set_extra_tags(self):
# The following is in view of the georeferencing
self.datafile.SetField("PageNumber", self.pageNumber)
self.datafile.SetField("ModelPixelScaleTag", self.modelPixelScale)
self.datafile.SetField("ModelTiepointTag", self.modelTiepoint)
self.datafile.SetField("ModelTransformationTag", self.ModelTransformation)
self.datafile.SetField("GeoKeyDirectoryTag", self.GeoKeyDirectory)
#self.datafile.SetField("GeoDoubleParamsTag", self.GeoDoubleParams)
self.datafile.SetField("GeoAsciiParamsTag", str(self.GeoAsciiParams))
self.datafile.SetField("GDAL_METADATA", str(self.GdalMetadata))
def get_extra_tags(self):
# The following is in view of the georeferencing
self.pageNumber = self.datafile.GetField("PageNumber")
self.modelPixelScale = self.datafile.GetField("ModelPixelScaleTag")
self.modelTiepoint = self.datafile.GetField("ModelTiepointTag")
self.ModelTransformation = self.datafile.GetField("ModelTransformationTag")
self.GeoKeyDirectory = self.datafile.GetField("GeoKeyDirectoryTag")
#self.GeoDoubleParams = self.datafile.GetField("GeoDoubleParamsTag")
self.GeoAsciiParams = self.datafile.GetField("GeoAsciiParamsTag")
self.GdalMetadata = self.datafile.GetField("GDAL_METADATA")
def copy_extra_tags(self, rg):
# Check input
if not isinstance(rg, BaseTiffRaster):
raise TypeError("Input is not a raster!")
# Ok, go ahead
if rg.pageNumber != None: self.pageNumber = rg.pageNumber
self.modelPixelScale = rg.modelPixelScale
self.modelTiepoint = rg.modelTiepoint
if rg.ModelTransformation != None: self.ModelTransformation = rg.ModelTransformation
self.GeoKeyDirectory = rg.GeoKeyDirectory
#self.GeoDoubleParams = rg.GeoDoubleParams
self.GeoAsciiParams = rg.GeoAsciiParams
self.GdalMetadata = rg.GdalMetadata
\ No newline at end of file
# -*- coding: utf-8 -*-
# Copyright (c) 2004-2015 Alterra, Wageningen-UR
# Steven Hoek (steven.hoek@wur.nl), April 2015
import os.path
import stat
import struct
from array import array
from math import fabs
from bandraster import BandRaster
try:
import numpy
HAS_NUMPY = True
except ImportError:
HAS_NUMPY = False
class BilRaster(BandRaster):
"A raster represented by 2 files, with extensions 'bil' and 'bil.hdr'"
# TODO: change some function calls so that it's more obvious how to use "with" statements
# Use of such statements are essential when files are written to disk!
# Data attributes - assign some dummy values for the mean time
DATAFILEXT = 'bil'
HEADEREXT = "hdr"
name = "dummy.bil";
WORLDEXT = "blw"
def __init__(self, filepath, *dataformat):
# Initialise super class instance
BandRaster.__init__(self, filepath, dataformat[0])
self.byteorder = self.INTEL
self.pixeltype = self.UNSIGNEDINT
def open(self, mode, ncols=1, nrows=1, nbands=3, xll=0, yll=0, cellsize=100, nodatavalue=256):
self.nbits = struct.calcsize(self.dataformat)*8
result = super(BilRaster, self).open(mode, ncols, nrows, nbands, xll, yll, cellsize, nodatavalue);
if (mode[0] == 'w'):
return result
else:
# Open the file
if os.path.exists(os.path.join(self.folder, self.name)):
# Check given format string is valid
bytesperpix = 2 #default
try:
bytesperpix = struct.calcsize(self.dataformat)
except:
raise ValueError, "Supplied data format " + str(self.dataformat) + " is invalid"
# end try
# Check file size matches with size attributes
fileinfo = os.stat(os.path.join(self.folder, self.name))
filesize = fileinfo[stat.ST_SIZE]
if (filesize == 0) and (mode[0] == 'w'):
print "Empty BIL file found. I'm going to overwrite it ..."
else: