Python osgeo.gdal.Open() Examples

The following are code examples for showing how to use osgeo.gdal.Open(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: lidar   Author: giswqs   File: slicing.py    MIT License 9 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)
    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# convert images in a selected folder to shapefiles 
Example 2
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def testRaster(dlocation, bandnum=1):
        """
        @return: a GDAL dataset object
        """
        success = True
        try:
            try:
                dataset = gdal.Open(str(dlocation), gdalconst.GA_ReadOnly)
            except Exception, e:
                raise LMError(currargs='Unable to open dataset {} with GDAL ({})'
                                    .format(dlocation, str(e)))
            try:
                band = dataset.GetRasterBand(1)
            except Exception, e:
                raise LMError(currargs='No band {} in dataset {} ({})'
                                    .format(band, dlocation, str(e))) 
Example 3
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def isValidDataset(self, dlocation=None):
        """
        @summary: Checks to see if the dataset at self.dlocations is a valid 
                     vector readable by OGR.
        @return: True if dataset is a valid OGR dataset; False if not
        """
        valid = False
        if dlocation is None:
            dlocation = self._dlocation
        if (dlocation is not None
             and (os.path.isdir(dlocation) 
                    or os.path.isfile(dlocation)) ):
            try:
                ds = ogr.Open(dlocation)
                ds.GetLayer(0)
            except Exception, e:
                pass
            else:
                valid = True 
Example 4
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def getShapefileRowHeaders(cls, shapefileFilename):
        """
        @summary: Get a (sorted by feature id) list of row headers for a shapefile
        @todo: Move this to a common Vector class in LmCommon or LmBackend and 
                 use with LmCompute.  This is a rough copy of what is now used for 
                 rasterIntersect.
        """
        ogr.RegisterAll()
        drv = ogr.GetDriverByName(LMFormat.getDefaultOGR().driver)
        ds = drv.Open(shapefileFilename)
        lyr = ds.GetLayer(0)
        
        rowHeaders = []
        
        for j in range(lyr.GetFeatureCount()):
            curFeat = lyr.GetFeature(j)
            siteIdx = curFeat.GetFID()
            x, y = curFeat.geometry().Centroid().GetPoint_2D()
            rowHeaders.append((siteIdx, x, y))
            
        return sorted(rowHeaders)

