Python osgeo.osr.SpatialReference() Examples

The following are code examples for showing how to use osgeo.osr.SpatialReference(). 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 7 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 writeSRS(self, srs):
        """
        @summary: Writes spatial reference system information to this raster file.
        @param srs: An osgeo.osr.SpatialReference object or
                        A WKT string describing the desired spatial reference system.  
                          Raster.populateStats will populate the Raster.srs 
                          attribute with a correctly formatted string.  Correctly 
                          formatted strings are also output by:
                             * osgeo.gdal.Dataset.GetProjection
                             * osgeo.osr.SpatialReference.ExportToWKT 
        @postcondition: The raster file is updated with new srs information.
        @raise LMError: on failure to open dataset or write srs
        """
        from LmServer.common.geotools import GeoFileInfo

        if (self._dlocation is not None and os.path.exists(self._dlocation)):
            geoFI = GeoFileInfo(self._dlocation, updateable=True)
            if isinstance(srs, osr.SpatialReference):
                srs = srs.ExportToWkt()
            geoFI.writeWktSRS(srs)
                
# ............................................... 
Example 3
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def _transformBBox(self, origEpsg=None, origBBox=None):
        if origEpsg is None:
            origEpsg = self._epsg
        if origBBox is None:
            origBBox = self._bbox
        minX, minY, maxX, maxY = origBBox
            
        if origEpsg != DEFAULT_EPSG:
            srcSRS = osr.SpatialReference()
            srcSRS.ImportFromEPSG(origEpsg)
            dstSRS = osr.SpatialReference()
            dstSRS.ImportFromEPSG(DEFAULT_EPSG)
            
            spTransform = osr.CoordinateTransformation(srcSRS, dstSRS)
            # Allow for return of either (x, y) or (x, y, z) 
            retvals = spTransform.TransformPoint(minX, minY)
            newMinX, newMinY = retvals[0], retvals[1]
            retvals = spTransform.TransformPoint(maxX, maxY)
            newMaxX, newMaxY = retvals[0], retvals[1]
            return (newMinX, newMinY, newMaxX, newMaxY)
        else:
            return origBBox
        
# ............................................... 
Example 4
Project: GeoPy   Author: aerler   File: simple_regrid.py    GNU General Public License v3.0 6 votes vote down vote up
def __init__(self, lon, lat):
    '''
    initialize projection dataset based on (regular) lat/lon vectors
    '''    
    epsg = 4326 # EPSG code for regular lat/long grid
    # size of dataset
    size = (len(lon), len(lat))
    # GDAL geotransform vector
    dx = lon[1]-lon[0]; dy = lat[1]-lat[0]
    ulx = lon[0]-dx/2.; uly = lat[0]-dy/2. # coordinates of upper left corner (same for source and sink)
    # GT(2) & GT(4) are zero for North-up; GT(1) & GT(5) are pixel width and height; (GT(0),GT(3)) is the top left corner
    geotransform = (ulx, dx, 0., uly, 0., dy) 
    # GDAL projection 
    projection = osr.SpatialReference()
    projection.ImportFromEPSG(epsg)
    # create GDAL projection object from parent instance
    super(LatLonProj,self).__init__(projection=projection, geotransform=geotransform, size=size)
    self.epsg = epsg # save projection code number
    
## function to reproject and resample a 2D array 
Example 5
Project: GeoPy   Author: aerler   File: gdal.py    GNU General Public License v3.0 6 votes vote down vote up
def __setstate__(self, pickle):
    ''' support pickling, necessary for multiprocessing: GDAL is not pickable '''
    # handle projection
    self.projection = osr.SpatialReference() 
    self.projection.SetWellKnownGeogCS('WGS84')           
    self.projection.ImportFromWkt(pickle['_projection'])  # from Well-Known-Text
    del pickle['_projection'] # not actually an attribute
    # handle axes
    self.geotransform = pickle['_geotransform']
    self.isProjected = pickle['_isProjected']
    xlon, ylat = getAxes(geotransform=self.geotransform, xlen=pickle['_xlon'], ylen=pickle['_ylat'], 
                         projected=self.isProjected)
    self.xlon = xlon; self.ylat = ylat
    del pickle['_geotransform'], pickle['_isProjected'], pickle['_xlon'], pickle['_ylat']
    # update instance dict with pickle dict
    self.__dict__.update(pickle) 
Example 6
Project: buzzard   Author: airware   File: _a_source.py    Apache License 2.0 6 votes vote down vote up
def __init__(self, back_ds, wkt_stored, rect, **kwargs):
        wkt_virtual = back_ds.virtual_of_stored_given_mode(
            wkt_stored, back_ds.wkt_work, back_ds.wkt_fallback, back_ds.wkt_forced,
        )

        if wkt_virtual is not None:
            sr_virtual = osr.SpatialReference(wkt_virtual)
        else:
            sr_virtual = None

        to_work, to_virtual = back_ds.get_transforms(sr_virtual, rect)

        self.back_ds = back_ds
        self.wkt_stored = wkt_stored
        self.wkt_virtual = wkt_virtual
        self.to_work = to_work
        self.to_virtual = to_virtual
        super().__init__(**kwargs) 
