Python osgeo.gdal.GA_ReadOnly() Examples

def GetGeoInfo(FileName):

    if exists(FileName) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')    
    
    
    SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
    if SourceDS == None:
        raise Exception("Unable to read the data file")
    
    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)
    
    return NDV, xsize, ysize, GeoT, Projection, DataType
#==============================================================================

#==============================================================================
# Function to read the original file's projection: 

def read_image():
    """ Read image ``f`` and return ``np.array`` of image values

    Image will be (nband, nrow, ncol)
    """
    def _read_image(f):
        ds = gdal.Open(f, gdal.GA_ReadOnly)
        dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
            ds.GetRasterBand(1).DataType)
        nrow, ncol, nband = ds.RasterYSize, ds.RasterYSize, ds.RasterCount
        img = np.empty((nband, nrow, ncol), dtype=dtype)
        for i_b in range(nband):
            b = ds.GetRasterBand(i_b + 1)
            img[i_b, ...] = b.ReadAsArray()
        return img
    return _read_image 

def read_naip(file_path, bands_to_use):
    """
    Read in a NAIP, based on www.machinalis.com/blog/python-for-geospatial-data-processing.

    Bands_to_use is an array like [0,0,0,1], designating whether to use each band (R, G, B, IR).
    """
    raster_dataset = gdal.Open(file_path, gdal.GA_ReadOnly)

    bands_data = []
    index = 0
    for b in range(1, raster_dataset.RasterCount + 1):
        band = raster_dataset.GetRasterBand(b)
        if bands_to_use[index] == 1:
            bands_data.append(band.ReadAsArray())
        index += 1
    bands_data = numpy.dstack(bands_data)

    return raster_dataset, bands_data 

def tag_with_locations(test_images, predictions, tile_size, state_abbrev):
    """Combine image data with label data, so info can be rendered in a map and list UI.

    Add location data for convenience too.
    """
    combined_data = []
    for idx, img_loc_tuple in enumerate(test_images):
        raster_filename = img_loc_tuple[2]
        raster_dataset = gdal.Open(os.path.join(NAIP_DATA_DIR, raster_filename), gdal.GA_ReadOnly)
        raster_tile_x = img_loc_tuple[1][0]
        raster_tile_y = img_loc_tuple[1][1]
        ne_lat, ne_lon = pixel_to_lat_lon_web_mercator(raster_dataset, raster_tile_x +
                                                       tile_size, raster_tile_y)
        sw_lat, sw_lon = pixel_to_lat_lon_web_mercator(raster_dataset, raster_tile_x,
                                                       raster_tile_y + tile_size)
        certainty = predictions[idx][0]
        formatted_info = {'certainty': certainty, 'ne_lat': ne_lat, 'ne_lon': ne_lon,
                          'sw_lat': sw_lat, 'sw_lon': sw_lon, 'raster_tile_x': raster_tile_x,
                          'raster_tile_y': raster_tile_y, 'raster_filename': raster_filename,
                          'state_abbrev': state_abbrev, 'country_abbrev': 'USA'
                          }
        combined_data.append(formatted_info)
    return combined_data 

def get_image_attribute(image_filename):
    """ Use GDAL to open image and return some attributes

    Args:
        image_filename (str): image filename

    Returns:
        tuple: nrow (int), ncol (int), nband (int), NumPy datatype (type)
    """
    try:
        image_ds = gdal.Open(image_filename, gdal.GA_ReadOnly)
    except Exception as e:
        logger.error('Could not open example image dataset ({f}): {e}'
                     .format(f=image_filename, e=str(e)))
        raise

    nrow = image_ds.RasterYSize
    ncol = image_ds.RasterXSize
    nband = image_ds.RasterCount
    dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
        image_ds.GetRasterBand(1).DataType)

    return (nrow, ncol, nband, dtype) 