# ............................................... 
Example 5
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def testVector(dlocation, driver=LMFormat.getDefaultOGR().driver):
        goodData = False
        featCount = 0
        if dlocation is not None and os.path.exists(dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(driver)
            try:
                ds = drv.Open(dlocation)
            except Exception, e:
                goodData = False
            else:
                try:
                    slyr = ds.GetLayer(0)
                except Exception, e:
                    goodData = False
                else: 
Example 6
Project: core   Author: lifemapper   File: extract_environment_values.py    GNU General Public License v3.0 6 votes vote down vote up
def get_metrics_for_layer(points, layer_filename, metricFunctions):
   """
   @summary: Get layer values for each point and then generate metrics
   """
   ds = gdal.Open(layer_filename)
   band = ds.GetRasterBand(1)
   data = np.array(band.ReadAsArray())
   gt = ds.GetGeoTransform()
   nodataVal = band.GetNoDataValue()
   
   values = []
   
   for x, y in points:
      px = int((x - gt[0]) / gt[1])
      py = int((y - gt[3]) / gt[5])
      try:
         val = data[py, px]
         if not is_close(val, nodataVal):
            values.append(data[py, px])
         else:
            print 'Could not append value at ({}, {}): {}'.format(px, py, val)
      except Exception, e:
            print 'Could not append value at ({}, {}): {}'.format(px, py, str(e)) 
Example 7
Project: core   Author: lifemapper   File: extract_environment_values.py    GNU General Public License v3.0 6 votes vote down vote up
def get_point_xys(points_filename, removeDuplicateLocations=True):
   """
   @summary: Get x,y pairs for each point in a shapefile
   """
   points = []
   
   ds = ogr.Open(points_filename)
   lyr = ds.GetLayer()
   
   for feat in lyr:
      geom = feat.GetGeometryRef()
      points.append((geom.GetX(), geom.GetY()))
   
   if removeDuplicateLocations:
      points = list(set(points))
      
   return points

# ............................................................................. 
Example 8
Project: core   Author: lifemapper   File: geotools.py    GNU General Public License v3.0 6 votes vote down vote up
def openDataSet(self, varpattern=None, updateable=False):
      """
      @summary: Open (or re-open) the dataset.
      @param varpattern: string to match for the final portion of a subdataset
             name, for those datasets where the data of interest is one of 
             multiple subdatasets.
      @param updateable: False if open in Read-Only mode, True if writeable. 
      """
      try:
         if updateable:
            self._dataset = None
            self._band = None
            self._dataset = gdal.Open(str(self.dlocation), gdalconst.GA_Update)
         elif self._dataset is None:
            self._dataset = gdal.Open(str(self.dlocation), gdalconst.GA_ReadOnly)
      except Exception,e:
         print('Exception raised trying to open layer=%s\n%s'
               % (self.dlocation, str(e)))
         raise LMError(['Unable to open %s' % self.dlocation]) 
Example 9
Project: spacesense   Author: spacesense-ai   File: datasets.py    GNU Lesser General Public License v3.0 6 votes vote down vote up
def __init__(self, folder_path):
        self.folder_path = folder_path
        self.level = self.folder_path.split('/')[-1].split('.')[0].split('_')[1][-2:]
        if self.level == '2A':
            self.band_files = glob(self.folder_path + '/GRANULE' + '/*/IM*/*10*/*B*')
            for ext in ('*B05*', '*B06*', '*B07*', '*B8A*', '*B11*', '*B12*'):
                self.band_files.extend(glob(os.path.join(self.folder_path + '/GRANULE' + '/*/IM*/*20*/', ext)))
            for ext in ('*B01*', '*B09*', '*B10*'):
                self.band_files.extend(glob(os.path.join(self.folder_path + '/GRANULE' + '/*/IM*/*60*/', ext)))
            self.band_files = np.array(sorted(self.band_files))
            self.band_details = [file.split('_')[-2] for file in self.band_files]

        else:
            self.band_files = np.array(sorted(glob(self.folder_path + '/GRANULE' + '/*/IM*/*B*.jp2')))
            self.band_details = [file.split('_')[-1].split('.')[0] for file in self.band_files]
        im = gdal.Open(self.band_files[1])
        # arf_base = im.ReadAsArray()
        # self.img_shp = arf_base.shape
        self.meta_data = im.GetMetadata()
        self.num_bands = len(self.band_files)
        self.save_folder = self.folder_path
        self.AOI = None
        self.data = None
        self.data_dictionary = None 
Example 10
Project: geophys_utils   Author: GeoscienceAustralia   File: _gdal_grid_utils.py    Apache License 2.0 6 votes vote down vote up
def get_gdal_wcs_dataset(wcs_url):
    clean_url = re.match('http[^?]+', wcs_url).group(0)
    temp_xml_path = os.path.join(tempfile.gettempdir(), re.sub('\W', '_', clean_url) + '.xml')
    
    wcs = WebCoverageService(wcs_url, version='1.0.0')
    variable_name = list(wcs.contents.keys())[0] # Only work with first variable
    
    xml_string = '''<WCS_GDAL>
  <ServiceURL>%s?</ServiceURL>
  <CoverageName>%s</CoverageName>
</WCS_GDAL>''' % (clean_url, variable_name)

    temp_xml_file = open(temp_xml_path, 'w') 
    temp_xml_file.write(xml_string)    
    temp_xml_file.close()                        

    return gdal.Open(temp_xml_path) 
Example 11
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    MIT License 6 votes vote down vote up
def getNoDataValue(rasterfn):
    """This gets the nodata value from the raster

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

    Returns:
        float: nodatavalue; the nodata value

    Author: SMM
    """
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.GetNoDataValue()
#==============================================================================

#============================================================================== 
Example 12
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    MIT License 6 votes vote down vote up
def setNoDataValue(rasterfn):
    """This sets the nodata value from the raster

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

    Returns:
        None

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.SetNoDataValue()
#==============================================================================

#============================================================================== 
Example 13
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMappingTools.py    MIT License 6 votes vote down vote up
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: 
Example 14
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_main.py    GNU General Public License v2.0 6 votes vote down vote up
def mtlBrowse (self):
        """
        Open browse dialog for metadata file selecting
        """
        filename = QtGui.QFileDialog.getOpenFileName(self, 'Open file', '', '*.txt')
        if filename:
            self.ui.satTabMTLPathLine.setText(filename)
        pass

    ##############################################################
    ################ END INTERFACE
    ##############################################################

    ##############################################################
    ################ CHECK INPUTS
    ############################################################## 
Example 15
Project: lidar   Author: giswqs   File: slicing.py    MIT License 6 votes vote down vote up
def writeRaster(arr, out_path, template):
    no_data = 0
    # First of all, gather some information from the template file
    data = gdal.Open(template)
    [cols, rows] = arr.shape
    trans = data.GetGeoTransform()
    proj = data.GetProjection()
    # nodatav = 0 #data.GetNoDataValue()
    # Create the file, using the information from the template file
    outdriver = gdal.GetDriverByName("GTiff")
    # http://www.gdal.org/gdal_8h.html
    # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,
    outdata   = outdriver.Create(str(out_path), rows, cols, 1, gdal.GDT_UInt32)
    # Write the array to the file, which is the original array in this example
    outdata.GetRasterBand(1).WriteArray(arr)
    # Set a no data value if required
    outdata.GetRasterBand(1).SetNoDataValue(no_data)
    # Georeference the image
    outdata.SetGeoTransform(trans)
    # Write projection information
    outdata.SetProjection(proj)
    return arr


# raster to vector 
Example 16
Project: lidar   Author: giswqs   File: filling.py    MIT License 6 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)

    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    # raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# extract sinks from dem 
Example 17
Project: eo_tools   Author: neel9102   File: reflectance_landsat_8_single_release.py    GNU General Public License v3.0 6 votes vote down vote up
def return_band(image_file_name, band_number):
    image = image_file_name
    dataset = gdal.Open(image,GA_ReadOnly)  
    if dataset is None:
        print "Could not open " + dataset
        sys.exit(1)
    geoTransform = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    rasterband = dataset.GetRasterBand(band_number)
    type(rasterband)
    ncol = dataset.RasterXSize
    nrow = dataset.RasterYSize
    band = rasterband.ReadAsArray(0,0,ncol,nrow)
    band = band.astype(numpy.uint16)
    return band,geoTransform,proj,ncol,nrow
    dataset = None
    band = None



# will return '/media/Arc/eo_archive_proc/VHR_SAT_IMAGE/SPOT6/20140704_SPOT/binned_SPOT6_20140704/B0.binned_SPOT6_20140704.tif'
# the function input defined in the beginining: out_put_dir, product just we have to change the product name..... 
Example 18
Project: eo_tools   Author: neel9102   File: reflectance_landsat_8_multiple_release.py    GNU General Public License v3.0 6 votes vote down vote up
def return_band(image_file_name, band_number):
    image = image_file_name
    dataset = gdal.Open(image,GA_ReadOnly)  
    if dataset is None:
        print "Could not open " + dataset
        sys.exit(1)
    geoTransform = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    rasterband = dataset.GetRasterBand(band_number)
    type(rasterband)
    ncol = dataset.RasterXSize
    nrow = dataset.RasterYSize
    band = rasterband.ReadAsArray(0,0,ncol,nrow)
    band = band.astype(numpy.uint16)
    return band,geoTransform,proj,ncol,nrow
    dataset = None
    band = None



# will return '/media/Arc/eo_archive_proc/VHR_SAT_IMAGE/SPOT6/20140704_SPOT/binned_SPOT6_20140704/B0.binned_SPOT6_20140704.tif'
# the function input defined in the beginining: out_put_dir, product just we have to change the product name..... 
Example 19
Project: eo_tools   Author: neel9102   File: reflectance_RE_multiple_release.py    GNU General Public License v3.0 6 votes vote down vote up
def return_band(image_file_name, band_number):
    image = image_file_name
    dataset = gdal.Open(image,GA_ReadOnly)  
    if dataset is None:
        print "Could not open " + dataset
        sys.exit(1)
    geoTransform = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    rasterband = dataset.GetRasterBand(band_number)
    type(rasterband)
    ncol = dataset.RasterXSize
    nrow = dataset.RasterYSize
    band = rasterband.ReadAsArray(0,0,ncol,nrow)
    band = band.astype(numpy.uint16)
    return band,geoTransform,proj,ncol,nrow
    dataset = None
    band = None


# will return '/media/Arc/eo_archive_proc/VHR_SAT_IMAGE/SPOT6/20140704_SPOT/binned_SPOT6_20140704/B0.binned_SPOT6_20140704.tif' 
Example 20
Project: eo_tools   Author: neel9102   File: mod09q1_lta_ndvi_computation_release.py    GNU General Public License v3.0 6 votes vote down vote up
def return_band(image_file_name, band_number):
    image = image_file_name
    dataset = gdal.Open(image,GA_ReadOnly)  
    if dataset is None:
        print "Could not open " + dataset
        sys.exit(1)
    geoTransform = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    rasterband = dataset.GetRasterBand(band_number)
    type(rasterband)
    ncol = dataset.RasterXSize
    nrow = dataset.RasterYSize
    band = rasterband.ReadAsArray(0,0,ncol,nrow)
    band = band.astype(numpy.uint16)
    return band,geoTransform,proj,ncol,nrow
    dataset = None
    band = None


# function to read the image bands 
Example 21
Project: eo_tools   Author: neel9102   File: mod09q1_lta_ndvi_computation_release.py    GNU General Public License v3.0 6 votes vote down vote up
def return_band(image_file_name, band_number):
    image = image_file_name
    dataset = gdal.Open(image,GA_ReadOnly)  
    if dataset is None:
        print "Could not open " + dataset
        sys.exit(1)
    geoTransform = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    rasterband = dataset.GetRasterBand(band_number)
    type(rasterband)
    ncol = dataset.RasterXSize
    nrow = dataset.RasterYSize
    band = rasterband.ReadAsArray(0,0,ncol,nrow)
    band = band.astype(numpy.float32)
    return band,geoTransform,proj,ncol,nrow
    dataset = None
    band = None 
Example 22
Project: eo_tools   Author: neel9102   File: reflectance_wv2_mul_release.py    GNU General Public License v3.0 6 votes vote down vote up
def return_band(image_file_name, band_number):
    image = image_file_name
    dataset = gdal.Open(image,GA_ReadOnly)  
    if dataset is None:
        print "Could not open " + dataset
        sys.exit(1)
    geoTransform = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    rasterband = dataset.GetRasterBand(band_number)
    type(rasterband)
    ncol = dataset.RasterXSize
    nrow = dataset.RasterYSize
    band = rasterband.ReadAsArray(0,0,ncol,nrow)
    band = band.astype(numpy.uint16)
    return band,geoTransform,proj,ncol,nrow
    dataset = None
    band = None

# dfunction to efine the output name 
Example 23
Project: eo_tools   Author: neel9102   File: MAGIC_ypm_and_statistics_release.py    GNU General Public License v3.0 6 votes vote down vote up
def extract_crop_type_info(in_shp_path):
    crop_list = []
    for inshape in find_shp_files(in_shp_path):
        ds = ogr.Open(inshape)
        lyr = ds.GetLayer(0)
        lyr.ResetReading()
        ft = lyr.GetNextFeature()
        while ft:
            field_name = ft.GetFieldAsString('field_name')
            field_name = field_name.replace('_', '')
            field_name = field_name.replace('.', '')
            crop_type = ft.GetFieldAsString('crop')
            x = [field_name, crop_type]
            if x not in crop_list:
                crop_list.append(x)
            ft = lyr.GetNextFeature()
        ds = None
    return crop_list 
Example 24
Project: unmixing   Author: arthur-e   File: utils.py    MIT License 6 votes vote down vote up
def as_array(path, band_axis=True):
    '''
    A convenience function for opening a raster as an array and accessing its
    spatial information; takes a single string argument. Arguments:
        path        The path of the raster file to open as an array
        band_axis   True to include a band axis, even for single-band rasters
    '''
    ds = gdal.Open(path)
    arr = ds.ReadAsArray()
    gt = ds.GetGeoTransform()
    wkt = ds.GetProjection()
    ds = None

    # Make sure that single-band rasters have a band axis
    if band_axis and len(arr.shape) < 3:
        shp = arr.shape
        arr = arr.reshape((1, shp[0], shp[1]))

    return (arr, gt, wkt) 
Example 25
Project: core   Author: lifemapper   File: create_mask.py    GNU General Public License v3.0 5 votes vote down vote up
def get_layer_dimensions(template_layer_filename):
    """
    @summary: Get layer information
    @param template_layer_filename: The template layer to get information for
    """
    ds = gdal.Open(template_layer_filename)

    num_cols = ds.RasterXSize
    num_rows = ds.RasterYSize
    
    minx, cell_size, _, maxy, _, _ = ds.GetGeoTransform()
    
    maxx = minx + (cell_size * num_cols)
    miny = maxy - (cell_size * num_rows)
    
    bbox = (minx, miny, maxx, maxy)

    srs = osr.SpatialReference()
    srs.ImportFromWkt(ds.GetProjectionRef())
    # TODO: Try to find a better way to fail
    try:
        epsg = int(srs.GetAttrValue('AUTHORITY', 1))
    except:
        epsg = DEFAULT_EPSG
    ds = None
    
    return bbox, cell_size, epsg

# ............................................................................. 
Example 26
Project: core   Author: lifemapper   File: raster_validator.py    GNU General Public License v3.0 5 votes vote down vote up
def validate_raster_file(raster_filename, raster_format=None):
    """Validates a raster file by seeing if it can be loaded by GDAL

    Args:
        raster_filename : A file path for a raster file to validate.
        raster_format : An optional LMFormat class to use for validation.  If
            not provided, the file extension will be used to determine how to
            validate the file.
    """
    msg = 'Valid'
    valid = False
    if os.path.exists(raster_filename):
        # If a raster format was not provided, try to get it from the file 
        #    extension
        if raster_format is None:
            _, ext = os.path.splitext(raster_filename)
            if ext == LMFormat.ASCII.ext:
                raster_format = LMFormat.ASCII
            elif ext == LMFormat.GTIFF.ext:
                raster_format = LMFormat.GTIFF
            else:
                msg = 'Extension {} did not map to a known raster format'.format(
                                                                           ext)

        if raster_format is not None:
            try:
                ds = gdal.Open(raster_filename)
                if ds is None:
                    msg = 'Could not open {}'.format(raster_filename)
                else:
                    lyr = ds.GetRasterBand(1)
                    valid = True
            except Exception, e:
                msg = str(e)
        else:
            msg = 'File does not exist' 
Example 27
Project: core   Author: lifemapper   File: layerset.py    GNU General Public License v3.0 5 votes vote down vote up
def _getRasterInfo(self, srcpath, getHisto=False):
        """
        @summary: Uses GDAL to retrieve the minimum and maximum values from a 
                  RASTER data source.  Note that for some types of data source 
                  (like ASCII grids), this process can be quite slow.
        @param srcpath: full path to the raster dataset
        @return: list of [min,max,nodata]
        """
        try:
            src = gdal.Open(srcpath, gdalconst.GA_ReadOnly)
        except Exception, e:
            print('Exception opening {} ({})'.format(srcpath, e))
            return (None, None, None, None) 
Example 28
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def _openWithGDAL(self, dlocation=None, bandnum=1):
        """
        @return: a GDAL dataset object
        """
        if dlocation is None:
            dlocation = self._dlocation
        try:
            dataset = gdal.Open(str(dlocation), gdalconst.GA_ReadOnly)
            band = dataset.GetRasterBand(1)
        except Exception, e:
            raise LMError(['Unable to open dataset or band {} with GDAL ({})'
                                .format(dlocation, str(e))]) 
Example 29
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def getSRS(self):
        if (self._dlocation is not None and os.path.exists(self._dlocation)):
            ds = gdal.Open(str(self._dlocation), gdalconst.GA_ReadOnly)
            wktSRS = ds.GetProjection()
            if wktSRS is not '':
                srs = osr.SpatialReference()
                srs.ImportFromWkt(wktSRS)
            else:
                srs = self.createSRSFromEPSG()
            ds = None
            return srs
        else:
            raise LMError(currargs='Input file %s does not exist' % self._dlocation)

# ............................................... 
Example 30
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def isValidDataset(self):
        """
        @summary: Checks to see if dataset is a valid raster
        @return: True if raster is a valid GDAL dataset; False if not
        """
        valid = True
        if (self._dlocation is not None and os.path.exists(self._dlocation)):
            try:
                self._dataset = gdal.Open(self._dlocation, gdalconst.GA_ReadOnly)
            except Exception, e:
                valid = False 
Example 31
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def verifyField(self, dlocation, ogrFormat, attrname):
        """
        @summary: Read OGR-accessible data and save the features and 
                     featureAttributes on the Vector object
        @param dlocation: Location of the data
        @param ogrFormat: OGR-supported data format code, available at
                             http://www.gdal.org/ogr/ogr_formats.html
        @return: boolean for success/failure 
        @raise LMError: on failure to read data.
        @note: populateStats calls this
        """
        if attrname is None:
            fieldOk = True
        else:
            fieldOk = False
            if dlocation is not None and os.path.exists(dlocation):
                ogr.RegisterAll()
                drv = ogr.GetDriverByName(ogrFormat)
                try:
                    ds = drv.Open(dlocation)
                except Exception, e:
                    raise LMError(['Invalid datasource' % dlocation, str(e)])
                
                lyrDef = ds.GetLayer(0).GetLayerDefn()
                # Make sure given field exists and is the correct type
                for i in range(lyrDef.GetFieldCount()):
                    fld = lyrDef.GetFieldDefn(i)
                    fldname = fld.GetNameRef()
                    if attrname == fldname:
                        fldtype = fld.GetType()
                        if fldtype in (ogr.OFTInteger, ogr.OFTReal, ogr.OFTBinary):
                            fieldOk = True
                        break 
Example 32
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def readWithOGR(self, dlocation, ogrFormat, featureLimit=None, doReadData=False):
        """
        @summary: Read OGR-accessible data and set the features and 
                     featureAttributes on the Vector object
        @param dlocation: Location of the data
        @param ogrFormat: OGR-supported data format code, available at
                             http://www.gdal.org/ogr/ogr_formats.html
        @return: boolean for success/failure 
        @raise LMError: on failure to read data.
        @note: populateStats calls this
        @todo: remove featureLimit, read subsetDLocation if there is a limit 
        """
        thisBBox = localIdIdx = geomIdx = feats = featAttrs = None
        if dlocation is not None and os.path.exists(dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(ogrFormat)
            try:
                ds = drv.Open(dlocation)
            except Exception, e:
                raise LMError(['Invalid datasource' % dlocation, str(e)])
                          
            self.clearFeatures() 
            try:
                slyr = ds.GetLayer(0)
            except Exception, e:
                raise LMError(currargs='#### Failed to GetLayer from %s' % dlocation,
                                  prevargs=e.args, doTrace=True)

            # Get bounding box 
Example 33
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def getSRS(self):
        if self._dlocation is not None and os.path.exists(self._dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(self._dataFormat)
            try:
                ds = drv.Open(self._dlocation)
            except Exception, e:
                raise LMError(['Invalid datasource' % self._dlocation, str(e)])

            vlyr = ds.GetLayer(0)
            srs = vlyr.GetSpatialRef()
            if srs is None:
                srs = self.createSRSFromEPSG()
            return srs 
Example 34
Project: core   Author: lifemapper   File: lmmap.py    GNU General Public License v3.0 5 votes vote down vote up
def _getRasterInfo(self, srcpath, getHisto=False):
      """
      @summary: Uses GDAL to retrieve the minimum and maximum values from a 
                RASTER data source.  Note that for some types of data source 
                (like ASCII grids), this process can be quite slow.
      @param srcpath: full path to the raster dataset
      @return: list of [min,max,nodata]
      """
      try:
         src = gdal.Open(srcpath, gdalconst.GA_ReadOnly)
      except Exception, e:
         print('Exception opening %s (%s)' % (srcpath,str(e)) )
         return (None, None, None, None) 
Example 35
Project: core   Author: lifemapper   File: create_blank_mask.py    GNU General Public License v3.0 5 votes vote down vote up
def createMaskRaster(inRasterFn, outRasterFn):
   """
   @summary: Generate a mask layer based on the input raster layer
   """
   ds = gdal.Open(inRasterFn)
   band = ds.GetRasterBand(1)
   
   cols = ds.RasterXSize
   rows = ds.RasterYSize
   
   data = band.ReadAsArray(0, 0, cols, rows)
   
   newData = np.ones(data.shape, dtype=int)
   
   drv = ds.GetDriver()
   outDs = drv.Create(outRasterFn, cols, rows, 1, gdalconst.GDT_Int32)
   outBand = outDs.GetRasterBand(1)
   
   # Write the data
   outBand.WriteArray(newData, 0, 0)
   outBand.FlushCache()
   outBand.SetNoDataValue(DEFAULT_NODATA)
   
   # Georeference the image and set the projection
   outDs.SetGeoTransform(ds.GetGeoTransform())
   outDs.SetProjection(ds.GetProjection())
   del newData

    
# ................................................................ 
Example 36
Project: core   Author: lifemapper   File: geotools.py    GNU General Public License v3.0 5 votes vote down vote up
def getArray(self, bandnum, doFlip=False, doShift=False):
      """
      @summary: Read the dataset into numpy array  
      @param bandnum: The band number to read.
      @param doFlip: True if data begins at the southern edge of the region
      @param doShift: True if the leftmost edge of the data should be shifted 
             to the center (and right half shifted around to the beginning) 
      """
      if 'numpy' in dir():
         inds = gdal.Open(self.dlocation, gdalconst.GA_ReadOnly)
         inband = inds.GetRasterBand(bandnum)
         arrtype = self._getNumpyType(self.gdalBandType)
         outArr = numpy.empty([self.ysize, self.xsize], dtype=arrtype)
         
         for row in range(self.ysize):
            scanline = inband.ReadAsArray(0, row, self.xsize, 1, self.xsize, 1)
            
            if doShift:
               scanline = self._cycleRow(scanline, arrtype, 0, self.xsize/2, 
                                         self.xsize)
            if doFlip:
               newrow = self.ysize-row-1
            else:
               newrow = row
   
            outArr[newrow] = scanline
   
         inds = None   
         return outArr
      else:
         raise LMError('numpy missing - unable to getArray')
      
# ............................................. 
Example 37
Project: core   Author: lifemapper   File: geotools.py    GNU General Public License v3.0 5 votes vote down vote up
def copyDataset(self, bandnum, outfname, format='GTiff', kwargs={}):
      """
      @summary: Copy the dataset into a new file.  
      @param bandnum: The band number to read.
      @param outfname: Filename to write this dataset to.
      @param format: GDAL-writeable raster format to use for new dataset. 
                     http://www.gdal.org/formats_list.html
      @param doFlip: True if data begins at the southern edge of the region
      @param doShift: True if the leftmost edge of the data should be shifted 
             to the center (and right half shifted around to the beginning) 
      @param nodata: Value used to indicate nodata in the new file.
      @param srs: Spatial reference system to use for the data. This is only 
                  necessary if the dataset does not have an SRS present.  This
                  will NOT project the dataset into a different projection.
      """
      driver = gdal.GetDriverByName(format)
      metadata = driver.GetMetadata()
      if not (metadata.has_key(gdal.DCAP_CREATECOPY) 
                and metadata[gdal.DCAP_CREATECOPY] == 'YES'):
         raise LMError(currargs='Driver %s does not support CreateCopy() method.' 
                       % format)
      inds = gdal.Open( self.dlocation )
      try:
         outds = driver.CreateCopy(outfname, inds, 0, **kwargs)
      except Exception, e:
         raise LMError(currargs='Creation failed for %s from band %d of %s (%s)'
                                % (outfname, bandnum, self.dlocation, str(e))) 
Example 38
Project: core   Author: lifemapper   File: geotools.py    GNU General Public License v3.0 5 votes vote down vote up
def loadBand(self, bandnum=1):
      """
      @summary: Open the dataset and save the band to a member attribute for 
                further examination.  
      @param bandnum: The band number to read.
      """
      if self._band is None:
         if self._dataset is None:
            self.openDataSet(None)
         self._band = self._dataset.GetRasterBand(bandnum)
    
# ............................................... 
Example 39
Project: core   Author: lifemapper   File: geotools.py    GNU General Public License v3.0 5 votes vote down vote up
def rasterSize(datasrc):
      '''
      @summary: Return [width, height] in pixels
      '''
      datasrc = gdal.Open(str(datasrc))
      return [datasrc.RasterXSize, datasrc.RasterYSize]      

# ............................................... 
Example 40
Project: core   Author: lifemapper   File: layer_encoder.py    GNU General Public License v3.0 5 votes vote down vote up
def _read_raster_layer(self, raster_filename):
        """Reads a raster layer for processing

        Args:
            raster_filename: The file path for the raster layer.

        Returns:
            A tuple containing a window function for returning a portion of the
            numpy array generated by the layer and the NODATA value to use with
            this layer.
        """
        ds = gdal.Open(raster_filename)
        band = ds.GetRasterBand(1)
        layer_array = band.ReadAsArray()
        nodata = band.GetNoDataValue()
        
        num_y, num_x = layer_array.shape
        min_x, x_res, _, max_y, _, y_res = ds.GetGeoTransform()
        max_x = min_x + (num_x * x_res)
        min_y = max_y + (y_res * num_y)
        layer_bbox = (min_x, min_y, max_x, max_y)
        window_func = self._get_window_function(
            layer_array, layer_bbox, self.shapegrid_resolution,
            num_cell_sides=self.shapegrid_sides)
        
        return (window_func, nodata)
    
    # ............................... 
Example 41
Project: core   Author: lifemapper   File: layer_encoder.py    GNU General Public License v3.0 5 votes vote down vote up
def _read_shapegrid(self, shapegrid_filename):
        """Read the shapegrid

        Args:
            shapegrid_filename: The file location of the shapegrid
        """
        shapegrid_dataset = ogr.Open(shapegrid_filename)
        self.shapegrid_layer = shapegrid_dataset.GetLayer()
        tmp = self.shapegrid_layer.GetExtent()
        self.shapegrid_bbox = (tmp[0], tmp[2], tmp[1], tmp[3])
        
        self.num_cells = self.shapegrid_layer.GetFeatureCount()

        feature_0 = self.shapegrid_layer.GetFeature(0)
        geom = feature_0.GetGeometryRef()
        
        # Get resolution and number of sides
        geom_wkt = geom.ExportToWkt()
        boundary_points = geom_wkt.split(',')
        if len(boundary_points) == 5:
            # Square
            envelope = geom.GetEnvelope()
            self.shapegrid_resolution = (envelope[1] - envelope[0],
                                         envelope[3] - envelope[2])
        else:
            # Hexagon
            center = geom.Centroid()
            x_cent = center.GetX()
            y_cent = center.GetY()
            x1, y1 = boundary_points[1].split(' ')
            self.shapegrid_resolution = np.sqrt(
                (x_cent - x1)**2 + (y_cent - y1)**2)
        self.shapegrid_sides = len(boundary_points) - 1
        #self.shapegrid_layer.ResetReading()
        self.shapegrid_layer = None

    # ............................... 
Example 42
Project: core   Author: lifemapper   File: layer_encoder.py    GNU General Public License v3.0 5 votes vote down vote up
def get_geojson(self):
        """Formats the encoded matrix as GeoJSON
        """
        ret = {
            'type' : 'FeatureCollection'
        }
        features = []
        
        column_headers = self.encoded_matrix.getColumnHeaders()
        
        column_enum = [(j, str(k)) for j, k in enumerate(column_headers)]
        
        shapegrid_dataset = ogr.Open(self.shapegrid_filename)
        shapegrid_layer = shapegrid_dataset.GetLayer()
        
        i = 0
        feat = shapegrid_layer.GetNextFeature()
        while feat is not None:
            ft_json = json.loads(feat.ExportToJson())
            # right_hand_rule(ft_json['geometry']['coordinates'])
            # TODO(CJ): Remove this if updated library adds first id correctly
            ft_json['id'] = feat.GetFID()
            ft_json['properties'] = dict(
                [(k, self.encoded_matrix.data[i, j].item()
                  ) for j, k in column_enum])
            features.append(ft_json)
            i+=1
            feat = shapegrid_layer.GetNextFeature()

        ret['features'] = features
        shapegrid_dataset = shapegrid_layer = None
        return ret 
Example 43
Project: spacesense   Author: spacesense-ai   File: datasets.py    GNU Lesser General Public License v3.0 5 votes vote down vote up
def __init__(self, hdf_data_path):
        self.hdf_data_path = hdf_data_path
        df = gdal.Open(self.hdf_data_path, gdal.GA_ReadOnly)
        self.meta_data = df.GetMetadata()
        self.sdf = df.GetSubDatasets()
        self.band_details = [self.sdf[i][1] for i in range(len(self.sdf))]
        self.num_bands = len(self.sdf)
        self.img_shp = gdal.Open(self.sdf[0][0]).ReadAsArray().shape
        self.AOI = None
        self.data = None 
Example 44
Project: spacesense   Author: spacesense-ai   File: datasets.py    GNU Lesser General Public License v3.0 5 votes vote down vote up
def get_data(self):
        self.data = np.zeros((self.img_shp[0], self.img_shp[1], self.num_bands))
        for i in range(self.num_bands):
            self.data[:, :, i] = gdal.Open(self.sdf[i][0]).ReadAsArray() 
Example 45
Project: fetchLandsatSentinelFromGoogleCloud   Author: vascobnunes   File: fels.py    MIT License 5 votes vote down vote up
def check_full_tile(image):
    gdalData = gdal.Open(image)
    if gdalData is None:
        sys.exit("ERROR: can't open raster")

    # get width and heights of the raster
    xsize = gdalData.RasterXSize
    ysize = gdalData.RasterYSize

    # process the raster
    band_i = gdalData.GetRasterBand(1)
    raster = band_i.ReadAsArray()

    # create dictionary for unique values count
    count = {}

    # count unique values for the given band
    for col in range(xsize):
        for row in range(ysize):
            cell_value = raster[row, col]

            # check if cell_value is NaN
            if cell_value == 0:
                # add cell_value to dictionary
                if cell_value in count:
                    count[cell_value] += 1
                else:
                    count[cell_value] = 1
                break
    for key in sorted(count.keys()):
        if count[key] is not None:
            return "Partial" 
Example 46
Project: radiometric_normalization   Author: planetlabs   File: gimage.py    Apache License 2.0 5 votes vote down vote up
def load(filename, nodata=None, last_band_alpha=False):
    logging.info('GImage: Loading {} as GImage'.format(filename))
    gdal_ds = gdal.Open(filename)
    if gdal_ds is None:
        raise Exception('Unable to open file "{}" with gdal.Open()'.format(
            filename))

    alpha, band_count = read_alpha_and_band_count(gdal_ds, last_band_alpha)
    bands = _read_all_bands(gdal_ds, band_count)
    metadata = read_metadata(gdal_ds)

    if nodata:
        alpha = alpha * _nodata_to_mask(bands, nodata)
    return GImage(bands, alpha, metadata) 
Example 47
Project: radiometric_normalization   Author: planetlabs   File: normalize_wrapper.py    Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example 48
Project: radiometric_normalization   Author: planetlabs   File: display_wrapper.py    Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example 49
Project: radiometric_normalization   Author: planetlabs   File: transformation_wrapper.py    Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example 50
Project: radiometric_normalization   Author: planetlabs   File: pif_wrapper.py    Apache License 2.0 5 votes vote down vote up
def _open_image_and_get_info(path, last_band_alpha):
    gdal_ds = gdal.Open(path)
    alpha_band, band_count = gimage.read_alpha_and_band_count(
        gdal_ds, last_band_alpha=last_band_alpha)
    return gdal_ds, alpha_band, band_count 
Example 51
Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    Apache License 2.0 5 votes vote down vote up
def test__read_all_bands(self):
        gdal_ds = gdal.Open(self.test_photometric_alpha_image)
        bands = gimage._read_all_bands(gdal_ds, 3)
        numpy.testing.assert_array_equal(bands[0], self.band)
        numpy.testing.assert_array_equal(bands[1], self.band)
        numpy.testing.assert_array_equal(bands[2], self.band) 
Example 52
Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    Apache License 2.0 5 votes vote down vote up
def test_read_single_band(self):
        gdal_ds = gdal.Open(self.test_photometric_alpha_image)
        band = gimage.read_single_band(gdal_ds, 3)
        numpy.testing.assert_array_equal(band, self.band) 
Example 53
Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    Apache License 2.0 5 votes vote down vote up
def test_read_alpha_and_band_count(self):
        gdal_ds = gdal.Open(self.test_photometric_alpha_image)
        alpha, band_count = gimage.read_alpha_and_band_count(gdal_ds)

        self.assertEqual(band_count, 3)
        numpy.testing.assert_array_equal(alpha, self.mask) 
Example 54
Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    Apache License 2.0 5 votes vote down vote up
def test__save_to_ds(self):
        output_file = 'test_save_to_ds.tif'

        test_band = numpy.array([[0, 1], [2, 3]], dtype=numpy.uint16)
        test_gimage = gimage.GImage([test_band], self.mask, self.metadata)
        output_ds = gdal.GetDriverByName('GTiff').Create(
            output_file, 2, 2, 2, gdal.GDT_UInt16,
            options=['ALPHA=YES'])
        gimage._save_to_ds(test_gimage, output_ds, nodata=3)

        # Required for gdal to write to file
        output_ds = None

        test_ds = gdal.Open(output_file)

        saved_number_of_bands = test_ds.RasterCount
        self.assertEquals(saved_number_of_bands, 2)

        saved_band = test_ds.GetRasterBand(1).ReadAsArray()
        numpy.testing.assert_array_equal(saved_band, self.band)

        saved_nodata = test_ds.GetRasterBand(1).GetNoDataValue()
        self.assertEqual(saved_nodata, 3)

        saved_alpha = test_ds.GetRasterBand(2).ReadAsArray()
        numpy.testing.assert_array_equal(saved_alpha, self.mask * 255)

        os.unlink(output_file) 
Example 55
Project: vmap   Author: dshean   File: vmap.py    MIT License 5 votes vote down vote up
def get_vel(fn, fill=True):
    ds = gdal.Open(fn)
    if fill:
        import dem_downsample_fill
        ds = dem_downsample_fill.gdalfill_ds(ds)
    u_b = ds.GetRasterBand(1)
    v_b = ds.GetRasterBand(2)
    u = iolib.b_getma(u_b)
    v = iolib.b_getma(v_b)
    m = np.ma.sqrt(u*u + v*v)
    return u, v, m 
Example 56
Project: sentinelloader   Author: flaviostutz   File: sentinel2loader.py    MIT License 5 votes vote down vote up
def cropRegion(self, geoPolygon, sourceGeoTiffs):
        """Returns an image file with contents from a bunch of GeoTiff files cropped to the specified geoPolygon.
           Pay attention to the fact that a new file is created at each request and you should delete it after using it"""
        logger.debug("Cropping polygon from %d files" % (len(sourceGeoTiffs)))
        desiredRegion = Polygon(geoPolygon)

#         #show tile images
#         for fn in tilesData:
#             ds = gdal.Open(fn).ReadAsArray()
#             plt.figure(figsize=(10,10))
#             plt.imshow(ds[0])

        source_tiles = ' '.join(sourceGeoTiffs)
        tmp_file = "%s/tmp/%s.tiff" % (self.dataPath, uuid.uuid4().hex)
        if not os.path.exists(os.path.dirname(tmp_file)):
            os.makedirs(os.path.dirname(tmp_file))
        
        #define output bounds in destination srs reference
        bounds = desiredRegion.bounds
        s1 = convertWGS84To3857(bounds[0], bounds[1])
        s2 = convertWGS84To3857(bounds[2], bounds[3])

        logger.debug('Combining tiles into a single image. sources=%s tmpfile=%s' % (source_tiles, tmp_file))
        os.system("gdalwarp -multi -srcnodata 0 -t_srs EPSG:3857 -te %s %s %s %s %s %s" % (s1[0],s1[1],s2[0],s2[1],source_tiles,tmp_file))
        
        return tmp_file 
Example 57
Project: sentinelloader   Author: flaviostutz   File: sentinel2loader.py    MIT License 5 votes vote down vote up
def _getBandDataFloat(self, geoPolygon, bandName, resolution, dateReference):
        bandFile = self.getRegionBand(geoPolygon, bandName, resolution, dateReference)
        
        gdalBand = gdal.Open(bandFile)
        geoTransform = gdalBand.GetGeoTransform()
        projection = gdalBand.GetProjection()
        
        data = gdalBand.ReadAsArray().astype(np.float)
        os.remove(bandFile)
        return data, geoTransform, projection 
Example 58
Project: usv_sim_lsa   Author: disaster-robotics-proalertas   File: windLBM.py    Apache License 2.0 5 votes vote down vote up
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) 
Example 59
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    MIT License 5 votes vote down vote up
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 
Example 60
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    MIT License 5 votes vote down vote up
def array2raster(rasterfn,newRasterfn,array,driver_name = "ENVI", noDataValue = -9999):
    """Takes an array and writes to a GDAL compatible raster. It needs another raster to map the dimensions.

    Args:
        FileName (str): The filename (with path and extension) of a raster that has the same dimensions as the raster to be written.
        newRasterfn (str): The filename (with path and extension) of the new raster.
        array (np.array): The array to be written
        driver_name (str): The type of raster to write. Default is ENVI since that is the LSDTOpoTools format
        noDataValue (float): The no data value

    Return:
        np.array: A numpy array with the data from the raster.

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName(driver_name)
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outRaster.GetRasterBand(1).SetNoDataValue( noDataValue )
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
#============================================================================== 
Example 61
Project: LSDMappingTools   Author: LSDtopotools   File: rotated_mapping_tools.py    MIT License 5 votes vote down vote up
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) 
Example 62
Project: LSDMappingTools   Author: LSDtopotools   File: rotated_mapping_tools.py    MIT License 5 votes vote down vote up
def ReadShapeFile(ShapeFile):

    """
    Open shapefile and create Polygon dictionary
    returns dictionary of shapely Polygons

    MDH

    """

    #open shapefile and read shapes
    Shapes = fiona.open(ShapeFile)

    # get the input coordinate system
    Input_CRS = Proj(Shapes.crs)

    # Create a dictionary of shapely polygons
    PolygonDict = {}

    # loop through shapes and add to dictionary
    for Feature in Shapes:
        if Feature['geometry']['type'] == 'MultiPolygon':
            Shape = MultiPolygon(shape(Feature['geometry']))
            Value = float(Feature['properties']['ID'])
        elif Feature['geometry']['type'] == 'Polygon':
            Shape = Polygon(shape(Feature['geometry']))
            Value = float(Feature['properties']['ID'])

        #check for multipolygons
        if Value in PolygonDict:
            Polygons = [Shape, PolygonDict[Value]]
            Shape = cascaded_union(Polygons)
        #update dictionary
        PolygonDict[Value] = Shape

    return PolygonDict, Input_CRS 