Example 7
Project: buzzard   Author: airware   File: _dataset_back_conversions.py    Apache License 2.0 6 votes vote down vote up
def __init__(self, wkt_work, wkt_fallback, wkt_forced, analyse_transformation, **kwargs):

        if wkt_work is not None:
            sr_work = osr.SpatialReference(wkt_work)
        else:
            sr_work = None
        if wkt_fallback is not None:
            sr_fallback = osr.SpatialReference(wkt_fallback)
        else:
            sr_fallback = None
        if wkt_forced is not None:
            sr_forced = osr.SpatialReference(wkt_forced)
        else:
            sr_forced = None

        self.wkt_work = wkt_work
        self.wkt_fallback = wkt_fallback
        self.wkt_forced = wkt_forced
        self.sr_work = sr_work
        self.sr_fallback = sr_fallback
        self.sr_forced = sr_forced
        self.analyse_transformations = analyse_transformation
        super(BackDatasetConversionsMixin, self).__init__(**kwargs) 
Example 8
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 9
Project: unmixing   Author: arthur-e   File: utils.py    MIT License 6 votes vote down vote up
def get_coord_transform(source_epsg, target_epsg):
    '''
    Creates an OGR-framework coordinate transformation for use in projecting
    coordinates to a new coordinate reference system (CRS). Used as, e.g.:
        transform = get_coord_transform(source_epsg, target_epsg)
        transform.TransformPoint(x, y)
    Arguments:
        source_epsg     The EPSG code for the source CRS
        target_epsg     The EPSG code for the target CRS
    '''
    # Develop a coordinate transformation, if desired
    transform = None
    source_ref = osr.SpatialReference()
    target_ref = osr.SpatialReference()
    source_ref.ImportFromEPSG(source_epsg)
    target_ref.ImportFromEPSG(target_epsg)
    return osr.CoordinateTransformation(source_ref, target_ref) 
Example 10
Project: unmixing   Author: arthur-e   File: lsma.py    MIT License 6 votes vote down vote up
def get_idx_as_shp(self, path, gt, wkt):
        '''
        Exports a Shapefile containing the locations of the extracted
        endmembers. Assumes the coordinates are in decimal degrees.
        '''
        coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)

        driver = ogr.GetDriverByName('ESRI Shapefile')
        ds = driver.CreateDataSource(path)
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)

        layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint)
        for pair in coords:
            feature = ogr.Feature(layer.GetLayerDefn())

            # Create the point from the Well Known Text
            point = ogr.CreateGeometryFromWkt('POINT(%f %f)' %  pair)
            feature.SetGeometry(point) # Set the feature geometry
            layer.CreateFeature(feature) # Create the feature in the layer
            feature.Destroy() # Destroy the feature to free resources

        # Destroy the data source to free resources
        ds.Destroy() 
Example 11
Project: surfclass   Author: Kortforsyningen   File: test_rasterreader.py    MIT License 6 votes vote down vote up
def test_rasterreader_bboxrounding(classraster_filepath):
    reader = RasterReader(classraster_filepath)
    assert reader
    assert isinstance(reader.srs, osr.SpatialReference)
    assert reader.bbox == Bbox(727000.0, 6171000.0, 728000.0, 6172000.0)
    pix_win = reader.bbox_to_pixel_window(reader.bbox)
    assert pix_win == (0, 0, 500, 500)
    # Try bboxes which are slightly off pixel borders
    pix_win = reader.bbox_to_pixel_window(
        Bbox(727001.9999, 6171000, 728000.0, 6172000.0)
    )
    assert pix_win == (1, 0, 499, 500)
    pix_win = reader.bbox_to_pixel_window(Bbox(727000, 6171000, 727999.99, 6171999.99))
    assert pix_win == (0, 0, 500, 500)
    pix_win = reader.bbox_to_pixel_window(
        Bbox(727000, 6171000, 728000.00001, 6172000.00001)
    )
    assert pix_win == (0, 0, 500, 500) 
Example 12
Project: geo-mapgen   Author: Gael-de-Sailly   File: map_transform.py    GNU Lesser General Public License v2.1 6 votes vote down vote up
def get_map_bounds(mapname):
	if mapname in maps:
		thismap = maps[mapname]
	else:
		print("Map", mapname, "does not exist.")
		return
	xsize, ysize = thismap.RasterXSize, thismap.RasterYSize
	gt = thismap.GetGeoTransform()
	proj = osr.SpatialReference()
	proj.ImportFromWkt(thismap.GetProjection())
	transform = osr.CreateCoordinateTransformation(proj, wgs)
	minp = gm.transform(gt, (0, 0))
	maxp = gm.transform(gt, (xsize, ysize))
	xmin, ymin, _ = transform.TransformPoint(minp[0], minp[1])
	xmax, ymax, _ = transform.TransformPoint(maxp[0], maxp[1])
	north = max(ymin, ymax) # Maps might be reversed, so we're not sure which one is actually the maximum
	east = max(xmin, xmax)
	south = min(ymin, ymax)
	west = min(xmin, xmax)
	return (north, east, south, west) 