def open_inumpyut(self,vrt):
        if vrt:
            self.in_ds = gdal.Open(vrt, gdal.GA_ReadOnly)
        else:
            raise Exception("No vrt file was specified")
            
        if self.options.verbose:
            print("Inumpyut file:", "( %sP x %sL - %s bands)" % (self.in_ds.RasterXSize, self.in_ds.RasterYSize, self.in_ds.RasterCount))

        if not self.in_ds:
            # Note: GDAL prints the ERROR message too
            self.error("It is not possible to open the inumpyut file '%s'." % vrt )
            
        if self.in_ds.RasterCount == 0:
            self.error( "Inumpyut file '%s' has no raster band" % vrt )
        
        self.out_ds = self.in_ds
        
        # Get alpha band (either directly or from NODATA value)
        self.alphaband = self.out_ds.GetRasterBand(1).GetMaskBand()
        if (self.alphaband.GetMaskFlags() & gdal.GMF_ALPHA) or self.out_ds.RasterCount==4 or self.out_ds.RasterCount==2:
            self.dataBandsCount = self.out_ds.RasterCount - 1
        else:
            self.dataBandsCount = self.out_ds.RasterCount

    # ------------------------------------------------------------------------- 

def read_pixel_timeseries(files, px, py):
    """ Returns NumPy array containing timeseries values for one pixel

    Args:
        files (list): List of filenames to read from
        px (int): Pixel X location
        py (int): Pixel Y location

    Returns:
        np.ndarray: Array (nband x n_images) containing all timeseries data
            from one pixel
    """
    nrow, ncol, nband, dtype = get_image_attribute(files[0])

    if px < 0 or px >= ncol or py < 0 or py >= nrow:
        raise IndexError('Row/column {r}/{c} is outside of image '
                         '(nrow/ncol: {nrow}/{ncol})'.
                         format(r=py, c=px, nrow=nrow, ncol=ncol))

    Y = np.zeros((nband, len(files)), dtype=dtype)

    for i, f in enumerate(files):
        ds = gdal.Open(f, gdal.GA_ReadOnly)
        for b in range(nband):
            Y[b, i] = ds.GetRasterBand(b + 1).ReadAsArray(px, py, 1, 1)

    return Y 

def get_image_metadata(img_name):
    dataset = gdal.Open(img_name, gdal.GA_ReadOnly)
    metadata = dataset.GetMetadata()
    rpc_data = dataset.GetMetadata('RPC')
    date_time = metadata['NITF_IDATIM']
    year = int(date_time[0:4])
    month = int(date_time[4:6])
    day = int(date_time[6:8])
    return metadata, month, day


# epipolar rectify an image pair and save the results 

def read_metadata(self, img_name):

        # Read the metadata
        dataset = gdal.Open(img_name, gdal.GA_ReadOnly)
        metadata = dataset.GetMetadata()
        rpc_data = dataset.GetMetadata('RPC')

        # read image and get size
        img = dataset.ReadAsArray()
        self.rows = img.shape[1]
        self.columns = img.shape[2]
        img = None

        # Extract RPC metadata fields
        self.lon_off = float(rpc_data['LONG_OFF'])
        self.lon_scale = float(rpc_data['LONG_SCALE'])
        self.lat_off = float(rpc_data['LAT_OFF'])
        self.lat_scale = float(rpc_data['LAT_SCALE'])
        self.height_off = float(rpc_data['HEIGHT_OFF'])
        self.height_scale = float(rpc_data['HEIGHT_SCALE'])
        self.line_off = float(rpc_data['LINE_OFF'])
        self.line_scale = float(rpc_data['LINE_SCALE'])
        self.samp_off = float(rpc_data['SAMP_OFF'])
        self.samp_scale = float(rpc_data['SAMP_SCALE'])
        self.line_num_coeff = np.asarray(rpc_data['LINE_NUM_COEFF'].split(), dtype=np.float64)
        self.line_den_coeff = np.asarray(rpc_data['LINE_DEN_COEFF'].split(), dtype=np.float64)
        self.samp_num_coeff = np.asarray(rpc_data['SAMP_NUM_COEFF'].split(), dtype=np.float64)
        self.samp_den_coeff = np.asarray(rpc_data['SAMP_DEN_COEFF'].split(), dtype=np.float64)

        # Close the dataset
        dataset = None

    # Forward project normalized WGS84 array to image coordinates
    # Rows and columns are computed with separate calls and coefficient arrays
    # Not currently checking the denominator for zero 