Example 63
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_common_lib.py    GNU General Public License v2.0 5 votes vote down vote up
def Landsat8_simple_temperature (mode, raster_path, channel_number, metadata_path, output_path=None, base_raster=None):
    # mode 0 - generate tiff
    # mode 1 - return np array

    """
    Return brigtness temperature for Landsat's 8 B10 or B11.
    :param mode: int; mode 0 - generate tiff; mode 1 - return np array
    :param raster_path: string; input raster path
    :param channel_number: int; number of channel to be processed (10 or 11)
    :param metadata_path: string; path to metadata file
    :param output_path: string; path for output brightness temperature file
    :param base_raster: base geotiff layer (needed if mode = 0)
    :return: string; path for output radiance file
    """
    try:
        int(channel_number)
    except:
        raise NameError('Channel number is not correct')

    if (channel_number < 10) or (channel_number > 11):
        raise NameError('Channel number is not correct')

    landsat_radiance_band_array = Landsat8_DN_to_radiance(1,raster_path,channel_number,metadata_path)
    K1_constant = float(read_metadata_parameter(metadata_path,'K1_CONSTANT_BAND_' + str(channel_number)))
    K2_constant = float(read_metadata_parameter(metadata_path,'K2_CONSTANT_BAND_' + str(channel_number)))

    landsat_temperature_array = (K2_constant / np.log((K1_constant/landsat_radiance_band_array)+1)) - 273.15
    del landsat_radiance_band_array
    if mode == 0:
        base_raster = gdal.Open(re.sub("^\s+|\n|\r|\s+$", '', raster_path))
        cols = base_raster.RasterXSize
        rows = base_raster.RasterYSize
        cell_type = gdal.GDT_Float32
        driver_name = 'GTiff'
        projection = base_raster.GetProjection()
        transform = base_raster.GetGeoTransform()
        save_nparray_as_raster(landsat_temperature_array,driver_name,cell_type,cols,rows,projection,transform,output_path)
        del landsat_temperature_array
    else:
        return landsat_temperature_array 