Example 13
Project: pangaea   Author: snowman2   File: xlsm.py    BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def to_projection(self, variable, projection):
        """Convert Grid to New Projection.

            Parameters
            ----------
            variable: :obj:`str`
                Name of variable in dataset.
            projection: :func:`osr.SpatialReference`
                Projection to convert data to.

            Returns
            -------
            :func:`xarray.Dataset`
        """
        new_data = []
        for band in range(self._obj.dims[self.time_dim]):
            arr_grid = ArrayGrid(in_array=self._obj[variable][band].values,
                                 wkt_projection=self.projection.ExportToWkt(),
                                 geotransform=self.geotransform)
            ggrid = arr_grid.to_projection(projection, gdalconst.GRA_Average)
            new_data.append(ggrid.np_array())

        self.to_datetime()
        return self._export_dataset(variable, np.array(new_data),
                                    ggrid) 
Example 14
Project: qgis-wps4server   Author: 3liz   File: UMN.py    GNU General Public License v3.0 6 votes vote down vote up
def getSpatialReference(self, output, datatype):
        """
        :param output: :class:`pywps.Process.InAndOutputs.ComplexOutput`
        :param datatype: String
        :return: projection of the output
        """

        sr = osr.SpatialReference()
        if datatype == "raster":
            wkt = self.dataset.GetProjection()
            res = sr.ImportFromWkt(wkt)
            if res == 0:
                return sr
        elif datatype == "vector":
            layer = self.dataset.GetLayer()
            ref = layer.GetSpatialRef()
            if ref:
                return ref
        return None 
Example 15
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 16
Project: core   Author: lifemapper   File: create_mask.py    GNU General Public License v3.0 5 votes vote down vote up
def write_tiff(out_filename, bbox, cell_size, data, epsg,
               nodata=DEFAULT_NODATA):
    """
    @summary: Write a GeoTiff raster layer from a numpy array
    @param out_filename: The file location to write the raster to
    @param bbox: (minx, miny, maxx, maxy) tuple of raster coordinates
    @param cell_size: The cell size, in map units, of each cell in the raster
    @param data: The data array to use for raster content
    @param epsg: The epsg code of the map projection to use for this raster
    @param nodata: A value to use for NODATA
    @todo: Probably should establish a common layer package in LmCommon
    """
    minx, miny, maxx, maxy = bbox
    (num_rows, num_cols) = data.shape
    
    readyFilename(out_filename, overwrite=True)
    
    drv = gdal.GetDriverByName(LMFormat.GTIFF.driver)
    ds = drv.Create(out_filename, num_cols, num_rows, 1, gdalconst.GDT_Byte)
    ds.SetGeoTransform([minx, cell_size, 0, maxy, 0, -cell_size])

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    ds.SetProjection(srs.ExportToWkt())
    
    outBand = ds.GetRasterBand(1)
    outBand.WriteArray(data)
    outBand.FlushCache()
    outBand.SetNoDataValue(nodata)

# .............................................................................
# Testing code 
Example 17
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 18
Project: core   Author: lifemapper   File: lmobj.py    GNU General Public License v3.0 5 votes vote down vote up
def createSRSFromEPSG(self, epsgcode=None):
        srs = SpatialReference()
        if epsgcode is not None:
            srs.ImportFromEPSG(epsgcode)
        else:
            srs.ImportFromEPSG(self._epsg)
        return srs

# ............................................... 
Example 19
Project: core   Author: lifemapper   File: lmobj.py    GNU General Public License v3.0 5 votes vote down vote up
def translatePoints(points, srcWKT=None, srcEPSG=None, dstWKT=None, dstEPSG=None):
        """
        @summary: Translates a list of (x, y) pairs from one coordinate system to 
                         another
        @param points: A list of points [(x, y)+]
        @param srcWKT: The Well-Known Text for the source coordinate system
        @param srcEPSG: The EPSG code of the source coordinate system
        @param dstWKT: The Well-Known Text for the destination coordinate system
        @param dstEPSG: The EPSG code of the destination coordinate system
        @note: Use either srcWKT or srcEPSG, srcWKT has priority
        @note: Use either dstWKT or dstEPSG, dstWKT has priority
        """
        srcSR = SpatialReference()
        if srcWKT is not None:
            srcSR.ImportFromWkt(srcWKT)
        elif srcEPSG is not None:
            srcSR.ImportFromEPSG(srcEPSG)
        else:
            raise Exception("Either srcWKT or srcEPSG must be specified")
    
        dstSR = SpatialReference()
        if dstWKT is not None:
            dstSR.ImportFromWkt(dstWKT)
        elif dstEPSG is not None:
            dstSR.ImportFromEPSG(dstEPSG)
        else:
            raise Exception("Either dstWKT or dstEPSG must be specified")
        
        transPoints = []
        
        trans = CoordinateTransformation(srcSR, dstSR)
        for pt in points:
            x, y, _ = trans.TransformPoint(pt[0], pt[1])
            transPoints.append((x, y))
        trans = None
        return transPoints
    