def getProjectionSystem(self):
        '''
        Testing with Greenland.
        '''
        if not self.demname:
            raise Exception('At least the DEM parameter must be set for geogrid')

        from osgeo import gdal, osr
        ds = gdal.Open(self.demname, gdal.GA_ReadOnly)
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection())
        srs.AutoIdentifyEPSG()
        ds = None
#        pdb.set_trace()

        if srs.IsGeographic():
            epsgstr = srs.GetAuthorityCode('GEOGCS')
        elif srs.IsProjected():
            epsgstr = srs.GetAuthorityCode('PROJCS')
        elif srs.IsLocal():
            raise Exception('Local coordinate system encountered')
        else:
            raise Exception('Non-standard coordinate system encountered')
        if not epsgstr:  #Empty string->use shell command gdalsrsinfo for last trial
            cmd = 'gdalsrsinfo -o epsg {0}'.format(self.demname)
            epsgstr = subprocess.check_output(cmd, shell=True)
#            pdb.set_trace()
            epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0]
#            pdb.set_trace()
        if not epsgstr:  #Empty string
            raise Exception('Could not auto-identify epsg code')
#        pdb.set_trace()
        epsgcode = int(epsgstr)
#        pdb.set_trace()
        return epsgcode 

def getProjectionSystem(self, filename):
        '''
        Testing with Greenland.
        '''
        if not filename:
            raise Exception('File {0} does not exist'.format(filename))

        from osgeo import gdal, osr
        ds = gdal.Open(filename, gdal.GA_ReadOnly)
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection())
        srs.AutoIdentifyEPSG()
        ds = None
#        pdb.set_trace()

        if srs.IsGeographic():
            epsgstr = srs.GetAuthorityCode('GEOGCS')
        elif srs.IsProjected():
            epsgstr = srs.GetAuthorityCode('PROJCS')
        elif srs.IsLocal():
            raise Exception('Local coordinate system encountered')
        else:
            raise Exception('Non-standard coordinate system encountered')
        if not epsgstr:  #Empty string->use shell command gdalsrsinfo for last trial
            cmd = 'gdalsrsinfo -o epsg {0}'.format(filename)
            epsgstr = subprocess.check_output(cmd, shell=True)
#            pdb.set_trace()
            epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0]
#            pdb.set_trace()
        if not epsgstr:  #Empty string
            raise Exception('Could not auto-identify epsg code')
#        pdb.set_trace()
        epsgcode = int(epsgstr)
#        pdb.set_trace()
        return epsgcode 

def __init__(self, tiff_details, symmetry=(False, False, True)):
        """

        Parameters
        ----------
        tiff_details : TiffDetails
        symmetry : Tuple[bool]
        """
        if isinstance(tiff_details, str):
            tiff_details = TiffDetails(tiff_details)
        if not isinstance(tiff_details, TiffDetails):
            raise TypeError('GdalTiffChipper input argument must be a filename '
                            'or TiffDetails object.')

        self._tiff_details = tiff_details
        # initialize our dataset - NB: this should close gracefully on garbage collection
        self._data_set = gdal.Open(tiff_details.file_name, gdal.GA_ReadOnly)
        if self._data_set is None:
            raise ValueError(
                'GDAL failed with unspecified error in opening file {}'.format(tiff_details.file_name))
        # get data_size information
        data_size = (self._data_set.RasterYSize, self._data_set.RasterXSize)
        self._bands = self._data_set.RasterCount
        # TODO: get data_type information - specifically how does complex really work?
        complex_type = ''
        super(GdalTiffChipper, self).__init__(data_size, symmetry=symmetry, complex_type=complex_type)
        # 5.) set up our virtual array using GetVirtualMemArray
        try:
            self._virt_array = self._data_set.GetVirtualMemArray()
        except Exception:
            logging.error(
                msg="There has been some error using the gdal method GetVirtualMemArray(). "
                    "Consider falling back to the base sarpy tiff reader implementation (use_gdal=False)")
            raise
        # TODO: this does not generally work should we clunkily fall back to dataset.band.ReadAsArray()?
        #   This doesn't support slicing... 