Example 64
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_main.py    GNU General Public License v2.0 5 votes vote down vote up
def openModisSettings(self):
        """
        Open settings for MOD09 Retrieving
        """
        self.ModisSettings = l8_lst_swa_settings.L8_lst_swaSettingsDlg()
        self.ModisSettings.show() 
Example 65
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_main.py    GNU General Public License v2.0 5 votes vote down vote up
def outputBrowse(self):
        """
        Open browse dialog for output file selecting
        """
        filename = QtGui.QFileDialog.getSaveFileName(self, 'Save file', '', '*.tif')
        if filename:
            self.ui.outputLSTLine.setText(filename)
        pass 
Example 66
Project: landuse   Author: mchristoffersen   File: gdal_merge.py    MIT License 5 votes vote down vote up
def init_from_name(self, filename):
        """
        Initialize file_info from filename

        filename -- Name of file to read.

        Returns 1 on success or 0 if the file can't be opened.
        """
        fh = gdal.Open( filename )
        if fh is None:
            return 0

        self.filename = filename
        self.bands = fh.RasterCount
        self.xsize = fh.RasterXSize
        self.ysize = fh.RasterYSize
        self.band_type = fh.GetRasterBand(1).DataType
        self.projection = fh.GetProjection()
        self.geotransform = fh.GetGeoTransform()
        self.ulx = self.geotransform[0]
        self.uly = self.geotransform[3]
        self.lrx = self.ulx + self.geotransform[1] * self.xsize
        self.lry = self.uly + self.geotransform[5] * self.ysize

        ct = fh.GetRasterBand(1).GetRasterColorTable()
        if ct is not None:
            self.ct = ct.Clone()
        else:
            self.ct = None

        return 1 