# ============================================================================ 
Example 20
Project: core   Author: lifemapper   File: createshape.py    GNU General Public License v3.0 5 votes vote down vote up
def _addUserFieldDef(self, newDataset):
        spRef = osr.SpatialReference()
        spRef.ImportFromEPSG(DEFAULT_EPSG)
        maxStrlen = LMFormat.getStrlenForDefaultOGR()
     
        newLyr = newDataset.CreateLayer('points', geom_type=ogr.wkbPoint, srs=spRef)
        if newLyr is None:
            raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                    'Layer creation failed')
        
        for idx, vals in self.op.columnMeta.iteritems():
            if vals is not None:
                fldname = str(vals[OccDataParser.FIELD_NAME_KEY])
                fldtype = vals[OccDataParser.FIELD_TYPE_KEY] 
                fldDef = ogr.FieldDefn(fldname, fldtype)
                if fldtype == ogr.OFTString:
                    fldDef.SetWidth(maxStrlen)
                returnVal = newLyr.CreateField(fldDef)
                if returnVal != 0:
                    raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                         'Failed to create field {}'.format(fldname))
                
        # Add wkt field
        fldDef = ogr.FieldDefn(LM_WKT_FIELD, ogr.OFTString)
        fldDef.SetWidth(maxStrlen)
        returnVal = newLyr.CreateField(fldDef)
        if returnVal != 0:
            raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                    'Failed to create field {}'.format(fldname))
        
        return newLyr
            

    # ............................................... 
Example 21
Project: EarthSim   Author: pyviz-topics   File: io.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def get_ccrs(filename):
    """
    Loads WKT projection string from file and return
    cartopy coordinate reference system.
    """
    inproj = osr.SpatialReference()
    proj = open(filename, 'r').readline()
    inproj.ImportFromWkt(proj)
    projcs = inproj.GetAuthorityCode('PROJCS')
    return ccrs.epsg(projcs) 
Example 22
Project: GeoPy   Author: aerler   File: gdal.py    GNU General Public License v3.0 5 votes vote down vote up
def getProjFromDict(projdict, name='', GeoCS='WGS84', convention='Proj4'):
  ''' Initialize a projected OSR SpatialReference instance from a dictionary using Proj4 conventions. 
      Valid parameters are documented here: http://trac.osgeo.org/proj/wiki/GenParms
      Projections are described here: http://www.remotesensing.org/geotiff/proj_list/ '''
  # initialize
  projection = osr.SpatialReference()
  # interpret dictionary according to convention
  if convention == 'Proj4':
    # start with projection, which is usually a string
    proj = projdict['proj']
    if proj == 'longlat' or proj == 'latlong':
      projstr = '+proj=latlong'; lproj = False
    else: 
      projstr = '+proj={0:s}'.format(proj); lproj = True 
    # loop over entries
    for key, value in projdict.items():
      if key is not 'proj':
        if not isinstance(key, str): raise TypeError
        if not isinstance(value, (float,np.inexact,int,np.integer)): raise TypeError
        # translate dict entries to string
        projstr = '{0:s} +{1:s}={2:f}'.format(projstr, key, value)
    # load projection from proj4 string
    projection.ImportFromProj4(projstr)
  else:
    raise NotImplementedError
  # more meta data
  if lproj: projection.SetProjCS(name)  # establish that this is a projected system
  # N.B.: note that the seemingly analogous method SetGeogCS() has a different function (and is not necessary)
  projection.SetWellKnownGeogCS(GeoCS)  # default reference datum/geoid
  # return finished projection object (geotransform can be inferred from coordinate vectors)
  return projection 
Example 23
Project: geophys_utils   Author: GeoscienceAustralia   File: _dem_utils.py    Apache License 2.0 5 votes vote down vote up
def get_pixel_size(self, index_tuple):
        """
        Returns X & Y sizes in metres of specified pixel as a tuple.
        N.B: Pixel ordinates are zero-based from top left
        """
        x, y = index_tuple
        logger.debug('(x, y) = (%f, %f)', x, y)
        
        native_spatial_reference = osr.SpatialReference()
        native_spatial_reference.ImportFromWkt(self.crs)
        
        latlong_spatial_reference = native_spatial_reference.CloneGeogCS()
        coord_transform_to_latlong = osr.CoordinateTransformation(native_spatial_reference, latlong_spatial_reference)

        # Determine pixel centre and edges in georeferenced coordinates
        xw = self.GeoTransform[0] + x * self.GeoTransform[1]
        yn = self.GeoTransform[3] + y * self.GeoTransform[5] 
        xc = self.GeoTransform[0] + (x + 0.5) * self.GeoTransform[1]
        yc = self.GeoTransform[3] + (y + 0.5) * self.GeoTransform[5] 
        xe = self.GeoTransform[0] + (x + 1.0) * self.GeoTransform[1]
        ys = self.GeoTransform[3] + (y + 1.0) * self.GeoTransform[5] 
        
        logger.debug('xw = %f, yn = %f, xc = %f, yc = %f, xe = %f, ys = %f', xw, yn, xc, yc, xe, ys)
        
        # Convert georeferenced coordinates to lat/lon for Vincenty
        lon1, lat1, _z = coord_transform_to_latlong.TransformPoint(xw, yc, 0)
        lon2, lat2, _z = coord_transform_to_latlong.TransformPoint(xe, yc, 0)
        logger.debug('For X size: (lon1, lat1) = (%f, %f), (lon2, lat2) = (%f, %f)', lon1, lat1, lon2, lat2)
        x_size, _az_to, _az_from = vinc_dist(earth.F, earth.A, 
                                             lat1 * RADIANS_PER_DEGREE, lon1 * RADIANS_PER_DEGREE, 
                                             lat2 * RADIANS_PER_DEGREE, lon2 * RADIANS_PER_DEGREE)
        
        lon1, lat1, _z = coord_transform_to_latlong.TransformPoint(xc, yn, 0)
        lon2, lat2, _z = coord_transform_to_latlong.TransformPoint(xc, ys, 0)
        logger.debug('For Y size: (lon1, lat1) = (%f, %f), (lon2, lat2) = (%f, %f)', lon1, lat1, lon2, lat2)
        y_size, _az_to, _az_from = vinc_dist(earth.F, earth.A, 
                                             lat1 * RADIANS_PER_DEGREE, lon1 * RADIANS_PER_DEGREE, 
                                             lat2 * RADIANS_PER_DEGREE, lon2 * RADIANS_PER_DEGREE)
        
        logger.debug('(x_size, y_size) = (%f, %f)', x_size, y_size)
        return (x_size, y_size) 
