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: 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 2
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 3
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 4
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 5
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 6
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 7
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 8
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 9
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 10
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 11
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 12
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 13
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 14
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 15
Project: lidar   Author: giswqs   File: slicing.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])
    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 16
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 17
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 18
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 19
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 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.float32)
    return band,geoTransform,proj,ncol,nrow
    dataset = None
    band = None 
Example 21
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 22
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 23
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 24
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 25
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 26
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 27
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 28
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 29
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 30
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 31
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 32
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 33
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 34
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 35
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 36
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 37
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 38
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 39
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 40
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 41
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 42
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 43
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 44
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 45
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 46
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 47
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 48
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 49
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 50
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)