Example 67
Project: eo_tools   Author: neel9102   File: MAGIC_ypm_and_statistics_release.py    GNU General Public License v3.0 5 votes vote down vote up
def return_raster(in_file, band_number):
    dataset = gdal.Open(in_file)
    band = dataset.GetRasterBand(band_number)
    ncol = dataset.RasterXSize # xsize
    nrow = dataset.RasterYSize # ysize
    block_sizes = band.GetBlockSize()
    x_block = block_sizes[0]
    y_block = block_sizes[1]
    return dataset, band, ncol, nrow, x_block, y_block 
Example 68
Project: eo_tools   Author: neel9102   File: MAGIC_ypm_and_statistics_release.py    GNU General Public License v3.0 5 votes vote down vote up
def get_raster_array(in_raster, band_number):
    dataset = gdal.Open(in_raster)
    band = dataset.GetRasterBand(band_number)
    ncol = dataset.RasterXSize
    nrow = dataset.RasterYSize
    noDataValue = band.GetNoDataValue()
    numarray = band.ReadAsArray(0,0,ncol,nrow)
    numarray = numpy.ma.masked_values (numarray, noDataValue)
    print "reading bands as array....."
    dataList = []
    for x in xrange(numarray.shape[0]):
      for y in xrange(numarray.shape[1]):
        dataList.append(numarray[x][y])
    dataList = remove_values_from_list(dataList, noDataValue)
    return dataList 