Example 24
Project: geophys_utils   Author: GeoscienceAustralia   File: _crs_utils.py    Apache License 2.0 5 votes vote down vote up
def get_spatial_ref_from_wkt(wkt_or_crs_name):
    '''
    Function to return SpatialReference object for supplied WKT
    @param wkt: Well-known text or CRS name for SpatialReference, including "EPSG:XXXX"
    @return spatial_ref: SpatialReference from WKT
    '''
    spatial_ref = SpatialReference()
    
    # Try to resolve WKT
    result = spatial_ref.ImportFromWkt(wkt_or_crs_name)
    if not result:
        return spatial_ref

    # Try to resolve CRS name - either mapped or original
    result = spatial_ref.SetWellKnownGeogCS(CRS_NAME_MAPPING.get(wkt_or_crs_name) or wkt_or_crs_name) 
    if not result:
        return spatial_ref

    # Try common formulations for UTM zones
    #TODO: Fix this so it works in the Northern hemisphere 
    modified_crs_name = re.sub('\s+', '', wkt_or_crs_name.strip().upper())
    utm_match = (re.match('(\w+)/MGAZONE(\d+)', modified_crs_name) or
                 re.match('(\w+)/(\d+)S', modified_crs_name) or
                 re.match('(EPSG:283)(\d{2})', modified_crs_name) 
                 )
    if utm_match:
        modified_crs_name = utm_match.group(1)
        utm_zone = int(utm_match.group(2))
        result = spatial_ref.SetWellKnownGeogCS(CRS_NAME_MAPPING.get(modified_crs_name) or modified_crs_name)
    if not result:
        spatial_ref.SetUTM(utm_zone, False) # Put this here to avoid potential side effects in downstream code
        return spatial_ref

    assert not result, 'Invalid WKT or CRS name' 
Example 25
Project: geophys_utils   Author: GeoscienceAustralia   File: _crs_utils.py    Apache License 2.0 5 votes vote down vote up
def get_utm_wkt(coordinate, from_wkt):
    '''
    Function to return CRS for UTM zone of specified coordinates.
    Used to transform coords to metres
    @param coordinate: single coordinate pair
    '''
    def utm_getZone(longitude):
        return (int(1 + (longitude + 180.0) / 6.0))

    def utm_isNorthern(latitude):
        if (latitude < 0.0):
            return 0
        else:
            return 1
        
    coordinate_array = np.array(coordinate).reshape((1,2))

    latlon_coord_trans = get_coordinate_transformation(
        from_wkt, 'EPSG:4326')
    latlon_coord = coordinate if latlon_coord_trans is None else latlon_coord_trans.TransformPoints(
        coordinate_array)[0][0:2]
        
    # Set UTM coordinate reference system
    utm_spatial_ref = SpatialReference()
    utm_spatial_ref.SetWellKnownGeogCS('WGS84')
    utm_spatial_ref.SetUTM(utm_getZone(
        latlon_coord[0]), utm_isNorthern(latlon_coord[1]))

    return utm_spatial_ref.ExportToWkt() 
Example 26
Project: PostgreSQL-GeoPackage   Author: EOX-A   File: gpkg-pg_dump.py    MIT License 5 votes vote down vote up
def create_gpkg(
    gpkg_name, proj_string, size=(1, 1), geotransform=[0, 1, 0, 0, 0, -1],
    creation_options=None
):
    if os.path.exists("%s.gpkg" % gpkg_name):
        sys.stderr.write(
            "ERROR: SQLite GeoPackage '%s.gpkg' already exists.\n" % gpkg_name
        )
        sys.exit(1)

    gdal.AllRegister()
    drv = gdal.GetDriverByName("GPKG")
    try:
        gpkg = drv.Create(
            "%s.gpkg" % gpkg_name, size[0], size[1], 1, gdal.GDT_Byte,
            creation_options
        )

        proj = osr.SpatialReference()
        res = proj.SetWellKnownGeogCS(proj_string)
        if res != 0:
            if proj_string[0:4] == 'EPSG':
                proj.ImportFromEPSG(int(proj_string[5:]))
        gpkg.SetProjection(proj.ExportToWkt())
        gpkg.SetGeoTransform(geotransform)
        gpkg = None
    except Exception as e:
        sys.stderr.write(
            "ERROR: Cannot create SQLite GeoPackage '%s.gpkg'. "
            "Error message was: '%s'.\n" % (gpkg_name, e.message)
        )
        sys.exit(1) 