def _dataset_from_dir(folder):
        files = set(f for f in os.scandir(folder) if f.is_file())
        for f in files:
            if op.splitext(f.name)[-1] == '':
                break
        else:
            raise ImportError('could not load dataset from %s' % folder)

        return gdal.Open(f.path, gdal.GA_ReadOnly) 

def _getLOS(self, filename, component):
        path = op.dirname(filename)
        fn = glob.glob(op.join(path, component))
        if len(fn) != 1:
            raise ImportError('Cannot find LOS vector file %s!' % component)

        dataset = gdal.Open(fn[0], gdal.GA_ReadOnly)
        return self._readBandData(dataset) 

def toGDAL(self):
		'''Get GDAL dataset'''
		return gdal.Open(self.path, gdal.GA_ReadOnly) 

def PlotRaster(RasterFile, Map, alpha=1.):

    """
    Description goes here...

    MDH
    """

    print "Plotting raster..."

    #Read data
    gdata = gdal.Open(RasterFile, gdal.GA_ReadOnly)
    geo = gdata.GetGeoTransform()
    Data = gdata.ReadAsArray()

    # make topodat a masked array, masking values lower than sea level.
    Data = ma.masked_where(Data < 0, Data)

    #setup meshgrid for raster plotting
    xres = geo[1]
    yres = geo[5]
    xmin = geo[0] + xres * 0.5
    xmax = geo[0] + (xres * gdata.RasterXSize) - xres * 0.5
    ymin = geo[3] + (yres * gdata.RasterYSize) + yres * 0.5
    ymax = geo[3] - yres * 0.5

    x,y = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
    x,y = Map(x,y)
    Map.pcolormesh(x, y, Data.T, cmap=plt.cm.Greys, alpha=alpha) 

def _fromGDAL(self):
		'''Use GDAL to extract raster infos and init'''
		if self.path is None or not self.fileExists:
			raise IOError("Cannot find file on disk")
		ds = gdal.Open(self.path, gdal.GA_ReadOnly)
		self.size = xy(ds.RasterXSize, ds.RasterYSize)
		self.format = ds.GetDriver().ShortName
		if self.format in ['JP2OpenJPEG', 'JP2ECW', 'JP2KAK', 'JP2MrSID'] :
			self.format = 'JPEG2000'
		self.nbBands = ds.RasterCount
		b1 = ds.GetRasterBand(1) #first band (band index does not count from 0)
		self.noData = b1.GetNoDataValue()
		ddtype = gdal.GetDataTypeName(b1.DataType)#Byte, UInt16, Int16, UInt32, Int32, Float32, Float64
		if ddtype == "Byte":
			self.dtype = 'uint'
			self.depth = 8
		else:
			self.dtype = ddtype[0:len(ddtype)-2].lower()
			self.depth = int(ddtype[-2:])
		#Get Georef
		self.georef = GeoRef.fromGDAL(ds)
		#Close (gdal has no garbage collector)
		ds, b1 = None, None

	#######################################
	# Dynamic properties
	####################################### 

def GetGeoInfo(FileName):
    """This gets information from the raster file using gdal

    Args:
        FileName (str): The filename (with path and extension) of the raster.

    Return:
        float: A vector that contains:
            * NDV: the nodata values
            * xsize: cellsize in x direction
            * ysize: cellsize in y direction
            * GeoT: the tranform (a string)
            * Projection: the Projection (a string)
            * DataType: The type of data (an int explaing the bits of each data element)

    Author: SMM
    """


    if exists(FileName) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')


    SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
    if SourceDS == None:
        raise Exception("Unable to read the data file")

    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)

    return NDV, xsize, ysize, GeoT, Projection, DataType
#==============================================================================

#==============================================================================
# This gets the UTM zone, if it exists 