Example 69
Project: eo_tools   Author: neel9102   File: MAGIC_clip_raster_by_poly_release.py    GNU General Public License v3.0 5 votes vote down vote up
def coord2pixelOffset(rasterfn,x,y):
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    xOffset = int((x - originX)/pixelWidth)
    yOffset = int((y - originY)/pixelHeight)
    return xOffset,yOffset 
Example 70
Project: unmixing   Author: arthur-e   File: utils.py    MIT License 5 votes vote down vote up
def as_raster(path):
    '''
    A convenience function for opening a raster and accessing its spatial
    information; takes a single string argument. Arguments:
        path    The path of the raster file to open as a gdal.Dataset
    '''
    ds = gdal.Open(path)
    gt = ds.GetGeoTransform()
    wkt = ds.GetProjection()
    return (ds, gt, wkt) 
Example 71
Project: unmixing   Author: arthur-e   File: utils.py    MIT License 5 votes vote down vote up
def array_to_raster(a, gt, wkt, xoff=None, yoff=None, dtype=None):
    '''
    Creates a raster from a given array, with optional x- and y-offsets
    if the array was clipped. Arguments:
        a       A NumPy array
        gt      A GDAL GeoTransform tuple
        wkt     Well-Known Text projection
        xoff    The offset in the x-direction; should be provided when clipped
        yoff    The offset in the y-direction; should be provided when clipped
        dtype   The data type to coerce on the array
    '''
    if dtype is not None:
        a = a.astype(dtype)

    try:
        rast = gdal_array.OpenNumPyArray(a)

    except AttributeError:
        # For backwards compatibility with older version of GDAL
        rast = gdal.Open(gdalnumeric.GetArrayFilename(a))

    except:
        rast = gdal_array.OpenArray(a)

    rast.SetGeoTransform(gt)
    rast.SetProjection(wkt)
    if xoff is not None and yoff is not None:
        # Bit of a hack; essentially, re-create the raster but with the
        #   correct X and Y offsets (don't know how to do this without the
        #   use of CopyDatasetInfo())
        return array_to_raster_clone(a, rast, xoff, yoff)

    return rast 