Example 27
Project: buzzard   Author: airware   File: _a_source.py    Apache License 2.0 5 votes vote down vote up
def proj4_virtual(self):
        if self.wkt_virtual is None:
            return None # pragma: no cover
        return osr.SpatialReference(self.wkt_virtual).ExportToProj4() 
Example 28
Project: buzzard   Author: airware   File: _a_source.py    Apache License 2.0 5 votes vote down vote up
def proj4_stored(self):
        if self.wkt_stored is None:
            return None # pragma: no cover
        return osr.SpatialReference(self.wkt_stored).ExportToProj4() 
Example 29
Project: buzzard   Author: airware   File: __init__.py    Apache License 2.0 5 votes vote down vote up
def wkt_same(a, b):
    """Are two wkt equivalent"""
    if a == b:
        return True
    sra = osr.SpatialReference(a)
    srb = osr.SpatialReference(b)
    return bool(sra.IsSame(srb)) 
Example 30
Project: buzzard   Author: airware   File: tools.py    Apache License 2.0 5 votes vote down vote up
def sreq(*items, verbose_on=None):
    """SR items are all equal"""
    v = tuple(map(int, gdal.__version__.split('.')))
    for a, b in itertools.combinations(items, 2):
        a = osr.SpatialReference(osr.GetUserInputAsWKT(a))
        if v < (3,):
            a.StripCTParms()
        a = a.ExportToProj4()
        a = osr.SpatialReference(osr.GetUserInputAsWKT(a))
        if v < (3,):
            a.StripCTParms()

        b = osr.SpatialReference(osr.GetUserInputAsWKT(b))
        if v < (3,):
            b.StripCTParms()
        b = b.ExportToProj4()
        b = osr.SpatialReference(osr.GetUserInputAsWKT(b))
        if v < (3,):
            b.StripCTParms()

        res = bool(a.IsSame(b))
        if res is verbose_on:
            print('')
            print(a.ExportToPrettyWkt())
            print('---------- vs ----------')
            print(b.ExportToPrettyWkt())
            print('')

        if not res:
            return False
    return True 
Example 31
Project: buzzard   Author: airware   File: _dataset.py    Apache License 2.0 5 votes vote down vote up
def proj4(self):
        """Dataset's work spatial reference in WKT proj4.
        Returns None if `mode 1`.
        """
        if self._back.wkt_work is None:
            return None
        return osr.SpatialReference(self._back.wkt_work).ExportToProj4() 
Example 32
Project: press-alt   Author: mrwojtek   File: elevation.py    Apache License 2.0 5 votes vote down vote up
def __init__(self, file, band=1):
        self.dataset = gdal.Open(file)
        self.geotransform = self.dataset.GetGeoTransform()
        
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.dataset.GetProjection())
        self.proj = pyproj.Proj(srs.ExportToProj4())
        
        band = self.dataset.GetRasterBand(band) 
Example 33
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 34
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 35
Project: surface   Author: africker   File: surface.py    MIT License 5 votes vote down vote up
def read(self, infile):
		# never want to change the original input raster so use read only constant
		self.infile = infile
		self.raster = gdal.Open(self.infile, GA_ReadOnly)
		self.band = self.raster.GetRasterBand(1)
		self.NDV = self.band.GetNoDataValue()
		self.x = self.band.XSize
		self.y = self.band.YSize
		DataType = self.band.DataType
		self.DataType = gdal.GetDataTypeName(DataType)
		self.GeoT = self.raster.GetGeoTransform()
		# Problem
		self.prj = osr.SpatialReference()
		self.prj.ImportFromWkt(self.raster.GetProjectionRef()) 
Example 36
Project: lidar   Author: giswqs   File: filling.py    MIT License 5 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 37
Project: eo_tools   Author: neel9102   File: extract_mod09q1_veg_index_by_coord.py    GNU General Public License v3.0 5 votes vote down vote up
def wgs84_to_modissinu ( lon, lat ):
	wgs84 = osr.SpatialReference( )  # Define a SpatialReference object
	wgs84.ImportFromEPSG( 4326 )
	# Define the modis data projection, sinusoidal #sinu_modis.ImportFromEPSG( 27700)
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(wgs84, modis_sinu)
	modis_x, modis_y, modis_z = tx.TransformPoint ( lon, lat )
	return modis_x, modis_y, modis_z 
Example 38
Project: eo_tools   Author: neel9102   File: extract_mod09q1_veg_index_by_coord.py    GNU General Public License v3.0 5 votes vote down vote up
def modissinu_to_wgs84 ( modis_x, modis_y ):
	from osgeo import osr
	wgs84 = osr.SpatialReference( )
	wgs84.ImportFromEPSG( 4326 )
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(modis_sinu, wgs84)
	lon, lat, modis_z = tx.TransformPoint ( modis_x, modis_y ) 
