'''
High level geo file image access
Base GeoImage class is a wrapper around gdal providing easy access to
image metadata, file read/write, and i/o. Higher level classes (such
as DGImage) can build on this class to provide data specific metadata and
spectral transformations.
'''
from __future__ import division
from osgeo import gdal, gdalconst, osr, ogr
import numpy as np
import os
import warnings
import collections
import textwrap
import tempfile
import logging
import math
import platform
from collections import Sequence
import tinytools as tt
# package import
import constants as const
# Module setup
gdal.UseExceptions()
ogr.UseExceptions()
logger = logging.getLogger(__name__)
# To get access to logging statements from the command line:
# import logging
# logging.basicConfig(level=logging.DEBUG) # or your desired level
class OverlapError(ValueError):
'''Raise when the window does not overlap the image. This can be
caught and passed when the window is expected to not overlap in
some cases.
'''
pass
class GeoImage(object):
"""
Base image class providing high-level access to image data and metadata
as well as methods to easily interact with the image and related
vector files.
Functionality is tested against .TIL, .VRT, and .TIF formats, but this
base class should support any GDAL format.
Parameters
----------
file_in : str
String describing a file on disk that is of a valid input format.
derived_dir : str
The location to store files created by the class.
Attributes
----------
files : tinytools.bunch.OrderedBunch
A collection of the image files used in the object. The GeoImage
class populates:
-derived_dir : The directory into which create files are stored.
-dfile : The input image file.
-dfile_tiles : The tiles of the virtual data set - if the data set
doesn't contain tiles, this should be equal to dfile.
Inherited classes can override/extend this listing.
meta : tinytools.bunch.OrderedBunch
Summary metadata of the base image.
shape : tuple
shape of the image in gdal format (bands,x,y).
resolutions : tuple
length 2 tuple with resolutions of x and y image dimensions.
"""
def __init__(self, file_in, derived_dir=None):
"""Initialize class with data and meta-data from file. __init__
class is in the class definition. """
## Search for files that are needed
if not os.path.isfile(file_in):
raise ValueError("The file that was passed in does not exist.")
## Start populating files dictionary for geoimage info
# ... variables created here:
# self.files_dict['fin']
# self.files_dict['gdal_file_list']
# self.files_dict['dfile']
# self.files_dict['dfile_tiles']
# self.files_dict['derived_dir']
# Get the file name and full file directory
ifile = os.path.abspath(file_in)
fname = os.path.basename(ifile)
fdir = os.path.dirname(ifile)
flist = os.listdir(fdir)
# Create files dictionary to populate - this will be bunched later
#!# self.files_dict = {}
self.files = tt.bunch.OrderedBunch({})
# Set the place to store/retrive derived files (i.e. spectral data)
if derived_dir:
if not os.path.isdir(derived_dir):
raise ValueError("The requested derived data directory does "
"not exist.")
if not os.access(os.path.join(derived_dir, ''), os.W_OK):
raise ValueError("Write access is required for the requested "
"location passed into derived_dir.")
self.files.derived_dir = os.path.join(derived_dir, '')
else:
if os.access(fdir, os.W_OK):
self.files.derived_dir = fdir
else:
self.files.derived_dir = fdir
warnings.warn("The input file location is not writable. "
"Derived file creation (i.e. spectral files) "
"will not be available. Either write permissions "
"need to be provided or the object can be "
"reinstantiated with a writable location passed "
"to the input variable dervied_store_dir.")
### Setup the dataset and subdataset variables
(tmpfile,tmptiles)=self._get_file_and_tiles(ifile)
#!# self.files_dict['dfile'] = tmpfile
#!# self.files_dict['dfile_tiles'] = tmptiles
self.files.dfile = tmpfile
self.files.dfile_tiles = tmptiles
# Create the files dictionary bunch
#!#self.files = tt.bunch.OrderedBunch(self.files_dict)
## Populate geoimage info
# ... variables created here:
# self._fobj
# self.meta_geoimg_dict
# self.meta
# Open the image in files.dfile
self._fobj = self._get_gdal_obj(self.files.dfile,
self.files.dfile_tiles)
# Populate metadata info from gdal
self._set_metadata()
def _get_file_and_tiles(self, ifile):
# If fname is a .til file then create .vrt
if tt.files.filter(ifile, '*.TIL', case_sensitive=False):
file_loc = ifile
tiles_loc = tt.pvl.read_from_pvl(ifile,'filename')
dname = os.path.dirname(ifile)
tiles_loc = [os.path.join(dname,x) for x in tiles_loc]
# If fname is a VRT, pull subdatasets from the file object
elif tt.files.filter(ifile, '*.VRT', case_sensitive=False):
file_loc = ifile
tmp = gdal.Open(file_loc) # Open to pull VRT file list
tmp_files = tmp.GetFileList()
tiles_loc = [x for x in tmp_files if not
tt.files.filter(x, '*.VRT', case_sensitive=False)]
tmp = None # Close the opened file from above
# If this is an ENVI file, then a file without an extension should
# exist and will be the root file.
elif os.path.isfile(os.path.splitext(ifile)[0]):
tmp = os.path.splitext(ifile)[0]
file_loc = tmp
tiles_loc = [tmp]
# If file input isn't a .TIL, .VRT, or ENVI file, just pass it into the
# object vars to pass to gdal.
else:
# Else, just pass it on to open in gdal
file_loc = ifile
tiles_loc = [ifile]
return (file_loc,tiles_loc)
def _get_gdal_obj(self, dfile, dfile_tiles):
'''Return gdal object for the GeoImage.'''
# Need to handle .TIL files specifically because gdal does not fully
# support them.
if tt.files.filter(dfile, '*.TIL', case_sensitive=False):
# If this is a .TIL file, then create a VRT, copy it to an
# "in memory" obj and then remove the VRT file from disk.
# Get a temporary file
file_temp = tempfile.NamedTemporaryFile(suffix=".VRT")
# Build the vrt
# If there is only one tile, use gdal_translate to ease the
# support of 1b files since gdalbuildvrt doesn't support non-
# projected files.
if len(dfile_tiles) == 1:
cmd = []
cmd.append("gdal_translate")
cmd.append("-of")
cmd.append("VRT")
cmd.append(dfile_tiles[0])
cmd.append(file_temp.name)
else:
cmd = []
cmd.append("gdalbuildvrt")
cmd.append(file_temp.name)
for i in dfile_tiles: cmd.append(i)
# Execute the buildvrt command and print the output via logging.
dump = tt.cmd_line.exec_cmd(cmd,ret_output=True)
logger.debug(dump)
del dump
# gdal does not issue errors to the command line,
# so try/except on tt.cmd_line.exec_cmd won't work...
# Check that file exists to see if the buildvrt succeeded.
if not os.path.isfile(file_temp.name):
raise StandardError("Creation of file " + file_temp.name + " "
"failed. This could possibly be a "
"write access problem?")
vvv = gdal.Open(file_temp.name)
# Create MEM copy of VRT - created "in memory" by passing an empty
# file name string to the VRT driver. If I used the "MEM" driver,
# the full dataset would be read into memory on creation.
drv = gdal.GetDriverByName('VRT')
obj = drv.CreateCopy('',vvv)
# Delete VRT
file_temp.close()
if os.path.isfile(file_temp.name):
raise StandardError("Removal of file " + file_temp.name + " "
"failed. There is something wrong with "
"the .TIL handling.")
else:
obj = gdal.Open(self.files.dfile, gdalconst.GA_ReadOnly)
return obj
def _set_metadata(self):
""" Get image metadata."""
meta_geoimg_dict = read_geo_file_info(self._fobj)
# Need to handle case of a .TIL that results an in memory VRT.
# In this case, file_name will be the VRT string when returned
# from the gdal driver above and does not print well or return the
# intended information.
if not os.path.isfile(meta_geoimg_dict['file_name']):
meta_geoimg_dict['file_name'] = self.files.dfile
# Add geoio class name to dictionary
meta_geoimg_dict['class_name'] = self.__class__.__name__
### OrderedBunch the metadata from the read_geo_file_info dictionary
self.meta = tt.bunch.OrderedBunch(meta_geoimg_dict)
# Set class members
self.shape = self.meta.shape
self.resolution = self.meta.resolution
def __repr__(self):
"""Human readable image summary similar to the R package 'raster'."""
sss = ''
su = self.meta
prefixes = collections.OrderedDict()
prefixes['Class Name'] = (['class_name'],'')
prefixes['Driver Name'] = (['driver_name'],'')
if 'product_level' in su:
prefixes['Product Level'] = (['product_level'],'')
prefixes['Data Type'] = (['gdal_dtype_name'],'')
prefixes['File Name'] = (['file_name'],'')
prefixes['File List'] = (['file_list'], '')
prefixes['Dimensions'] = (['shape'],
' (nlayers, nrows, ncols)')
prefixes['Resolution'] = (['resolution'],' (x,y)')
# The following is inverted because xstart, xend, etc are calculated
# in pixel space and this is inverted from the North = max, South = min
# paradigm.
prefixes['Extent'] = (['extent'],' (ul_x, ul_y, lr_x, lr_y)')
prefixes['Projection String'] = (['pprint_proj_string'],'')
prefixes['Geo Transform'] = (['geo_transform'],'')
if 'authority' in su:
prefixes['Authority'] = (['authority'], '')
if 'band_centers' in su:
prefixes['Band Centers (nm)'] = (['band_centers'],'')
### Loop through prefixes and su to print data to screen
# Gen max length of labels to set prefix length
prelen = max([len(x) for x in prefixes])
# Loop through each prefix to put together wrapped string
for x in prefixes:
prefix = x+' '*(prelen-len(x))+' : '
width_set = 80
wrapper = textwrap.TextWrapper(initial_indent=prefix,
width=width_set,
replace_whitespace=False,
subsequent_indent=' '*len(prefix))
message = ', '.join([str(su[y]) for y in prefixes[x][0]])
message = message+prefixes[x][1]
# Handle different message formats:
# If message contains new lines, just pass those through to print.
if message.find('\n') != -1:
sss = sss + prefix + message.replace('\n','\n'+' '*prelen)+'\n'
# Else if message is not empty, pass to wrapper.fill
elif message:
sss = sss + wrapper.fill(message) + '\n'
# Else just pass the empty message along with prefix
else:
sss = sss + prefix + '\n'
return sss
def print_img_summary(self):
"""Echo the object's __repr__ method."""
print(self.__repr__())
def __iter__(self):
'''Yield from default iter_window iterator.'''
for x in self.iter_window():
yield x
def iter_base(self, xoff, yoff, win_xsize, win_ysize, **kwargs):
'''
Base iterator function to yield data from array-like window parameters.
Parameters
----------
xoff : array_like
x offset(s) for the image regions to be read.
yoff : array_like
y offset(s) for the image regions to be read.
win_xsize : array_like
window x-dim size(s) for the image regions to be read.
win_ysize : array_like
window y-dim size(s) for the image regions to be read.
kwargs : optional
keyword arguments to be passed to get_data.
Yields
------
ndarray
Three dimensional numpy array of data from the requested region
of the image.
'''
logger.debug('*** begin iter_base ***')
# Broadcast array_like
windows = np.broadcast(xoff,yoff,win_xsize,win_ysize)
# Iterate through windows generated from input parameters
for w in windows:
logger.debug('window parameters: xoff %s, yoff %s, '
'win_xsize %s, win_ysize %s',
w[0], w[1], w[2], w[3])
yield self.get_data(window=w,**kwargs)
def iter_window(self, win_size=None, stride=None, **kwargs):
'''
Window iterator that yields data from the image based on win_size
and stride.
win_size and stride are both optional arguments. If neither are
passed, then the method pulls win_size from GDAL GetBlockSize() and
uses that to step through the image. If only win_size is provided,
the method yields adjoining windows of the size requested. If only
stride is provided, an error is rasied.
Parameters
----------
win_size : array-like, length 2, optional
The size of the requested image chip in x and y.
stride : array-like, length 2, optional
The size of the step between each yielded chip in x and y.
kwargs: optional
Arguments for get_data().
Yields
------
ndarray
Three dimensional numpy array of pixel values from the
requested region of the image.
'''
logger.debug('*** begin iter_window ***')
# Check input values
if win_size:
if any(x <= 0 for x in win_size):
raise ValueError('No value in win_size can be equal '
'to or less than zero.')
if stride:
if any(x <= 0 for x in stride):
raise ValueError('No value in stride can be equal '
'to or less than zero.')
# if NOT win_size and NOT stride
# use gdal to figure out block size and then continue on below.
if not win_size and not stride:
# Get block size from gdal
b = self._fobj.GetRasterBand(1)
win_size = b.GetBlockSize()
logger.debug('win_size is: %s, stride is: %s', win_size, stride)
# if win_size and NOT stride
# set stride to make windows adjoining
if win_size and not stride:
# Set vars for easy access below
xs = self.meta.shape[1]
ys = self.meta.shape[2]
xsize, ysize = win_size
# Find starting offsets by identifying the pixels that don't fit in
# the requested window blocks and then split the different
# between ends of the image using floor (int) to reduce fractions.
x_extra_pixels = xs % win_size[0]
xoff = int(x_extra_pixels / 2.0)
y_extra_pixels = ys % win_size[1]
yoff = int(y_extra_pixels / 2.0)
# Use while True to loop through get_data until outside the image
xoff_start = xoff
xsize, ysize = win_size
while True:
logger.debug(' xoff is %s,\tyoff is %s', xoff, yoff)
yield self.get_data(window=[xoff, yoff, xsize, ysize],**kwargs)
xoff += xsize
if xoff > self.meta.shape[1]:
xoff = xoff_start
yoff += ysize
if yoff > self.meta.shape[2]:
break
# if NOT win_size and stride, raise error
elif not win_size and stride:
raise ValueError('Setting stride and not setting win_size is not '
'allowed because there is no resonable value to '
'set win_size to. In this case stride can be '
'even or odd which could result in alternative '
'size return blocks around the center pixel '
'(or fractional pixel).')
# if win_size and stride
# just do it
elif win_size and stride:
# Set vars for easy access below
xs = self.meta.shape[1]
ys = self.meta.shape[2]
xsize, ysize = win_size
xstride, ystride = stride
# Find starting offset by identifying pixels that don't fit in
# the requested size/stride and then split the different between
# ends of the image using floor (int) to reduce fractions.
x_extra_pixels = (xs - xsize) % xstride
xoff = int(x_extra_pixels/2.0)
y_extra_pixels = (ys - ysize) % ystride
yoff = int(y_extra_pixels/2.0)
# Start the yield loop
xoff_start = xoff
while True:
logger.debug(' xoff is %s,\tyoff is %s', xoff, yoff)
yield self.get_data(window=[xoff, yoff, xsize, ysize], **kwargs)
xoff += xstride
if xoff > self.meta.shape[1]:
xoff = xoff_start
yoff += ystride
if yoff > self.meta.shape[2]:
break
def iter_window_random(self, win_size=None, no_chips=1000, **kwargs):
"""Random chip iterator.
Parameters
----------
win_size : array-like, length 2
The size of the requested image chip in x and y.
no_chips : int, optional
Number of chips to generate.
kwargs: optional
Arguments for get_data().
Yields
------
ndarray
Three dimensional numpy array of pixel values from the
requested region of the image.
"""
# Check input values
if win_size:
if any(x <= 0 for x in win_size):
raise ValueError('No value in win_size can be equal '
'to or less than zero.')
counter = no_chips
xs = self.meta.shape[1]
ys = self.meta.shape[2]
xsize, ysize = win_size
while True:
# select random offset
xoff = np.random.randint(xs-xsize+1)
yoff = np.random.randint(ys-ysize+1)
yield self.get_data(window=[xoff, yoff, xsize, ysize], **kwargs)
counter -= 1
if counter == 0: break
def iter_components(self, **kwargs):
"""This is a convenience method that iterataes (via yield) through
the components in the image object. Any kwargs valid for get_data
can be passed through.
kwargs can be any valid arugment for get_data
Parameters
----------
None
Yields
------
ndarray
Three dimensional numpy array of pixel values from the
requested region of the image.
"""
for c in xrange(len(self.files.dfile_tiles)):
yield self.get_data(component=c, **kwargs)
def iter_vector(self, vector=None, properties=False, filter=None, **kwargs):
"""This method iterates (via yeild) through a vector object or file.
Any kwargs valid for get_data can be passed through."""
if 'window' in kwargs.keys():
raise ValueError("The window argument is not valid for this " \
"method. They both define a retrieval " \
"geometry. Pass one or the other.")
if 'geom' in kwargs.keys():
raise ValueError("The geom argument is not valid for this " \
"method. The vector file passed in defines " \
"the retrieval geometry.")
# ToDo Test for overlap of geom and image data?
obj = ogr.Open(vector)
lyr = obj.GetLayer(0)
lyr_sr = lyr.GetSpatialRef()
img_proj = self.meta.projection_string
img_trans = self.meta.geo_transform
img_sr = osr.SpatialReference()
img_sr.ImportFromWkt(img_proj)
coord_trans = osr.CoordinateTransformation(lyr_sr, img_sr)
for feat in lyr:
# Return feature properties data is requested
if properties is True:
prop_out = feat.items()
elif properties:
if isinstance(properties, (list, tuple, str)):
if isinstance(properties, str):
properties = [properties]
it = feat.items()
if not all(x for x in properties if x in it.keys()):
raise ValueError("One or more of the requested "
"properties are not in the vector "
"feature.")
prop_out = {x: it[x] for x in properties if x in it.keys()}
if not prop_out:
prop_out = None
warnings.warn("No properties value found matching "
"request.")
else:
raise ValueError("Invalid properties argument.")
# Determine if the feature should be returned based on value of
# filter and if the value exists in the feature properties.
if filter:
# The filter should be either a list of dictionary key/value
# paris of length one, or of list of key/value pairs. The
# idea is that you can filter against more than one value
# of a key when you can pass a list of pairs.
if isinstance(filter, dict) & (len(filter) != 1):
raise ValueError("Filters should be passed in as a " \
"list of dictionaries that will " \
"be used to filter against the " \
"feature property values.")
# If filter is a dict of len 1, convert to a list of len 1 for
# the looping code below.
if isinstance(filter,dict):
filter = [filter]
# Get feature properties to check against.
prop_test = feat.items()
# raise warning if a filter item key is not in properties
if any(f.keys()[0] not in prop_test.keys() for f in filter):
warnings.warn("Requested filter key is not present in "
"vector properties.")
# If any filter is caught, pass on, otherwise return
if any(prop_test.get(d.keys()[0], None) ==
d.values()[0] for d in filter):
pass
else:
if properties:
yield (None, None)
continue
else:
yield None
continue
# Get the geometry to pass to get_data
geom = feat.geometry()
# Use transform from above to put geom in image space
geom.Transform(coord_trans)
# Catch and pass OverlapError for the iterator
try:
data = self.get_data(geom=geom, **kwargs)
except OverlapError:
data = None
# Yield the data
if properties:
yield (data, prop_out)
else:
yield data
def get_data_from_coords(self, coords, **kwargs):
"""Method to get pixel data given polygon coordintes in the same projection as
the image. kwargs can be anything accepted by get_data.
Parameters
----------
coords : list of lists. polygon coordinates formatted as follows:
- lat/long (EPSG:4326) projection: [[lon_1, lat_1], [lon_2, lat_2], ...]
- UTM projection: [[x_1, y_1], [x_2, y_2], ...]
Returns
------
ndarray
Three dimensional numpy array of data from region of the image specified
in the feature geometry.
"""
# create ogr geometry ring from feature geom
ring = ogr.Geometry(ogr.wkbLinearRing)
for point in coords:
lon, lat = point[0], point[1]
ring.AddPoint(lon, lat)
ring.AddPoint(coords[0][0], coords[0][1]) # close geom ring with first point
geom = ogr.Geometry(ogr.wkbPolygon)
geom.AddGeometry(ring)
return self.get_data(geom=geom, **kwargs)
def get_data_from_vec_extent(self, vector=None, **kwargs):
"""This is a convenience method to find the extent of a vector and
return the data from that extent. kwargs can be anything accepted
by get_data."""
if vector is None:
raise ValueError("Requires a vector to read. The vector can be " \
"a string that describes a vector object or a " \
"path to a valid vector file.")
if 'window' in kwargs.keys():
raise ValueError("The window argument is not valid for this " \
"method. The vector file passed in defines " \
"the retrieval geometry.")
if 'geom' in kwargs.keys():
raise ValueError("The geom argument is not valid for this " \
"method. The vector file passed in defines " \
"the retrieval geometry.")
if 'mask' in kwargs.keys():
raise ValueError("A mask request is not valid for this method " \
"because it retrives data from the full extent " \
"of the vector. You might want a rasterize " \
"method or iter_vector?")
# ToDo Test for overlap of geom and image data?
obj = ogr.Open(vector)
lyr = obj.GetLayer(0)
lyr_sr = lyr.GetSpatialRef()
img_proj = self.meta.projection_string
img_trans = self.meta.geo_transform
img_sr = osr.SpatialReference()
img_sr.ImportFromWkt(img_proj)
coord_trans = osr.CoordinateTransformation(lyr_sr,img_sr)
extent = self.ogr_extent_to_extent(lyr.GetExtent())
window = self.extent_to_window(extent,coord_trans)
[xoff, yoff, win_xsize, win_ysize] = window
return self.get_data(window = window, **kwargs)
def ogr_extent_to_extent(self,ogr_extent):
return [ogr_extent[0],ogr_extent[2],ogr_extent[1],ogr_extent[3]]
def extent_to_window(self,extent,coord_trans=None):
"""
Convert an extent in the form [ul_x, ul_y, lr_x, lr_y] in projection
space to a window of the form [xoff, yoff, win_xsize, win_ysize].
If a coord_trans is provided, it is used to convert between projections
before the window calculation is made. This is provided since the
creation of this object can be somewhat expensive.
Parameters
----------
extent : list
extent in the form [ul_x, ul_y, lr_x, lr_y] to be converted to
a window list.
coord_trans : osr object
osr object that converts between two image projections.
Returns
-------
list
list that defines a pixel window of the form
[xoff, yoff, win_xsize, win_ysize].
"""
# Set corner vectors
ul_vec = extent[:2]
lr_vec = extent[2:]
if coord_trans:
# Online dobumentation says that there could be a bug in
# TransformPoints - seems to be working with 1.11.4 so
# I'll leave it in for now. When working with a bag vector file
# TransformPoints returned inf values that exploded below whereas
# the sequental TransformPoint returned values that triggered
# the overlap error below. Leaving the single call for now
# as it is faster and cleaner.
[ul_img, lr_img] = coord_trans.TransformPoints([ul_vec, lr_vec])
#ul_img = coord_trans.TransformPoint(*ul_vec)
#lr_img = coord_trans.TransformPoint(*lr_vec)
# Cut off third returned dimension (should be zero)
ul_img = ul_img[:-1]
lr_img = lr_img[:-1]
else:
ul_img = ul_vec
lr_img = lr_vec
xs,ys = self.proj_to_raster(*zip(*[ul_img,lr_img]))
xoff = int(math.floor(min(xs)))
yoff = int(math.floor(min(ys)))
xmax = int(math.ceil(max(xs)))
ymax = int(math.ceil(max(ys)))
win_xsize = xmax-xoff
win_ysize = ymax-yoff
window = [xoff, yoff, win_xsize, win_ysize]
# Logging info if needed
logger.debug('vector geo extent...\n\t%s\n\t%s',ul_vec,lr_vec)
logger.debug('image geo extent...\n\t%s\n\t%s',ul_img,lr_img)
logger.debug('geo transform...\n\t%s', self.meta.geo_transform)
logger.debug('raster xy extent...\n\t%s,\n\t%s',xs,ys)
logger.debug('requested window...\n\t%s',window)
if ((xoff + win_xsize <= 0) or (yoff + win_ysize <= 0) or
(xoff > self.meta.shape[1]) or (yoff > self.meta.shape[2])):
raise OverlapError("The requested data window has no " \
"content. Perhaps the image and vector " \
"do not overlap or the projections may " \
"not be able to be automatically reconciled?")
return window
def proj_to_raster(self, projx, projy):
'''
Method to convert points in projection space to points in raster space.
Input can be in a variety of types as long as both the input parameters
are of the same type. The method will attempt to return a data type
as similar as possible to the input type.
Parameters
----------
projx : int, float, list, tuple, or numpy.ndarray
Input point in projected space.
projy : int, float, list, tuple, or numpy.ndarray
Input point in projected space.
Returns
-------
x : float, list, tuple, or numpy.ndarray
raster x value calculated from projx, projy, and the object's
geo_transform
y : float, list, tuple, or numpy.ndarray
raster y value calculated from projx, projy, and the object's
geo_transform
'''
# This method can't handle mixed types
if type(projx) != type(projy):
raise ValueError('The type of x and y should be the same and '
'either integers, float, tuples, lists, or '
'numpy arrays.')
# Determine input type so that the conversion back can be done
# on return
intype = 0
if not hasattr(projx, '__iter__'):
intype = 1
elif isinstance(projx, list):
intype = 2
elif isinstance(projx, tuple):
intype = 3
elif isinstance(projx, np.ndarray):
intype = 4
else:
raise ValueError("The input type was not recognized.")
# Convert to numpy arrays for the calculation
projx = np.asarray(projx)
projy = np.asarray(projy)
# Get geo_transform from object
gm = self.meta.geo_transform
# Transform per inverse of http://www.gdal.org/gdal_datamodel.html
x = (gm[5] * (projx - gm[0]) - gm[2] * (projy - gm[3])) / \
(gm[5] * gm[1] + gm[4] * gm[2])
y = (projy - gm[3] - x * gm[4]) / gm[5]
# Return to input type
if intype == 1:
return float(x), float(y)
elif intype == 2:
return list(x), list(y)
elif intype == 3:
return tuple(x), tuple(y)
elif intype == 4:
return x, y
def raster_to_proj(self, x, y):
'''
Method to convert points in raster space to points in projection space.
Input can be in a variety of types as long as both the input parameters
are of the same type. The method will attempt to return a data type
as similar as possible to the input type.
Parameters
----------
x : int, float, list, tuple, or numpy.ndarray
Input point in raster space.
y : int, float, list, tuple, or numpy.ndarray
Input point in raster space.
Returns
-------
projx : float, list, tuple, or numpy.ndarray
projection x value calculated from x, y, and the object's
geo_transform
projy : float, list, tuple, or numpy.ndarray
projection y value calculated from x, y, and the object's
geo_transform
'''
# This method can't handle mixed types
if type(x) != type(y):
raise ValueError('The type of x and y should be the same and '
'either integers, float, tuples, lists, or '
'numpy arrays.')
# Determine input type so that the conversion back can be done
# on return
intype = 0
if not hasattr(x, '__iter__'):
intype = 1
elif isinstance(x, list):
intype = 2
elif isinstance(x, tuple):
intype = 3
elif isinstance(x, np.ndarray):
intype = 4
else:
raise ValueError("The input type was not recognized.")
# Convert to numpy arrays for the calculation
x = np.asarray(x)
y = np.asarray(y)
# Get geo_transform from object
gm = self.meta.geo_transform
# Transform per http://www.gdal.org/gdal_datamodel.html
projx = gm[0] + gm[1] * x + gm[2] * y
projy = gm[3] + gm[4] * x + gm[5] * y
# Return to input type
if intype == 1:
return float(projx), float(projy)
elif intype == 2:
return list(projx), list(projy)
elif intype == 3:
return tuple(projx), tuple(projy)
elif intype == 4:
return projx, projy
def get_data(self, component = None,
bands = None,
window = None,
buffer = None,
geom = None,
mask = False,
mask_all_touched = False,
return_location = False,
boundless = True,
virtual = False):
"""
The central method for extracting pixel data from an image. It
provides several options that define the region of the image from
which pixels are pulled and how the data is returned to the users.
The returned numpy array follows the gdal convention using a
bands-first format (bands, x, y). If a single band is requested,
the first dimension will be singular. Numpy squeeze can be used
to quickly remove the singular dimension. The bands-first format
can be quickly converted to a bands-last format to more easily work
with other image processing packages using
tinytools.np_img.conv_to_bandslast (or the short block of numpy code
the function calls).
There isn't generally a reason to call both component and window,
but it isn't explicitly disallowed. If window and component are both
specified, the resulting data window is relative to the coordinates
of the specified component.
Spectral processing (i.e. at sensor radiance and TOA reflectance) and
band aliasing are availble in the higher level classes that
understand the meta data of specific sensors (i.e. DGImage).
Parameters
----------
component : int
Retrive data from a specific tile (base 1). If the file in not
a tiled format, there is only a single component. Order of the
tiles is stored in GeoImage.files.dfile_tiles.
bands : list of ints
The band numbers (base 1) to be retrieved.
window : list of ints
The window from which to retrive data in the format
[xoff, yoff, x_window_size, y_window_size]. An error will be
raise in no valid pixels are requested. If the window is
partially valid, the requested window will be padded.
buffer : int or length two list of ints
Number of pixels to add as a buffer around the requested pixels.
If an int, the same number will be added to both dimensions. If
a list of two ints, they will be interprested as xpad and ypad.
geom : valid ogr or fiona geometry
A shape geometry to define the pixel retrieval region. Must be
in the same projection as the image as no projection checking
can be done on a geometry object.
mask : bool
Should the returned array be a masked numpy array? If requested,
invalid pixels will be masked as defined by pixels either
outside the image or outside the geometry, if passed.
mask_all_touched : bool
Should the mask operation use an all-touched algorithm. Default
is center-touched.
virtual : bool
Should the returned numpy array use the gdal virtual raster
driver? This is currently only an option on linux.
return_location : bool
Should the returned numpy array be returned in a tuple along
with the location information of the requested image region.
boundless : bool
Should pixels outside the image be returned as masked/padded
values (True, default) or should an error be raised (False).
Returns
------
ndarray
Three dimensional numpy array of data from the requested region
of the image.
tuple
If return_location is True, a tuple will be returned with the
image numpy array as well as a dictionary specifying the location
information.
"""
if component is not None:
if component == 0:
raise ValueError("Component should be specified as based 1.")
if component > len(self.files.dfile_tiles):
raise ValueError("You've requested a component value greater "
"than the number available.")
y = GeoImage(self.files.dfile_tiles[component-1])
logger.debug('returning data from: '+
str(self.files.dfile_tiles[component-1]))
obj = y._fobj
else:
obj = self._fobj
# If bands not requested, set to pull all bands
if not bands:
bands = range(1,obj.RasterCount+1)
# Set window to pull
if window and geom:
raise ValueError("The arguments window and geom are mutually " \
"exclusive. They both define an image " \
"extent. Pass either one or the other.")
if window:
# Set extent parameters based on window if provided
if len(window) == 4:
[xoff, yoff, win_xsize, win_ysize] = window
else:
raise ValueError("Window must be length four and will be read" \
"as: xoff, yoff, win_xsize, win_ysize")
elif geom:
# Set window size based on a geom object in image space
# ToDo - Add all_touched option to this and mask function.
g = self._instantiate_geom(geom)
extent = self.ogr_extent_to_extent(g.GetEnvelope())
window = self.extent_to_window(extent)
[xoff, yoff, win_xsize, win_ysize] = window
else:
# Else use extent of image to set extent params
xoff = 0
yoff = 0
win_xsize = self.meta.shape[1]
win_ysize = self.meta.shape[2]
# Add buffer
if buffer:
if isinstance(buffer,int):
buffer = [buffer]
if len(buffer) == 1:
xbuff = buffer[0]
ybuff = buffer[0]
elif len(buffer) == 2:
xbuff = buffer[0]
ybuff = buffer[1]
else:
raise ValueError("Buffer must be either length one or two.")
# Apply the buffer to the readasarray parameters
xoff = xoff-xbuff
yoff = yoff-ybuff
win_xsize = win_xsize+2*xbuff
win_ysize = win_ysize+2*ybuff
# Handle out-of-bounds cases
# (i.e. xoff = 0; buffer = 3; xoff - buffer)
if boundless:
# initialize buffer vars
np_xoff_buff = 0
np_yoff_buff = 0
np_xlim_buff = 0
np_ylim_buff = 0
if xoff < 0:
np_xoff_buff = xoff
win_xsize = win_xsize + xoff
xoff = 0
if yoff < 0 :
np_yoff_buff = yoff
win_ysize = win_ysize + yoff
yoff = 0
xpos = xoff+win_xsize
xlim = self.meta.shape[1]
ypos = yoff+win_ysize
ylim = self.meta.shape[2]
if xpos > xlim:
np_xlim_buff = xpos-xlim
win_xsize = win_xsize-np_xlim_buff
if ypos > self.meta.shape[2]:
np_ylim_buff = ypos-ylim
win_ysize = win_ysize-np_ylim_buff
else:
# else set the buffer vars to zero
np_xoff_buff = 0
np_yoff_buff = 0
np_xlim_buff = 0
np_ylim_buff = 0
# Check platform if virtual is requested
if virtual:
if platform.system() != 'Linux':
raise ValueError('The gdal virtual driver can only run '
'on the linux platform.')
###############################
##### Read the image data #####
###############################
if not virtual:
# Read data one band at a time with ReadAsArray
zt = len(bands)
dt = const.DICT_GDAL_TO_NP[self.meta.gdal_dtype]
data = np.empty([zt, win_ysize, win_xsize], dtype=dt)
for i, b in enumerate(bands):
bobj = obj.GetRasterBand(b)
data[i, :, :] = bobj.ReadAsArray(xoff=xoff,
yoff=yoff,
win_xsize=win_xsize,
win_ysize=win_ysize)
elif virtual:
data = obj.GetVirtualMemArray(gdalconst.GF_Write,
xoff=xoff,
yoff=yoff,
xsize=win_xsize,
ysize=win_ysize,
bufxsize=win_xsize,
bufysize=win_ysize,
band_list=bands)
# Since the virtual api, returns 2d if only one band is requested,
# expand to the geoio standard 3d image.
if data.ndim == 2:
data = data[np.newaxis, :, :]
# Convert numpy array to masked numpy array if requested.
if mask and geom:
# Set image parameters
xres = self.meta.resolution[0]
yres = self.meta.resolution[1]
(xmin, xmax, ymin, ymax) = g.GetEnvelope()
ul_env = [xmin, ymax]
ul_raster = self.proj_to_raster(*ul_env)
ul_corner = [math.floor(ul_raster[0]),math.floor(ul_raster[1])]
xmin_corner,ymax_corner = self.raster_to_proj(*ul_corner)
# Create temporary raster to burn
drv = gdal.GetDriverByName('MEM')
tds = drv.Create('', win_xsize, win_ysize, 1, gdal.GDT_Byte)
tds.SetGeoTransform((xmin_corner, xres, 0, ymax_corner, 0, -yres))
tds.SetProjection(self.meta.projection_string)
# Create ogr layr from geom
odrv = ogr.GetDriverByName('Memory')
ds = odrv.CreateDataSource('')
ltype = g.GetGeometryType()
lsrs = osr.SpatialReference(self.meta.projection_string)
lyr = ds.CreateLayer('burnshp', lsrs, ltype)
feat = ogr.Feature(lyr.GetLayerDefn())
feat.SetGeometryDirectly(g)
lyr.CreateFeature(feat)
# Run the burn
if mask_all_touched:
err = gdal.RasterizeLayer(tds, [1], lyr, burn_values=[1],
options=['ALL_TOUCHED=TRUE'])
else:
err = gdal.RasterizeLayer(tds, [1], lyr, burn_values=[1])
# build the masked array
m = tds.ReadAsArray().astype('bool')
data = np.ma.array(data, mask=np.tile(~m, (data.shape[0], 1, 1)))
if mask and not geom:
# This code will mask values outside the image as well as zeros
# inside the image.
data = np.ma.array(data,mask=~data.astype('bool'))
# The line below will only mask values outside the image.
# data = np.ma.array(data,mask=np.zeros(data.shape).astype('bool'))
# Pad the output array if needed
pad_tuples = ((0,0),
(np.abs(np_yoff_buff), np.abs(np_ylim_buff)),
(np.abs(np_xoff_buff), np.abs(np_xlim_buff)))
if not mask and np.any([np.any(x) for x in pad_tuples]):
data = np.pad(data, pad_tuples, 'constant',constant_values=0)
elif mask:
mpad = np.pad(data.mask, pad_tuples, 'constant', constant_values=1)
data = np.pad(data, pad_tuples, 'constant', constant_values=0)
data = np.ma.array(data,mask=mpad)
# if "y" is open, close it
try:
y = None
except NameError:
pass
if return_location:
# Find upper-left pixel location including buffers is necessary.
ul_pix_x = xoff+np_xoff_buff
ul_pix_y = yoff+np_yoff_buff
# Create new geotransform tuple.
gt = self.meta.geo_transform
(lx,ly) = self.raster_to_proj(ul_pix_x,ul_pix_y)
gt_new = (lx, gt[1], gt[2], ly, gt[4], gt[5])
# Set the dictionary to be returned below
location_dict = {}
location_dict['upper_left_pixel'] = [ul_pix_x,ul_pix_y]
location_dict['geo_transform'] = gt_new
return data, location_dict
else:
return data
def _instantiate_geom(self,g):
"""Attempt to convert the geometry pass in to and ogr Geometry
object. Currently implements the base ogr.CreateGeometryFrom*
methods and will reform fiona geometry dictionaries into a format
that ogr.CreateGeometryFromJson will correctly handle.
"""
if isinstance(g,ogr.Geometry):
# If the input geometry is already an ogr object, create a copy
# of it. This is requred because of a subtle issue that causes
# gdal to crash if the input geom is used elsewhere. The working
# theory is that the geometry is freed when going out of a scope
# while it is needed in the upper level loop. In this code, the
# problem comes about between self.iter_vector and self.get_data
# with mask=True.
return ogr.CreateGeometryFromJson(str(g.ExportToJson()))
# Handle straight ogr GML
try:
return ogr.CreateGeometryFromGML(g)
except:
pass
# Handle straight ogr Wkb
try:
return ogr.CreateGeometryFromWkb(g)
except:
pass
# Handle straight ogr Json
try:
return ogr.CreateGeometryFromJson(g)
except:
pass
# Handle straight ogr Wkt
try:
return ogr.CreateGeometryFromWkt(g)
except:
pass
# Handle fiona Json geometry format
try:
gjson = str(g).replace('(','[').replace(')',']')
return ogr.CreateGeometryFromJson(gjson)
except:
pass
raise ValueError("A geometry object was not able to be created from " \
"the value passed in.")
def downsample(self,arr=None,
shape = None,
factor = None,
extent = None,
method = 'aggregate',
no_data_value = None,
source = None):
"""
Downsample array with function from geoio.downsample.py. The methods
available are 'aggregate' and 'nearest'. If the data arr is not
specified, all bands will be retrieved from the image object. For
detailed documentation, see geoio.downsample.downsample.
Parameters
----------
arr : array_like
Image data in a two or three dimension array in band-first order.
If it is not provided, the data will be pulled via get_data from
the current image object.
shape : list
Shape of the desired output array.
factor : integer, float, or length two iterable
Factor by which to scale the image (must be less than one).
extent : length four array-like
List of upper-left and lower-right corner coordinates for the edges
of the resample. Coordinates should be expressed in pixel space.
i.e. [-1,-2,502,501]
method : strings
Method to use for the downsample - 'aggregate', 'nearest', 'max', or
'min.
no_data_value : int
Data value to treat as no_data
source : strings
Package to use for algorithm - opencv ('cv2') or custom 'numba' below.
Returns
-------
ndarray
Three dimensional numpy array of downsampled data
"""
# importing downsample trigger numba compile - which slows down the
# whole package import if it is done at import. So, only pull in
# and compile when the code is necessary.
import downsample
if arr is None:
logger.debug('Array not passed in, retrieving all bands...')
arr = self.get_data()
return downsample.downsample(arr,
shape = shape,
factor = factor,
extent = extent,
method = method,
no_data_value = no_data_value,
source = source)
def downsample_like_that(self,ext_img,arr = None,**kwargs):
"""Pull downsampling parameters from another geoio image object and
use them to run downsampling. Currently assumes that they are in the
same projection. kwargs can be the downsample method requested,
no_data_value, or source arguments that are valide for
geoio.downsample.downsample.
"""
# importing downsample trigger numba compile - which slows down the
# whole package import if it is done at import. So, only pull in
# and compile when the code is necessary.
import downsample
####
# Eventually add code to check and unify projection
####
gt = self.meta.geo_transform
# Set info for corners call
ext_gt = ext_img.meta.geo_transform
ext_res = ext_img.meta.resolution
ext_shape = ext_img.meta.shape[1:]
ul_corner = ext_img.meta.extent[:2]
#ul_corner = [ext_img.meta.extent[x] for x in [0,2]]
ul_corner_pix = self.proj_to_raster(ul_corner[0],ul_corner[1])
lr_corner = ext_img.meta.extent[2:]
#lr_corner = [ext_img.meta.extent[x] for x in [1,3]]
lr_corner_pix = self.proj_to_raster(lr_corner[0],lr_corner[1])
ext_corners = [ul_corner_pix, lr_corner_pix, ext_res]
if arr is None:
arr = self.get_data()
if kwargs.has_key('no_data_value'):
no_data_value = kwargs['no_data_value']
del kwargs['no_data_value']
elif self.meta.no_data_value:
no_data_value = self.meta.no_data_value
else:
no_data_value = None
return downsample.downsample(arr,extent=ul_corner_pix+lr_corner_pix,
shape=ext_shape,
**kwargs)
def downsample_to_grid(self,arr,x_steps,y_steps,no_data_value=None,
method='aggregate',source=None):
"""Direct call to grid downsample for geoio.downsample.
Parameters
----------
arr : array_like or None
Image data in a two or three dimension array in band-first order.
If None is provided, the data will be pulled via get_data from
the current image object.
x_steps : array_like
The x-dim resample edges in image pixel space, can be outside the
dimension of the current image. i.e. [-1.5,0,1.5,3.0]
y_steps : array_like
The y-dim resampled edges inimage pixel space as x_steps.
no_data_value : int or float
The no_data_value to be handled in the downsample method
method : string
Requested downsample method - can be 'aggregate', 'nearest', 'max',
or 'min'. Not all options are available will all x_steps and
y_steps requests.
source : string
What source to pull the downsample algorithm from. If left blank,
the is automatically chosen based on the env and request
parameters. The options are 'numba' or 'cv2'.
Returns
-------
ndarray
Three dimensional numpy array of downsampled data
"""
# importing downsample trigger numba compile - which slows down the
# whole package import if it is done at import. So, only pull in
# and compile when the code is necessary.
import downsample
if arr is None:
logger.debug('Array not passed in, retrieving all bands...')
arr = self.get_data()
return downsample.downsample_to_grid(arr,
x_steps,
y_steps,
no_data_value=no_data_value,
method=method,
source=source)
def upsample(self, arr=None,
shape=None,
factor=None,
extent=None,
method='bilinear',
no_data_value=None):
"""Method to call gdal upsample operations. The output of this
should be at or above the resolution of the input data. Otherwise,
the downsampling code should be used. Alternatively, the resample
method can be used to automatically choose reasonable defaults.
Parameters
----------
arr : array_like
Image data in a two or three dimension array in band-first order.
If it is not provided, the data will be pulled via get_data from
the current image object.
shape : list
Shape of the desired output array.
factor : integer, float, or length two iterable
Factor by which to scale the image (must be greater than one).
extent : length four array-like
List of upper-left and lower-right corner coordinates for the edges
of the resample. Coordinates should be expressed in pixel space.
i.e. [-1,-2,502,501]
method : strings
Method to use for the upsample - 'nearest', 'bilinear', 'cubic',
or 'average'.
no_data_value : int
Data value to treat as no_data
Returns
-------
ndarray
Three dimensional numpy array of downsampled data
"""
if factor is not None and factor<1:
raise ValueError('factor should be equal to or greater than one '
'for resampling/upsampling operations.')
if factor:
if isinstance(factor,(int,float)):
factor = [factor,factor]
shape = [int(math.ceil(self.shape[1]*factor[0])),
int(math.ceil(self.shape[1]*factor[1]))]
elif shape:
pass
if extent is not None:
ul_projx, ul_projy = self.raster_to_proj(extent[0], extent[1])
lr_projx, lr_projy = self.raster_to_proj(extent[2], extent[3])
else:
ul_projx, ul_projy = self.meta.extent[:2]
lr_projx, lr_projy = self.meta.extent[2:]
# extent is:
# ul_x, ul_y, lr_x, lr_y
# geo_transform is
# ul_projx, xres, xangle, ul_proj, yangle, yres
# (487546.8, 1.2, 0.0, 4443553.199999999, 0.0, -1.2)
gt = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
gt[0] = ul_projx
gt[3] = ul_projy
gt[1] = (lr_projx - ul_projx) / shape[1]
gt[5] = (lr_projy - ul_projy) / shape[0]
#import ipdb; ipdb.set_trace()
# Set up destination image in memory
drv = gdal.GetDriverByName('MEM')
dst = drv.Create('',
shape[1],
shape[0],
self.shape[0],
self.meta.gdal_dtype)
dst.SetGeoTransform(gt)
dst.SetProjection(self.meta.projection_string)
if no_data_value is not None:
blist = [x + 1 for x in range(ext_img.shape[0])]
[dst.GetRasterBand(x).SetNoDataValue(no_data_value) for x in blist]
# Run resample, pull data, and del the temp gdal object.
data = self._upsample_from_gdalobj(self.get_gdal_obj(), dst, method)
del dst
return data
def upsample_like_that(self, ext_img, method='bilinear', no_data_value=None):
"""Use gdal.ReprojectImage to upsample the object so that it looks
like the geoio image object passed in as the ext_img argument. Method
can be those listed in GeoImage.upsample.
"""
if ext_img.meta.resolution[0]>self.meta.resolution[0]:
raise ValueError('The requested resolution is not at or higher '
'than the base object. Use downsample or '
'resample methods.')
# Set up destination image in memory
drv = gdal.GetDriverByName('MEM')
dst = drv.Create('',
ext_img.shape[2],
ext_img.shape[1],
ext_img.shape[0],
ext_img.meta.gdal_dtype)
dst.SetGeoTransform(ext_img.meta.geo_transform)
dst.SetProjection(ext_img.meta.projection_string)
if no_data_value is not None:
blist = [x + 1 for x in range(ext_img.shape[0])]
[dst.GetRasterBand(x).SetNoDataValue(no_data_value) for x in blist]
# Run resample, pull data, and del the temp gdal object.
data = self._upsample_from_gdalobj(self.get_gdal_obj(),dst,method)
del dst
return data
def _upsample_from_gdalobj(self,src,dst,method='bilinear'):
"""Hidden to run the actual reprojection gdal code that is called
from two higher level methods."""
# Set reprojection method
if isinstance(method,int):
pass
elif method == "nearest":
method = gdal.GRA_NearestNeighbour
elif method == "bilinear":
method = gdal.GRA_Bilinear
elif method == "cubic":
method = gdal.GRA_Cubic
elif method == "average":
method = gdal.GRA_Average
else:
raise ValueError("requested method is not understood.")
# Do the reprojection
gdal.ReprojectImage(src,
dst,
self.meta.projection_string,
dst.GetProjection(),
method)
# Return data and free the temp image.
return dst.ReadAsArray()
def resample_like_that(self,ext_img):
"""Method to choose resonable default for up or downsample operations.
"""
if (ext_img.meta.resolution[0]<self.meta.resolution[0]) and \
(ext_img.meta.resolution[1]<self.meta.resolution[1]):
return self.upsample_like_that(ext_img)
elif (ext_img.meta.resolution[0]>self.meta.resolution[0]) and \
(ext_img.meta.resolution[1]>self.meta.resolution[1]):
return self.downsample_like_that(ext_img)
else:
raise ValueError('Mixed resampling between the x and y '
'dimensions is not supported.')
def write_img_like_this(self,new_fname,np_array,return_obj=False,
gdal_driver_name=None,options=[],
vrt_fallback="GTiff"):
"""Write the data passed in the variable "np_array" as a new image
with image parameters (projection, etc.) pulled from this object.
The new file is created as the data type of the numpy array
"np_array". "np_array" should have the same line/sample size as
the object data, but the number of bands can vary. This method uses
the create_geo_image function in this module. That function uses
gdal create and has minimal metadata copy since what should or
should not be copied is application dependent.
"""
# If the data array is 2D, add a dimension for a single band
if np_array.ndim == 2:
np_array = np_array[np.newaxis, :, :]
elif np_array.ndim == 3:
pass
else:
raise ValueError("This can't handle Arrays larger than three "
"dimensions.")
if not ((np_array.shape[1] == self.shape[2]) and \
(np_array.shape[2] == self.shape[1])):
raise ValueError("Data in is not the same size as that in the "
"current image object. Data in is shape: %s, "
"Object shape is: %s" %
(np_array.shape, self.shape))
if gdal_driver_name is None:
gdal_driver_name = self.meta.driver_name
## Write the geo image using my geoio function
# NDV set to 0 by default in create_geo_image
new_img_name = create_geo_image(new_file_name = new_fname,
data_np_array = np_array,
gdal_driver_name = gdal_driver_name,
#gdal_driver_name = 'GTiff',
gdal_geo_t = self.meta.geo_transform,
gdal_projection = self.meta.projection_string,
data_type = np_array.dtype,
NDV=self.meta.no_data_value,
options=options,
vrt_fallback=vrt_fallback)
# Output information about the new file
f = GeoImage(new_img_name)
logger.debug('')
logger.debug("######################")
logger.debug("New image file has been created at: "+new_img_name)
logger.debug("Data type is: "+str(np_array.dtype))
if self.meta.no_data_value:
logger.debug("No data value set to: " +
str(self.meta.no_data_value))
else:
logger.debug("No data value was not set.")
logger.debug("######################")
if return_obj:
return f
def write_img_replace_this(self,np_array):
"""Replace the data in the current object image with the data passed
in the variable "np_array". This method uses gdal to replace the
data in place - so all metadata is intact.
"""
# All dims must be exact
if not ((np_array.shape[0] == self.shape[0]) and \
(np_array.shape[1] == self.shape[2]) and \
(np_array.shape[2] == self.shape[1])):
raise ValueError("Data in is not the same shape as that in the "
"current image object. Data in is shape: %s, "
"Object shape is: %s" %
(np_array.shape, self.shape))
# Check that the incoming array is the same dtype as the image.
gi_dtype = self.meta.gdal_dtype
if not (np_array.dtype == const.DICT_GDAL_TO_NP[gi_dtype]):
raise ValueError("Data types must match exactly to do a "
"data replace.")
# Reload gdal object in update mode
reload_fname = self.meta_fname
self._fobj = None
self._fobj = gdal.Open(reload_fname, gdalconst.GA_Update)
# Load new data into gdal object file
nbands = self.shape[0]
for i in xrange(nbands):
b = self._fobj.GetRasterBand(i+1).WriteArray(np_array[i,:,:])
b = None
# Dereference gdal object and reopen as read only
self._fobj = None
self._fobj = gdal.Open(reload_fname, gdalconst.GA_ReadOnly)
# ToDo Add ability to repoject data during file write.
def get_stretch_values(self,**kwargs):
# Call the stretch_values function in this module.
return get_img_stretch_vals(self._fobj,**kwargs)
def get_gdal_obj(self):
"""Return gdal object hidden under geoio object."""
return self._fobj
def read_geo_file_info(fname_or_fobj):
""" Get image metadata."""
# class : RasterBrick
# dimensions : 3191, 921, 2938911, 11 (nrow, ncol, ncell, nlayers)
# resolution : 30, 30 (x, y)
# extent : 630885, 658515, 2796285, 2892015 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=40 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
# data source : /path/to/file/desertEO1H1580422003134110PZ.TIFF
# names : desertEO1H1580422003134110PZ.1, desertEO1H1580422003134110PZ.2, desertEO1H1580422003134110PZ.3, desertEO1H1580422003134110PZ.4, desertEO1H1580422003134110PZ.5, desertEO1H1580422003134110PZ.6, desertEO1H1580422003134110PZ.7, desertEO1H1580422003134110PZ.8, desertEO1H1580422003134110PZ.9, desertEO1H1580422003134110PZ.10, desertEO1H1580422003134110PZ.11
# If the filename passed in is a string, try to open the file, else you
# can assume that it is a file object already opened and passed in.
if isinstance(fname_or_fobj,str):
fobj = gdal.Open(fname_or_fobj, gdalconst.GA_ReadOnly)
else:
fobj = fname_or_fobj
summary = {}
summary['file_name'] = fobj.GetDescription()
summary['file_list'] = fobj.GetFileList()
summary['driver_name'] = fobj.GetDriver().ShortName
summary['no_data_value'] = fobj.GetRasterBand(1).GetNoDataValue()
summary['gdal_dtype'] = fobj.GetRasterBand(1).DataType
summary['gdal_dtype_name'] = gdal.GetDataTypeName(summary['gdal_dtype'])
### Get Image basics dimensions
x = fobj.RasterXSize
y = fobj.RasterYSize
summary['pixels'] = x*y
nbands = fobj.RasterCount
summary['shape'] = (nbands, x, y)
### Get image resolution and extents
# Pulled from GDALInfoReportCorner() at:
# http://www.gdal.org/gdalinfo_8c.html
gt = fobj.GetGeoTransform()
summary['geo_transform'] = gt
xstart = gt[0]
xend = gt[0] + gt[1]*x + gt[2]*y
ystart = gt[3]
yend = gt[3] + gt[4]*x + gt[5]*y
# Pixel Resolutions
xres = abs(gt[1])
yres = abs(gt[5])
# Add to summary
summary['resolution'] = (xres, yres)
summary['extent'] = (xstart, ystart, xend, yend)
### Get image projection and datum
try:
summary['projection_string'] = fobj.GetProjection()
# Get pretty printable projection_string
sr = osr.SpatialReference(fobj.GetProjectionRef())
summary['pprint_proj_string'] = sr.ExportToPrettyWkt()
summary['authority'] = sr.GetAttrValue("AUTHORITY",0)+ ':' + \
sr.GetAttrValue("AUTHORITY",1)
except:
# Used to handle 1b files without projections
pass
return summary
# Function to write a new file.
def create_geo_image(new_file_name, data_np_array, gdal_driver_name,
gdal_geo_t, gdal_projection, data_type, NDV=0,
options=[],vrt_fallback="GTiff"):
""" Create a new image file with all the parameters pass in. This
function is implemented in the GeoImage Class. It will autopopulate the
parameters so is a bit easier to use from there.
"""
if gdal_driver_name == "VRT":
# Change driver name to something that will write to disk. This
# defaults to GTiff but is configurable.
logger.debug("You have requested a VRT. Creating a single new real "
"file using %s format. This format can be ovridden "
"in the call using 'GTiff', 'ENVI', etc. Eventually, "
"this could be replaced with a file-for-file VRT "
"rewrite, but that does not always make sense in VRTs "
"since there can be pixels in the original image "
"(overlaps) that aren't in the "
"new data array.",vrt_fallback)
gdal_driver_name = vrt_fallback
# Strip ".vrt" extension if it exists and rename according to the
# driver that has been substituted.
# --- Only supports ENVI and GTiff now, I can't figure out how to do
# this dynamically in gdal. (as of 141218)
if tt.files.filter(new_file_name, '*.VRT', case_sensitive=False):
if gdal_driver_name == "GTiff":
new_file_name = os.path.splitext(new_file_name)[0] + ".TIF"
elif gdal_driver_name == "ENVI":
new_file_name = os.path.splitext(new_file_name)[0]
else:
new_file_name = os.path.splitext(new_file_name)[0]
## Convert data types to the appropriate value. The gdal create
## routine can pass either the gdal data type int or the alias (such as
## gdal.GDT_Float32).
# If the data type passed is of type numpy, convert it to GdalDataType
if data_type in const.DICT_NP_TO_GDAL:
data_type = const.DICT_NP_TO_GDAL[data_type]
# If gdal dtype class then pass
elif data_type in const.DICT_GDAL_TO_NP:
# Maybe convert to the integer form?
pass
# if gdal dtype name, then convert to data type integer alias
elif isinstance(data_type,str):
data_type = gdal.GetDataTypeByName(data_type)
# If the data array is 2D, add a dimension for a single band
if data_np_array.ndim == 2:
data_np_array = data_np_array[np.newaxis, :, :]
elif data_np_array.ndim == 3:
pass
else:
raise ValueError("This can't handle Arrays larger than three "
"dimensions.")
(n_bands,y_size,x_size) = data_np_array.shape
# Create driver and data set object
driver = gdal.GetDriverByName(gdal_driver_name)
# Check that the driver supports the data_np_array.dtype
dstr = gdal.GetDataTypeName(const.DICT_NP_TO_GDAL[data_np_array.dtype])
dlist = driver.GetMetadata()['DMD_CREATIONDATATYPES'].split()
if not dstr in dlist:
raise ValueError("The data type of the provided numpy array is not "
"supported by the requested file format.")
if options:
if not isinstance(options,list):
raise ValueErrror("The options list is malformed. It should be a "
"list of strings")
dst_ds = driver.Create(new_file_name, x_size, y_size, n_bands,
data_type, options)
# Set projection and geotransform
dst_ds.SetGeoTransform(gdal_geo_t)
dst_ds.SetProjection(gdal_projection)
### Write the new data
# Set nans to the original No Data Value
if NDV is not None:
data_np_array[np.isnan(data_np_array)] = NDV
data_np_array[np.isinf(data_np_array)] = NDV
### Other???
# Write the bands
for b in range(n_bands):
dst_ds.GetRasterBand(b+1).WriteArray( data_np_array[b,:,:] )
if NDV is not None:
dst_ds.GetRasterBand(b+1).SetNoDataValue( NDV )
# Once we're done, close properly the dataset
dst_ds = None
# Return the new file name in case it was changed due to vrt fallback.
return new_file_name
def get_img_stretch_vals(imgfname_or_gdalobj,stretch=[0.02,0.98],approx_ok=True):
"""Read image, mask very small values, and return max and min
imgfname_or_gdalobj: Either a file name string of a gdal object
stretch: stretch range to return in a list
approx_ok: Is it ok to approximate the max/min values (bool)
hist_nbins: How many bins the histogram is made of to find the
appropriate cut value (more bins is more accurate,
but slower).
"""
# ToDo Add tests for get_img_strech_vals
### Cases to consider:
# Lots of fill 0 values with much larger data
# Lots of fill 0 values with data near zero
# Large integers
# Floating point 0-1
# Floating point larger than 1
# Negative values to postivie values (a la NDVI)
#... basically this needs to respect the No Data Value and let the math fly
if isinstance(imgfname_or_gdalobj,str):
f = gdal.Open(imgfname_or_gdalobj)
elif isinstance(imgfname_or_gdalobj,gdal.Dataset):
f = imgfname_or_gdalobj
else:
raise ValueError("The variable passed in as imgfname_or_gdalobj is "
"neither an image filename string or a gdal object.")
# Setup for query of each band
outofrange = 0 # value zero means no, don't include out of range
hist_nbins = 1000
if approx_ok:
approx = 1
else:
approx = 0
for i in range(1,f.RasterCount+1):
# ToDo Add respect for NoDataValue
# Get gdal band object
b = f.GetRasterBand(i)
# Get band max/min
tmp = b.ComputeRasterMinMax(approx)
imgmax = tmp[1]
imgmin = tmp[0]
# Use max/min to get histogram
hist = b.GetHistogram(imgmin,imgmax,hist_nbins,outofrange,approx_ok)
hist = np.asarray(hist)
bins = np.linspace(imgmin, imgmax, len(hist))
# Get stretch values
cdf = hist.cumsum()
top = cdf.max() * stretch[1]
bottom = cdf.max() * stretch[0]
topcut = bins[np.where(cdf > top)[0][0]]
bottomcut = bins[np.where(cdf < bottom)[0][-1]]
return (bottomcut,topcut)
class GeoSet(Sequence):
def __init__(self,*args,**kwargs):
self.imgs = []
for x in args:
self.imgs.append(GeoImage(x))
self._set_set_meta()
def __repr__(self):
svals = [x.meta.class_name for x in self.imgs]
return ', '.join(svals)
def __len__(self):
return len(self.imgs)
def __getitem__(self,i):
return self.imgs[i]
def _set_set_meta(self):
# TBD - time difference, angle difference, etc.
pass
def get_data(self):
pass