Example 72
Project: unmixing   Author: arthur-e   File: utils.py    MIT License 5 votes vote down vote up
def array_to_raster_clone(a, proto, xoff=None, yoff=None):
    '''
    Creates a raster from a given array and a prototype raster dataset, with
    optional x- and y-offsets if the array was clipped. Arguments:
        a       A NumPy array
        proto   A prototype dataset
        xoff    The offset in the x-direction; should be provided when clipped
        yoff    The offset in the y-direction; should be provided when clipped
    '''
    try:
        rast = gdal_array.OpenNumPyArray(a)

    except AttributeError:
        # For backwards compatibility with older version of GDAL
        rast = gdal.Open(gdalnumeric.GetArrayFilename(a))

    except:
        rast = gdal_array.OpenArray(a)

    kwargs = dict()
    if xoff is not None and yoff is not None:
        kwargs = dict(xoff=xoff, yoff=yoff)

    # Copy the projection info and metadata from a prototype dataset
    if type(proto) == str:
        proto = gdal.Open(proto)

    gdalnumeric.CopyDatasetInfo(proto, rast, **kwargs)
    return rast 
Example 73
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 4 votes vote down vote up
def zipShapefiles(self, baseName=None):
        """
        @summary: Returns a wrapper around a tar gzip file stream
        @param baseName: (optional) If provided, this will be the prefix for the 
                                  names of the shape file's files in the zip file.
        """
        fnames = self.getShapefiles()
        tgStream = StringIO()        
        zipf = zipfile.ZipFile(tgStream, mode="w", 
                                      compression=zipfile.ZIP_DEFLATED, allowZip64=True)
        if baseName is None:
            baseName = os.path.splitext(os.path.split(fnames[0])[1])[0]
        
        for fname in fnames:
            ext = os.path.splitext(fname)[1]
            zipf.write(fname, "%s%s" % (baseName, ext))
        zipf.close()
        
        tgStream.seek(0)
        ret = ''.join(tgStream.readlines())
        tgStream.close()
        return ret

# # ...............................................
#     def populateStats(self, dlocation, dataFormat):
#         """
#         @summary: Sets self.ogrType and self._bbox by opening dataset
#         """
#         msgs = []
#         if dlocation is not None:
#             if not os.path.exists(dlocation):
#                 msgs.append('Vector file does not exist: {}'.format(dlocation))
#             else:
#                 try:
#                     ds = ogr.Open(dlocation)
#                     ds = None
#                 except Exception, e:
#                     print('Unable to open vector dlocation {}: {}'
#                             .format(dlocation, str(e)) )
#                 try:
#                     self.readData(dlocation=dlocation, dataFormat=dataFormat)              
#                 except LMError, e:
#                     raise
#                 except Exception, e:
#                     raise LMError(currargs='Unable to read vector data {}: {}'
#                                        .format(dlocation, str(e)))
                         
# ............................................... 
Example 74
Project: core   Author: lifemapper   File: create_mask_tiff.py    GNU General Public License v3.0 4 votes vote down vote up
def createMaskRaster(inRasterFn, pointsFn, outRasterFn):
   pts = getXYsForShapefile(pointsFn)
   ds = gdal.Open(inRasterFn)
   band = ds.GetRasterBand(1)
   
   cols = ds.RasterXSize
   rows = ds.RasterYSize
   
   transform = ds.GetGeoTransform()
   xOrigin = transform[0]
   yOrigin = transform[3]
   pixelWidth = transform[1]
   pixelHeight = -transform[5]
   
   data = band.ReadAsArray(0, 0, cols, rows)
   
   vals = set([])
   for x,y in pts:
      col = int((x - xOrigin) / pixelWidth)
      row = int((yOrigin - y) / pixelHeight)
      try:
         vals.add(data[row][col])
      except:
         pass
   
   listVals = list(vals)
   if len(listVals) == 0:
      raise Exception, 'No intersection between points and raster'
   newData = DEFAULT_NODATA * np.ones(data.shape, dtype=int)
   
   for i in range(data.shape[0]):
      for j in range(data.shape[1]):
         if data[i,j] in listVals:
            newData[i,j] = 1
   
   drv = ds.GetDriver()
   outDs = drv.Create(outRasterFn, cols, rows, 1, gdalconst.GDT_Int32)
   outBand = outDs.GetRasterBand(1)
   
   # Write the data
   outBand.WriteArray(newData, 0, 0)
   outBand.FlushCache()
   outBand.SetNoDataValue(DEFAULT_NODATA)
   
   # Georeference the image and set the projection
   outDs.SetGeoTransform(ds.GetGeoTransform())
   outDs.SetProjection(ds.GetProjection())
   del newData

    
# ................................................................ 
Example 75
Project: core   Author: lifemapper   File: geotools.py    GNU General Public License v3.0 4 votes vote down vote up
def writeBand(self, bandnum, outfname, format='GTiff', doFlip=False, 
                doShift=False, nodata=None, srs=None):
      """
      @summary: Write the dataset into a new file, line by line.  
      @param bandnum: The band number to read.
      @param outfname: Filename to write this dataset to.
      @param format: GDAL-writeable raster format to use for new dataset. 
                     http://www.gdal.org/formats_list.html
      @param doFlip: True if data begins at the southern edge of the region
      @param doShift: True if the leftmost edge of the data should be shifted 
             to the center (and right half shifted around to the beginning) 
      @param nodata: Value used to indicate nodata in the new file.
      @param srs: Spatial reference system to use for the data. This is only 
                  necessary if the dataset does not have an SRS present.  This
                  will NOT project the dataset into a different projection.
      """
      driver = gdal.GetDriverByName(format)
      metadata = driver.GetMetadata()
      if not (metadata.has_key(gdal.DCAP_CREATE) 
              and metadata[gdal.DCAP_CREATE] == 'YES'):
         raise LMError(currargs='Driver %s does not support Create() method.' 
                                % format)
         
      outds = driver.Create(outfname, self.xsize, self.ysize, 1, self.gdalBandType)
      if outds is None:
         raise LMError('Creation failed for %s from band %d of %s'
                        % (outfname, bandnum, self.dlocation))
      
      outds.SetGeoTransform(self.geoTransform)
      
      inds = gdal.Open(self.dlocation, gdalconst.GA_ReadOnly)
      inband = inds.GetRasterBand(bandnum)
      outband = outds.GetRasterBand(1)
      
      if nodata is None:
         nodata = inband.GetNoDataValue()
      if nodata is not None:
         outband.SetNoDataValue(nodata)
      if srs is None:
         srs = self.srs
      outds.SetProjection(srs)
      
      for row in range(self.ysize):
         scanline = inband.ReadAsArray(0, row, self.xsize, 1, self.xsize, 1)
         
         if doShift:
            arrType = self._getNumpyType(self.gdalBandType)
            scanline = self._cycleRow(scanline, arrType, 0, self.xsize/2, 
                                      self.xsize)

         if doFlip:
            outband.WriteArray(scanline, 0, self.ysize-row-1)
         else:
            outband.WriteArray(scanline, 0, row)
      # Close new dataset to flush to disk
      outds = None
      inds = None
      