Example 39
Project: eo_tools   Author: neel9102   File: interpolate_mod09q1_veg_index_by_coord_release.py    GNU General Public License v3.0 5 votes vote down vote up
def wgs84_to_modissinu ( lon, lat ):
	wgs84 = osr.SpatialReference( )
	wgs84.ImportFromEPSG( 4326 )
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(wgs84, modis_sinu)
	modis_x, modis_y, modis_z = tx.TransformPoint ( lon, lat )
	return modis_x, modis_y, modis_z 
Example 40
Project: eo_tools   Author: neel9102   File: interpolate_mod09q1_veg_index_by_coord_release.py    GNU General Public License v3.0 5 votes vote down vote up
def modissinu_to_wgs84 ( modis_x, modis_y ):
	from osgeo import osr
	wgs84 = osr.SpatialReference( )
	wgs84.ImportFromEPSG( 4326 )
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(modis_sinu, wgs84)
	lon, lat, modis_z = tx.TransformPoint ( modis_x, modis_y ) 
Example 41
Project: eo_tools   Author: neel9102   File: L2_LAC_SST_preprocessing.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array):
  
    driver = gdal.GetDriverByName("GTiff")
    outDataset = driver.Create(output_name, ncols, nrows, 1, gdal.GDT_Float32)
    outDataset.SetGeoTransform(geo_transform_1)
    #srs = osr.SpatialReference()
    #srs.ImportFromEPSG(4326) #WGS84 lat long.
    outDataset.SetProjection(proj_1)
    
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    
    outBand.SetNoDataValue(np.nan)
    outBand.FlushCache() 
Example 42
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 5 votes vote down vote up
def preprocess_sen2_images(l2_dir, out_dir, l1_dir, cloud_threshold=60, buffer_size=0, epsg=None,
                           bands=("B08", "B04", "B03", "B02"), out_resolution=10):
    """For every .SAFE folder in in_dir, stacks band 2,3,4 and 8  bands into a single geotif, creates a cloudmask from
    the combined fmask and sen2cor cloudmasks and reprojects to a given EPSG if provided"""
    safe_file_path_list = [os.path.join(l2_dir, safe_file_path) for safe_file_path in os.listdir(l2_dir)]
    for l2_safe_file in safe_file_path_list:
        with TemporaryDirectory() as temp_dir:
            log.info("----------------------------------------------------")
            log.info("Merging 10m bands in SAFE dir: {}".format(l2_safe_file))
            temp_path = os.path.join(temp_dir, get_sen_2_granule_id(l2_safe_file)) + ".tif"
            log.info("Output file: {}".format(temp_path))
            stack_sentinel_2_bands(l2_safe_file, temp_path, bands=bands, out_resolution=out_resolution)

            log.info("Creating cloudmask for {}".format(temp_path))
            l1_safe_file = get_l1_safe_file(l2_safe_file, l1_dir)
            mask_path = get_mask_path(temp_path)
            create_mask_from_sen2cor_and_fmask(l1_safe_file, l2_safe_file, mask_path, buffer_size=buffer_size)
            log.info("Cloudmask created")

            out_path = os.path.join(out_dir, os.path.basename(temp_path))
            out_mask_path = os.path.join(out_dir, os.path.basename(mask_path))

            if epsg:
                log.info("Reprojecting images to {}".format(epsg))
                proj = osr.SpatialReference()
                proj.ImportFromEPSG(epsg)
                wkt = proj.ExportToWkt()
                reproject_image(temp_path, out_path, wkt)
                reproject_image(mask_path, out_mask_path, wkt)
            else:
                log.info("Moving images to {}".format(out_dir))
                shutil.move(temp_path, out_path)
                shutil.move(mask_path, out_mask_path) 
Example 43
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    GNU General Public License v3.0 5 votes vote down vote up
def reproject_geotransform(in_gt, old_proj_wkt, new_proj_wkt):
    """
    Reprojects a geotransform from the old projection to a new projection. See
    [https://gdal.org/user/raster_data_model.html]

    Parameters
    ----------
    in_gt
        A six-element numpy array, usually an output from gdal_image.GetGeoTransform()
    old_proj_wkt
        The projection of the old geotransform in well-known text.
    new_proj_wkt
        The projection of the new geotrasform in well-known text.
    Returns
    -------
    out_gt
        The geotransform in the new projection

    """
    old_proj = osr.SpatialReference()
    new_proj = osr.SpatialReference()
    old_proj.ImportFromWkt(old_proj_wkt)
    new_proj.ImportFromWkt(new_proj_wkt)
    transform = osr.CoordinateTransformation(old_proj, new_proj)
    (ulx, uly, _) = transform.TransformPoint(in_gt[0], in_gt[3])
    out_gt = (ulx, in_gt[1], in_gt[2], uly, in_gt[4], in_gt[5])
    return out_gt 