def loadColisions():
	global obstacle_image, heightMapFile, array, width, height, minHeight, maxHeight, minColisionHeight, maxColisionHeight
	global barrierN, barrierS, barrierE, barrierW, barrierNE, barrierNW, barrierSE, barrierSW, barrier
	# loading height map image
	dataset = gdal.Open(heightMapFile, gdal.GA_ReadOnly)
	try:
		srcband = dataset.GetRasterBand(1)
	except RuntimeError, e:
        	print 'No band %i found' % band_num
        	print e
	        sys.exit(1) 

def _init_attrs(self, filenames):
        self.filenames = filenames
        self.datasets = [gdal.Open(f, gdal.GA_ReadOnly) for
                         f in self.filenames]

        self.n_image = len(filenames)
        self.n_band = self.datasets[0].RasterCount
        self.n_col = int(self.datasets[0].RasterXSize)
        self.datatype = gdal_array.GDALTypeCodeToNumericTypeCode(
            self.datasets[0].GetRasterBand(1).DataType)

        self.dataset_bands = [
            [ds.GetRasterBand(i + 1) for i in range(self.n_band)]
            for ds in self.datasets
        ] 

def test_unw_contains_same_data_as_numpy_array(self):
        from datetime import time
        temp_unw = tempfile.mktemp(suffix='.unw')
        temp_tif = tempfile.mktemp(suffix='.tif')

        # setup some header files for use in write_geotif
        dem_header_file = common.SML_TEST_DEM_HDR_GAMMA
        dem_header = gamma.parse_dem_header(dem_header_file)

        header = gamma.parse_epoch_header(
            os.path.join(common.SML_TEST_GAMMA, '20060828_slc.par'))
        header.update(dem_header)

        # insert some dummy data so we are the dem in write_fullres_geotiff is not
        # not activated and ifg write_fullres_geotiff operation works
        header[ifc.PYRATE_TIME_SPAN] = 0
        header[ifc.SLAVE_DATE] = 0
        header[ifc.DATA_UNITS] = 'degrees'
        header[ifc.DATA_TYPE] = ifc.ORIG
        header[ifc.SLAVE_TIME] = time(10)

        # now create aritrary data
        data = np.random.rand(dem_header[ifc.PYRATE_NROWS],
                              dem_header[ifc.PYRATE_NCOLS])

        # convert numpy array to .unw
        shared.write_unw_from_data_or_geotiff(geotif_or_data=data,
                                              dest_unw=temp_unw,
                                              ifg_proc=1)
        # convert the .unw to geotif
        shared.write_fullres_geotiff(header=header, data_path=temp_unw,
                                     dest=temp_tif, nodata=np.nan)

        # now compare geotiff with original numpy array
        ds = gdal.Open(temp_tif, gdal.GA_ReadOnly)
        data_lv_theta = ds.ReadAsArray()
        ds = None
        np.testing.assert_array_almost_equal(data, data_lv_theta)
        try:
            os.remove(temp_tif)
        except PermissionError:
            print("File opened by another process.")

        try:
            os.remove(temp_unw)
        except PermissionError:
            print("File opened by another process.") 

def _set_special_model_input(self, field, dims):
        ''' Special processing for setting certain input fields. Only relevant for image processing
        mode.

        Parameters
        ----------
        field : string
            The name of the input field for which the special processing is needed.
        dims : int list
            The dimensions of the output parameter array.

        Returns
        -------
        success : boolean
            True is the parameter was succefully set, false otherwise.
        array : float array
            The set parameter array.
        '''
        if field in ["flux_LR", "flux_LR_ancillary"]:
            # Low resolution data in case disaggregation is to be used.
            inputs = {}
            fid = gdal.Open(self.p[field], gdal.GA_ReadOnly)
            if fid is None:
                print("ERROR: Low resolution data for disaggregation is not avaiable.")
                return False, None
            prj_LR = fid.GetProjection()
            geo_LR = fid.GetGeoTransform()
            subset = []
            if "subset" in self.p:
                subset, geo_LR = self._get_subset(self.p["subset"], prj_LR, geo_LR)
                inputs[field] = fid.GetRasterBand(1).ReadAsArray(subset[0],
                                                                 subset[1],
                                                                 subset[2],
                                                                 subset[3])
            else:
                inputs[field] = fid.GetRasterBand(1).ReadAsArray()
            inputs['scale'] = [geo_LR, prj_LR, self.geo, self.prj]
            success = True
        else:
            success = False
            inputs = None

        return success, inputs 