# ............................................... 
Example 76
Project: spacesense   Author: spacesense-ai   File: datasets.py    GNU Lesser General Public License v3.0 4 votes vote down vote up
def get_data(self, AOI='all', resample=False, resize_raster_source='B02',
                 interpolation=Resampling.cubic_spline, save_as_npy=False):
        """
        :return: (a dictionary of 13 bands),(single array of dimension (13,x_img,y_img) in the 'order of bands')
        """

        if AOI == 'all':
            self.band_files_new = self.band_files
            print('no AOI selected')
        else:
            # check for existing clipped raster files in tiff for given aoi
            aoi_name = AOI.split('/')[-1].split('.')[0]
            self.band_files_new = np.array(sorted(glob(self.save_folder + '/' + aoi_name + '*.tiff')))

        if len(self.band_files_new) == 0:
            if AOI != 'all':
                start = time.time()
                for source_raster in self.band_files:
                    clip_to_aoi(source_raster, AOI, self.save_folder)

                print('clipped to AOI in: ', time.time() - start, ' seconds')
                self.band_files_new = np.array(sorted(glob(self.save_folder + '/' + aoi_name + '*.tiff')))

        data_dictionary = {}
        print("loading data...")
        for band in self.band_files_new:
            if self.level =='2A':
                name = '_'.join(band.split('_')[-3:-1])
                key = band.split('_')[-2]
            else:
                name = band.split('/')[-1].split('.')[0]
                key = band.split('_')[-1].split('.')[0]
            print('loading ',name)
            arf = gdal.Open(band)
            data = arf.ReadAsArray()
            if len(data.shape)==3:
                data = data[0]

            if save_as_npy:
                self.save_as_npy(data, save_path=self.save_folder + '/' + name)

            data_dictionary[key] = data
        self.data_dictionary = data_dictionary
        print("all 13 bands loaded")
        if resample:
            print("resampling data started...")
            start = time.time()
            self.data_resampled = self.sentinel_2_remap(ref_raster=resize_raster_source, interpolation=interpolation)
            print('resampled in: ', time.time() - start, ' seconds')
        else:
            print('resampling not selected ')
            self.data_resampled =[]
        return self.data_dictionary, self.data_resampled 
Example 77
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    MIT License 4 votes vote down vote up
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

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


#============================================================================== 
Example 78
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_common_lib.py    GNU General Public License v2.0 4 votes vote down vote up
def Landsat8_DN_to_radiance (mode, raster_path, channel_number, metadata_path, output_path = None):
    """
    Convert raw Landsat8 scene from DN to radiance
    :param mode: int; mode 0 - generate tiff; mode 1 - return np array
    :param raster_path: string; input raster path
    :param channel_number: int; number of channel to be processed
    :param metadata_path: string; path to metadata file
    :param output_path: string; path for output radiance file
    :return:
    """
    try:
        int(channel_number)
    except:
        raise NameError('Channel number is not correct')

    if (int(channel_number)) < 1 or (int(channel_number) > 11):
        raise NameError('Channel number is not correct')

    metadata_channel_rad_maximum_str = 'RADIANCE_MAXIMUM_BAND_' + str(channel_number)
    metadata_channel_rad_minimum_str = 'RADIANCE_MINIMUM_BAND_' + str(channel_number)

    metadata_channel_quantize_max_str = 'QUANTIZE_CAL_MAX_BAND_' + str(channel_number)
    metadata_channel_quantize_min_str = 'QUANTIZE_CAL_MIN_BAND_' + str(channel_number)

    rad_maximum = float(read_metadata_parameter(metadata_path, metadata_channel_rad_maximum_str))
    rad_minimum = float(read_metadata_parameter(metadata_path, metadata_channel_rad_minimum_str))
    quantize_maximum = float(read_metadata_parameter(metadata_path, metadata_channel_quantize_max_str))
    quantize_minimum = float(read_metadata_parameter(metadata_path, metadata_channel_quantize_min_str))

    #print raster_path.replace('/','\\')
    landsat_dn_band = gdal.Open(re.sub("^\s+|\n|\r|\s+$", '', raster_path))
    #print landsat_dn_band
    landsat_dn_band_array = np.array(landsat_dn_band.GetRasterBand(1).ReadAsArray().astype(np.float32))
    print 'a???'
    landsat_radiance_band_array = ((rad_maximum - rad_minimum)/(quantize_maximum - quantize_minimum))*(landsat_dn_band_array - quantize_minimum) + rad_minimum
    print 'b???'
    if mode == 0:
        # Write result
        cols = landsat_dn_band.RasterXSize
        rows = landsat_dn_band.RasterYSize
        cell_type = gdal.GDT_Float32
        driver_name = 'GTiff'
        projection = landsat_dn_band.GetProjection()
        transform = landsat_dn_band.GetGeoTransform()
        save_nparray_as_raster(landsat_radiance_band_array,driver_name,cell_type,cols,rows,projection,transform,output_path)
        del landsat_radiance_band_array
        return
    else:
        return landsat_radiance_band_array 
Example 79
Project: unmixing   Author: arthur-e   File: utils.py    MIT License 4 votes vote down vote up
def pixel_to_xy(pixel_pairs, gt=None, wkt=None, path=None, dd=False):
    '''
    Modified from code by Zachary Bears (zacharybears.com/using-python-to-
    translate-latlon-locations-to-pixels-on-a-geotiff/).
    This method translates given pixel locations into longitude-latitude
    locations on a given GeoTIFF. Arguments:
        pixel_pairs The pixel pairings to be translated in the form
                    [[x1, y1],[x2, y2]]
        gt          [Optional] A GDAL GeoTransform tuple
        wkt         [Optional] Projection information as Well-Known Text
        path        The file location of the GeoTIFF
        dd          True to use decimal degrees for longitude-latitude (False
                    is the default)

    NOTE: This method does not take into account pixel size and assumes a
            high enough image resolution for pixel size to be insignificant.
    '''
    assert path is not None or (gt is not None and wkt is not None), 'Function requires either a reference dataset or a geotransform and projection'

    pixel_pairs = map(list, pixel_pairs)
    srs = osr.SpatialReference() # Create a spatial ref. for dataset

    if path is not None:
        ds = gdal.Open(path) # Load the image dataset
        gt = ds.GetGeoTransform() # Get geotransform of the dataset

        if dd:
            srs.ImportFromWkt(ds.GetProjection()) # Set up coord. transform.

    else:
        srs.ImportFromWkt(wkt)

    # Will use decimal-degrees if so specified
    if dd:
        ct = osr.CoordinateTransformation(srs, srs.CloneGeogCS())

    # Go through all the point pairs and translate them to pixel pairings
    xy_pairs = []
    for point in pixel_pairs:
        # Translate the pixel pairs into untranslated points
        lon = point[0] * gt[1] + gt[0]
        lat = point[1] * gt[5] + gt[3]
        if dd:
            (lon, lat, holder) = ct.TransformPoint(lon, lat) # Points to space

        xy_pairs.append((lon, lat)) # Add the point to our return array

    return xy_pairs 
Example 80
Project: unmixing   Author: arthur-e   File: utils.py    MIT License 4 votes vote down vote up
def xy_to_pixel(xy_pairs, gt=None, wkt=None, path=None, dd=False):
    '''
    Modified from code by Zachary Bears (zacharybears.com/using-python-to-
    translate-latlon-locations-to-pixels-on-a-geotiff/).
    This method translates given longitude-latitude pairs into pixel
    locations on a given GeoTIFF. Arguments:
        xy_pairs    The X-Y coordinate pairs (e.g., decimal lat/lon pairs) to
                    be translated in the form [[lon1, lat1], [lon2, lat2]]
        gt          [Optional] A GDAL GeoTransform tuple
        wkt         [Optional] Projection information as Well-Known Text
        path        The file location of the GeoTIFF
        dd          True to use decimal degrees for longitude-latitude (False
                    is the default)

    NOTE: This method does not take into account pixel size and assumes a
    high enough image resolution for pixel size to be insignificant.
    Essentially, xy_pairs assumed to be pixel centroids but may be assigned
    to the nearest neighbor of the intended pixel. For this reason, for
    relevant applications, target pixels should be embedded in a matrix
    of similar pixels, e.g., a 90-by-90 meter area, for Landsat pixels.
    '''
    assert path is not None or (gt is not None and wkt is not None), \
        'Function requires either a reference dataset or a geotransform and projection'

    xy_pairs = map(list, xy_pairs)
    srs = osr.SpatialReference() # Create a spatial ref. for dataset

    if path is not None:
        ds = gdal.Open(path) # Load the image dataset
        gt = ds.GetGeoTransform() # Get geotransform of the dataset

        if dd:
            srs.ImportFromWkt(ds.GetProjection()) # Set up coord. transform.

    else:
        srs.ImportFromWkt(wkt) # Set up coord. transform.

    # Will use decimal-degrees if so specified
    if dd:
        ct = osr.CoordinateTransformation(srs.CloneGeogCS(), srs)

    # Go through all the point pairs and translate them to lng-lat pairs
    pixel_pairs = []
    for point in xy_pairs:
        if dd:
            # Change the point locations into the GeoTransform space
            (point[0], point[1], holder) = ct.TransformPoint(point[0], point[1])

        # Translate the x and y coordinates into pixel values
        x = (point[0] - gt[0]) / gt[1]
        y = (point[1] - gt[3]) / gt[5]
        pixel_pairs.append((int(x), int(y))) # Add point to our return array

    return pixel_pairs