Example 44
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    GNU General Public License v3.0 5 votes vote down vote up
def write_geometry(geometry, out_path, srs_id=4326):
    """
    Saves the geometry in an ogr.Geometry object to a shapefile.

    Parameters
    ----------
    geometry
        An ogr.Geometry object
    out_path
        The location to save the output shapefile
    srs_id
        The projection of the output shapefile. Can be an EPSG number or a WKT string.

    Notes
    -----
    The shapefile consists of one layer named 'geometry'.


    """
    # TODO: Fix this needing an extra filepath on the end
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    if type(srs_id) is int:
        srs.ImportFromEPSG(srs_id)
    if type(srs_id) is str:
        srs.ImportFromWkt(srs_id)
    layer = data_source.CreateLayer(
        "geometry",
        srs,
        geom_type=geometry.GetGeometryType())
    feature_def = layer.GetLayerDefn()
    feature = ogr.Feature(feature_def)
    feature.SetGeometry(geometry)
    layer.CreateFeature(feature)
    data_source.FlushCache()
    data_source = None 
Example 45
Project: gis   Author: xiaoyingpu   File: readerclass.py    MIT License 5 votes vote down vote up
def __str__(self):
        """
        returns the projection spec.
        """
        pszProjection = self.data_set.GetProjectionRef()
        if pszProjection is not None:
            hSRS = osr.SpatialReference()
            if hSRS.ImportFromWkt(pszProjection ) == gdal.CE_None:
                pszPrettyWkt = hSRS.ExportToPrettyWkt(False)
                return pszPrettyWkt
            else:
                return pszProjection 
Example 46
Project: surfclass   Author: Kortforsyningen   File: rasterio.py    MIT License 5 votes vote down vote up
def srs(self):
        """SpatialReference: Spatial reference system used by the loaded raster file."""
        if self._srs is None:
            srs = osr.SpatialReference()
            srs.ImportFromWkt(self._ds.GetProjection())
            self._srs = srs
        return self._srs 
Example 47
Project: surfclass   Author: Kortforsyningen   File: options.py    MIT License 5 votes vote down vote up
def srs_handler(ctx, param, value):
    out_srs = osr.SpatialReference()
    out_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    if out_srs.SetFromUserInput(value) != 0:
        raise ValueError("Failed to process SRS definition: %s" % value)
    return out_srs 
Example 48
Project: surfclass   Author: Kortforsyningen   File: test_rasterreader.py    MIT License 5 votes vote down vote up
def test_maskedrasterreader(classraster_filepath):
    reader = MaskedRasterReader(classraster_filepath)
    assert reader
    assert isinstance(reader.srs, osr.SpatialReference)
    assert reader.bbox == Bbox(727000.0, 6171000.0, 728000.0, 6172000.0)
    pix_win = reader.bbox_to_pixel_window(reader.bbox)
    assert pix_win == (0, 0, 500, 500)
    pix_win = reader.bbox_to_pixel_window(
        Bbox(727500.0, 6171600.0, 727600.0, 6171750.0)
    )
    assert pix_win == (250, 125, 50, 75)
    # Try using a point
    pnt = ogr.Geometry(ogr.wkbPoint)
    pnt.AddPoint(727500.0, 6171600.0)
    pnt.AssignSpatialReference(reader.srs)
    data = reader.read_2d(pnt)
    np.testing.assert_array_equal(data, np.ma.empty(shape=(0, 0)))
    # Now use a proper polygon
    poly = pnt.Buffer(100)
    data = reader.read_2d(poly)
    assert data.shape == (100, 100)
    # Check that we have a mask
    assert int(np.sum(data.mask)) == 2140
    # Check data without mask
    assert int(np.sum(data.data)) == 25412
    # Check data with mask. (Must be less than unmasked)
    assert int(np.sum(data.compressed())) == 20432

    # Test read_flattened
    flat_data = reader.read_flattened(poly)
    # read_2d returns a 100x100 of which 2140 are outside the poly
    assert flat_data.shape == (100 * 100 - 2140,)
    assert np.ma.is_masked(flat_data)
    assert int(np.ma.sum(flat_data)) == 20432 
Example 49
Project: bitcoin-arbitrage-trading-bot   Author: mammuth   File: jqvmap.py    MIT License 5 votes vote down vote up
def __init__(self, config):
    default_config = {
      "projection": "merc",
      "longitude0": 0
    }
    default_config.update(config)
    self.config = default_config

    self.spatialRef = osr.SpatialReference()
    projString = '+proj='+str(self.config['projection'])+' +a=6381372 +b=6381372 +lat_0=0'
    #if 'emulate_longitude0' in self.config and not self.config['emulate_longitude0']:
    projString += ' +lon_0='+str(self.config['longitude0'])
    self.spatialRef.ImportFromProj4(projString) 
Example 50
Project: geo-mapgen   Author: Gael-de-Sailly   File: map_transform.py    GNU Lesser General Public License v2.1 5 votes vote down vote up
def update_map(mapname, newfilepath, get_proj=None):
	if mapname in maps_paths and maps_paths[mapname] == newfilepath:
		return
	dataset = gdal.Open(newfilepath)
	if dataset:
		maps[mapname] = dataset
		maps_paths[mapname] = newfilepath
		proj = osr.SpatialReference()
		proj.ImportFromWkt(dataset.GetProjection())
		if proj.Validate() != 0 and get_proj:
			epsg = get_proj(mapname)
			proj.ImportFromEPSG(epsg)
			dataset.SetProjection(proj.ExportToWkt())