def _set_param_array(self, parameter, dims, band=1):
        '''Set model input parameter as an array.

        Parameters
        ----------
        parameter : float or string
            If float it should be the value of the parameter. If string of a number it is also
            the value of the parameter. Otherwise it is the name of the parameter and the value, or
            path to raster file which contains the values, is read from the parameters dictionary.
        dims : int list
            The dimensions of the output parameter array.
        band : int (default = 1)
            Band (in GDAL convention) of raster file to be read, if parameter is to be read from a
            raster file.

        Returns
        -------
        success : boolean
            True is the parameter was succefully set, false otherwise.
        array : float array
            The set parameter array.
        '''

        success = True
        array = None

        # See if the parameter is a number
        try:
            array = np.zeros(dims) + float(parameter)
            return success, array
        except ValueError:
            pass

        # Otherwise see if the parameter is a parameter name
        try:
            inputString = self.p[parameter]
        except KeyError:
            success = False
            return success, array
        # If it is then get the value of that parameter
        try:
            array = np.zeros(dims) + float(inputString)
        except ValueError:
            try:
                fid = gdal.Open(inputString, gdal.GA_ReadOnly)
                if self.subset:
                    array = fid.GetRasterBand(band).ReadAsArray(self.subset[0],
                                                                self.subset[1],
                                                                self.subset[2],
                                                                self.subset[3])
                else:
                    array = fid.GetRasterBand(band).ReadAsArray()
            except AttributeError:
                print("%s image not present for parameter %s"%(inputString, parameter))
                success = False
            finally:
                fid = None

        return success, array 

def read(self, filename, **kwargs):
        dataset = gdal.Open(filename, gdal.GA_ReadOnly)
        georef = dataset.GetGeoTransform()

        llLon = georef[0]
        llLat = georef[3] + dataset.RasterYSize * georef[5]

        c = self.container

        c.frame.spacing = 'degree'
        c.frame.llLat = llLat
        c.frame.llLon = llLon
        c.frame.dE = georef[1]
        c.frame.dN = abs(georef[5])

        displacement = self._readBandData(dataset)
        c.displacement = -displacement / (4*num.pi) * LAMBDA_SENTINEL

        try:
            los_n = self._getLOS(filename, '*.geo.N.tif')
            los_e = self._getLOS(filename, '*.geo.E.tif')
            los_u = self._getLOS(filename, '*.geo.U.tif')
        except ImportError:
            self._log.warning(
                'Cannot find LOS angle files *.geo.[NEU].tif,'
                ' using static Sentinel-1 descending angles.')

            heading = 83.
            incident = 50.

            un = num.sin(d2r*incident) * num.cos(d2r*heading)
            ue = num.sin(d2r*incident) * num.sin(d2r*heading)
            uz = num.cos(d2r*incident)

            los_n = num.full_like(c.displacement, un)
            los_e = num.full_like(c.displacement, ue)
            los_u = num.full_like(c.displacement, uz)

        c.phi = num.arctan2(los_n, los_e)
        c.theta = num.arcsin(los_u)

        c.meta.title = dataset.GetDescription()

        return c 

def CheckNoData(FileName):
    """This looks through the head file of an ENVI raster and if it doesn't find the nodata line it rewrites the file to include the nodata line.

    Args:
        FileName (str): The filename (with path and extension) of the raster.

    Return:
        int: The total number of pixels (although what it is really doing is updating the header file. The return is just to check if it is working and yes I know this is stupid. )

    Author: SMM
    """


    if exists(FileName) is False:
        raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')

    # read the file, and check if there is a no data value
    SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
    if SourceDS == None:
        raise Exception("Unable to read the data file")
    NoDataValue = SourceDS.GetRasterBand(1).GetNoDataValue()

    print("In the check nodata routine. Nodata is: ")
    print(NoDataValue)

    if NoDataValue == None:
        print("This raster does not have no data. Updating the header file")
        header_name = FileName[:-4]
        header_name = header_name+".hdr"

        # read the header
        if exists(header_name) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + header_name + '\'')
        else:
            this_file = open(header_name, 'r')
            lines = this_file.readlines()
            lines.append("data ignore value = -9999")
            this_file.close()

            this_file = open(header_name, 'w')
            for item in lines:
                this_file.write("%s" % item)        # no newline since a newline command character comes with the lines

            this_file.close()

    NDV, xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)

    return xsize*ysize

#==============================================================================


#============================================================================== 

def read_image(image_filename, bands=None, dtype=None):
    """ Return raster image bands as a sequence of NumPy arrays

    Args:
        image_filename (str): Image filename
        bands (iterable, optional): A sequence of bands to read from image.
            If `bands` is None, function returns all bands in raster. Note that
            bands are indexed on 1 (default: None)
        dtype (np.dtype): NumPy datatype to use for image bands. If `dtype` is
            None, arrays are kept as the image datatype (default: None)

    Returns:
        list: list of NumPy arrays for each band specified

    Raises:
        IOError: raise IOError if bands specified are not contained within
            raster
        RuntimeError: raised if GDAL encounters errors
    """
    try:
        ds = gdal.Open(image_filename, gdal.GA_ReadOnly)
    except:
        logger.error('Could not read image {i}'.format(i=image_filename))
        raise

    if bands:
        if not all([b in range(1, ds.RasterCount + 1) for b in bands]):
            raise IOError('Image {i} ({n} bands) does not contain bands '
                          'specified (requested {b})'.
                          format(i=image_filename, n=ds.RasterCount, b=bands))
    else:
        bands = range(1, ds.RasterCount + 1)

    if not dtype:
        dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
            ds.GetRasterBand(1).DataType)

    output = []
    for b in bands:
        output.append(ds.GetRasterBand(b).ReadAsArray().astype(dtype))

    return output 

def changemap(ctx, map_type, start_date, end_date, output,
              root, result, image, date_frmt, ndv, gdal_frmt, out_date_frmt,
              warn_on_empty, magnitude):
    """
    Examples: TODO
    """
    gdal_frmt = str(gdal_frmt)  # GDAL GetDriverByName doesn't work on Unicode

    frmt = '%Y%j'
    start_txt, end_txt = start_date.strftime(frmt), end_date.strftime(frmt)
    start_date, end_date = start_date.toordinal(), end_date.toordinal()

    try:
        image_ds = gdal.Open(image, gdal.GA_ReadOnly)
    except:
        logger.error('Could not open example image for reading')
        raise

    if map_type in ('first', 'last'):
        changemap, magnitudemap, magnitude_indices = get_change_date(
            start_date, end_date, result, image_ds,
            first=map_type == 'first', out_format=out_date_frmt,
            magnitude=magnitude,
            ndv=ndv, pattern='yatsm_r*', warn_on_empty=warn_on_empty
        )

        band_names = ['ChangeDate_s{s}-e{e}'.format(s=start_txt, e=end_txt)]
        write_output(changemap, output, image_ds, gdal_frmt, ndv,
                     band_names=band_names)

        if magnitudemap is not None:
            band_names = (['Magnitude Index {}'.format(i) for
                           i in magnitude_indices])
            name, ext = os.path.splitext(output)
            output = name + '_mag' + ext
            write_output(magnitudemap, output, image_ds, gdal_frmt, ndv,
                         band_names=band_names)

    elif map_type == 'num':
        changemap = get_change_num(
            start_date, end_date, result, image_ds,
            ndv=ndv, pattern='yatsm_r*', warn_on_empty=warn_on_empty
        )

        band_names = ['NumChanges_s{s}-e{e}'.format(s=start_txt, e=end_txt)]
        write_output(changemap, output, image_ds, gdal_frmt, ndv,
                     band_names=band_names)

    image_ds = None