Python QgsCoordinateReferenceSystem.createFromEpsg方法代码示例

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
def viewport_geo_array(map_canvas):
    """Obtain the map canvas current extent in EPSG:4326.

    :param map_canvas: A map canvas instance.
    :type map_canvas: QgsMapCanvas

    :returns: A list in the form [xmin, ymin, xmax, ymax] where all
        coordinates provided are in Geographic / EPSG:4326.
    :rtype: list

    .. note:: Delegates to extent_to_geo_array()
    """

    # get the current viewport extent
    myRect = map_canvas.extent()

    if map_canvas.hasCrsTransformEnabled():
        myCrs = map_canvas.mapRenderer().destinationCrs()
    else:
        # some code duplication from extentToGeoArray here
        # in favour of clarity of logic...
        myCrs = QgsCoordinateReferenceSystem()
        myCrs.createFromEpsg(4326)

    return extent_to_geo_array(myRect, myCrs)

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
def extentToGeoArray(theExtent, theSourceCrs):
    """Convert the supplied extent to geographic and return as as array.

    Args:
        theExtent: QgsRectangle to be transformed to geocrs.
        theSourceCrs: QgsCoordinateReferenceSystem representing the original
            extent's CRS.

    Returns:
        list: Transformed extents in EPSG:4326 in the form
              [xmin, ymin, xmix, ymax]

    Raises:
        None
    """

    myGeoCrs = QgsCoordinateReferenceSystem()
    myGeoCrs.createFromEpsg(4326)
    myXForm = QgsCoordinateTransform(
        theSourceCrs,
        myGeoCrs)

    # Get the clip area in the layer's crs
    myTransformedExtent = myXForm.transformBoundingBox(theExtent)

    myGeoExtent = [myTransformedExtent.xMinimum(),
                   myTransformedExtent.yMinimum(),
                   myTransformedExtent.xMaximum(),
                   myTransformedExtent.yMaximum()]
    return myGeoExtent

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
    def _extentsToLayer(self):
        """Memory layer for aggregation by using canvas extents as feature.

        We do this because the user elected to use no aggregation layer so we
        make a 'dummy' one which covers the whole study area extent.

        This layer is needed when postprocessing because we always want a
        vector layer to store aggregation information in.

        Returns:
            QgsMapLayer - a memory layer representing the extents of the clip.
        """

        # Note: this code duplicates from Dock.viewportGeoArray - make DRY. TS

        myRect = self.iface.mapCanvas().extent()
        myCrs = QgsCoordinateReferenceSystem()
        myCrs.createFromEpsg(4326)
        myGeoExtent = extent_to_geo_array(myRect, myCrs)

        if not self.layer.isValid():
            myMessage = self.tr(
                'An exception occurred when creating the entire area layer.')
            raise (Exception(myMessage))

        myProvider = self.layer.dataProvider()

        myAttrName = self.tr('Area')
        myProvider.addAttributes(
            [QgsField(myAttrName, QtCore.QVariant.String)])

        self.layer.startEditing()
        # add a feature the size of the impact layer bounding box
        myFeature = QgsFeature()
        # noinspection PyCallByClass,PyTypeChecker,PyArgumentList
        myFeature.setGeometry(QgsGeometry.fromRect(
            QgsRectangle(
                QgsPoint(myGeoExtent[0], myGeoExtent[1]),
                QgsPoint(myGeoExtent[2], myGeoExtent[3]))))
        myFeature.setAttributeMap({0: QtCore.QVariant(
            self.tr('Entire area'))})
        myProvider.addFeatures([myFeature])
        self.layer.commitChanges()

        try:
            self.keywordIO.update_keywords(
                self.layer,
                {self.defaults['AGGR_ATTR_KEY']: myAttrName})
        except InvalidParameterError:
            self.keywordIO.write_keywords(
                self.layer,
                {self.defaults['AGGR_ATTR_KEY']: myAttrName})
        except KeywordDbError, e:
            raise e

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
 def zoomTo(self, name):
     for res in self.results:
         if unicode(res.description) == unicode(name):
             dest_crs = self.canvas.mapRenderer().destinationCrs()
             src_crs = QgsCoordinateReferenceSystem()
             src_crs.createFromEpsg(res.epsg)
             transform = QgsCoordinateTransform(src_crs, dest_crs)
             new_point = transform.transform(res.x, res.y)
             x = new_point.x()
             y = new_point.y()
             self.canvas.setExtent(QgsRectangle(x,y,x,y))
             self.canvas.zoomScale(res.zoom)
             self.canvas.refresh()
             self.marker.setCenter(new_point)
             self.marker.show()
             return

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
def setCanvasCrs(theEpsgId, theOtfpFlag=False):
    """Helper to set the crs for the CANVAS before a test is run.

    Args:

        * theEpsgId  - Valid EPSG identifier (int)
        * theOtfpFlag - whether on the fly projections should be enabled
                        on the CANVAS. Default to False.
    """
        # Enable on-the-fly reprojection
    CANVAS.mapRenderer().setProjectionsEnabled(theOtfpFlag)

    # Create CRS Instance
    myCrs = QgsCoordinateReferenceSystem()
    myCrs.createFromEpsg(theEpsgId)  # google mercator

    # Reproject all layers to WGS84 geographic CRS
    CANVAS.mapRenderer().setDestinationCrs(myCrs)

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
 def visible(self,value):
     if value == True:
         if self.x is not None and self.y is not None:
             self._visible = True
             dest_crs = self.canvas.mapRenderer().destinationCrs()
             src_crs = QgsCoordinateReferenceSystem()
             src_crs.createFromEpsg(self.epsg)
             transform = QgsCoordinateTransform(src_crs, dest_crs)
             new_point = transform.transform(self.x, self.y)
             self._xtrans = new_point.x()
             self._ytrans = new_point.y()
             self.marker.setCenter(new_point)
             self.marker.show()
         else:
             self._visible = False
             raise ValueError("Can't show marker without x and y coordinates.")
     else:
         self._visible = False
         self.marker.hide()

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
def getWGS84resolution(theLayer, theGeoExtent=None):
    """Return resolution of raster layer in EPSG:4326

    Input
        theLayer: Raster layer
        theGeoExtent: Bounding box in EPSG:4326
        # FIXME (Ole), the second argumunt should be obtained within
                       this function to make it independent
    Output
        resolution.

    If input layer is already in EPSG:4326, simply return the resolution
    If not, work it out based on EPSG:4326 representations of its extent
    """

    msg = tr('Input layer to getWGS84resolution must be a raster layer. '
           'I got: %s' % str(theLayer.type())[1:-1])
    if not theLayer.type() == QgsMapLayer.RasterLayer:
        raise RuntimeError(msg)

    if theLayer.crs().authid() == 'EPSG:4326':
        # If it is already in EPSG:4326, simply use the native resolution
        myCellSize = theLayer.rasterUnitsPerPixel()
    else:
        # Otherwise, work it out based on EPSG:4326 representations of
        # its extent

        # Reproject extent to EPSG:4326
        myGeoCrs = QgsCoordinateReferenceSystem()
        myGeoCrs.createFromEpsg(4326)

        # Estimate cellsize
        # FIXME (Ole): Get geoextent from layer
        myColumns = theLayer.width()
        myGeoWidth = abs(theGeoExtent[3] - theGeoExtent[0])
        myCellSize = myGeoWidth / myColumns

    return myCellSize

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
    def _createPolygonLayer(self, crs=None, fields=None):
        """Creates an empty shape file layer"""

        if crs is None:
            crs = QgsCoordinateReferenceSystem()
            crs.createFromEpsg(4326)

        if fields is None:
            fields = {}

        myTempdir = temp_dir(sub_dir='preprocess')
        myOutFilename = unique_filename(suffix='.shp',
                                        dir=myTempdir)
        mySHPWriter = QgsVectorFileWriter(myOutFilename,
                                          'UTF-8',
                                          fields,
                                          QGis.WKBPolygon,
                                          crs)
        #flush the writer to write to file
        del mySHPWriter
        myName = self.tr('Entire area')
        myLayer = QgsVectorLayer(myOutFilename, myName, 'ogr')
        LOGGER.debug('created' + myLayer.name())
        return myLayer

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
 def coordRefSys(self, mapCoordSys):
     epsg = self.epsgList[0]
     coordRefSys = QgsCoordinateReferenceSystem()
     if QGis.QGIS_VERSION_INT >= 10900:
         idEpsgRSGoogle = "EPSG:%d" % epsg
         createCrs = coordRefSys.createFromOgcWmsCrs(idEpsgRSGoogle)
     else:
         idEpsgRSGoogle = epsg
         createCrs = coordRefSys.createFromEpsg(idEpsgRSGoogle)
     if not createCrs:
         google_proj_def = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 "
         google_proj_def += "+units=m [email protected] +wktext +no_defs"
         isOk = coordRefSys.createFromProj4(google_proj_def)
         if not isOk:
             return None
     return coordRefSys

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
 def coordRefSys(self, mapCoordSys):
     epsg = self.epsgList[0]
     coordRefSys = QgsCoordinateReferenceSystem()
     if QGis.QGIS_VERSION_INT >= 10900:
         idEpsgRSGoogle = "EPSG:%d" % epsg
         createCrs = coordRefSys.createFromOgcWmsCrs(idEpsgRSGoogle)
     else:
         idEpsgRSGoogle = epsg
         createCrs = coordRefSys.createFromEpsg(idEpsgRSGoogle)
     if not createCrs:
         proj_def =  "+proj=tmerc +lat_0=38 +lon_0=127.5 +k=0.9996 +x_0=1000000 +y_0=2000000 +ellps=GRS80 "
         proj_def += "+towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
         isOk = coordRefSys.createFromProj4(proj_def)
         if not isOk:
             return None
     return coordRefSys

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
class ToGrid(Operator):
    [STAT_AREA_IDX, STAT_COUNT_IDX] = range(2)

    def __init__(self, options=None, name="Base Abstract Grid Operator"):
        super(ToGrid, self).__init__(options, name)
        self._tmp_dir = options["tmp_dir"]

        self.mercator_crs = QgsCoordinateReferenceSystem()
        self.mercator_crs.createFromEpsg(3395)

    @property
    def output_types(self):
        return [OperatorDataTypes.Grid, OperatorDataTypes.Shapefile]

    @property
    def output_names(self):
        return ["Grid with Zone and Count Attributes", "Grid Shapefile"]

    output_descriptions = output_names

    # protected methods override
    ###########################
    def _verify_inputs(self, inputs):
        pass

    def _verify_outputs(self, outputs):
        pass

    # protected methods
    ###########################
    def _outputGeometryType(self):
        return QGis.WKBPolygon

    def _outputGeometryFromLatLon(self, lat, lon):
        # conversion lat/lon to grid_id is required
        # it acts as rounding function to lat/lon
        return self._outputGeometryFromGridId(latlon_to_grid(lat, lon))

    def _outputGeometryFromGridId(self, grid_id):
        [lat, lon] = grid_to_latlon(int(grid_id))
        return QgsGeometry.fromRect(
            QgsRectangle(
                lon - DEFAULT_HALF_GRID_SIZE,
                lat - DEFAULT_HALF_GRID_SIZE,
                lon + DEFAULT_HALF_GRID_SIZE,
                lat + DEFAULT_HALF_GRID_SIZE,
            )
        )

    def _create_grid(self, grid_name, grid_file, x_min, y_min, x_max, y_max, x_off, y_off):
        x_off2, y_off2 = x_off / 2.0, y_off / 2.0
        x_min = floor(x_min / x_off) * x_off
        x_max = ceil(x_max / x_off) * x_off
        y_min = floor(y_min / y_off) * y_off
        y_max = ceil(y_max / y_off) * y_off

        xtotal = int((x_max - x_min) / x_off) + 1
        ytotal = int((y_max - y_min) / y_off) + 1

        logAPICall.log(
            "x_min %f x_max %f y_min %f y_max %f x_off %f y_off %f xtotal %d, ytotal %d"
            % (x_min, x_max, y_min, y_max, x_off, y_off, xtotal, ytotal),
            logAPICall.DEBUG_L2,
        )
        fields = {0: QgsField("GRID_GID", QVariant.String)}
        writer = QgsVectorFileWriter(grid_file, "utf-8", fields, QGis.WKBPolygon, self._crs, "ESRI Shapefile")
        f = QgsFeature()
        for x in range(xtotal):
            for y in range(ytotal):
                lon = x_min + (x * x_off) + (x_off2)
                lat = y_min + (y * y_off) + (y_off2)
                # out_geom = QgsGeometry.fromRect(QgsRectangle(lon-x_off2, lat-y_off2,
                #                                             lon+x_off2, lat+y_off2))
                f.setGeometry(self._outputGeometryFromLatLon(lat, lon))
                f.addAttribute(0, QVariant(latlon_to_grid(lat, lon)))
                writer.addFeature(f)
        del writer
        return load_shapefile(grid_file, grid_name)

    def _create_zone_statistics(self, zone_layer, zone_field, count_field, zone_stat, zone_names):
        # project geometry into mercator and get area in m2
        mercator_transform = QgsCoordinateTransform(zone_layer.crs(), self.mercator_crs)
        zone_gid_idx = layer_field_index(zone_layer, GID_FIELD_NAME)
        zone_field_idx = layer_field_index(zone_layer, zone_field)
        count_field_idx = layer_field_index(zone_layer, count_field)
        for _f in layer_features(zone_layer):
            # project into mercator and get area in m2
            geom = _f.geometry()
            geom.transform(mercator_transform)
            area = geom.area()

            gid = _f.attributeMap()[zone_gid_idx].toString()
            # if count field is not defined, then set count to 0
            if count_field_idx >= 0:
                count = _f.attributeMap()[count_field_idx].toDouble()[0]
            else:
                count = 0

            self._update_stat(zone_stat, gid, count, area)
            zone_names[gid] = _f.attributeMap()[zone_field_idx]
#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
    def do_operation(self):
        """ perform create mapping scheme operation """

        # validate inputs
        fp_layer = self.inputs[0].value
        zone_layer = self.inputs[1].value
        zone_field = self.inputs[2].value
        count_field = self.inputs[3].value
        area_field = self.inputs[4].value

        # make sure input is correct
        # NOTE: these checks cannot be performed at set input time
        #       because the data layer maybe is not loaded yet
        self._test_layer_loaded(fp_layer)
        self._test_layer_loaded(zone_layer)
        self._test_layer_field_exists(zone_layer, GID_FIELD_NAME)
        self._test_layer_field_exists(zone_layer, zone_field)
        # count_field is not required
        # if count field is not defined, then generate building count from footprints
        # area_field is not required

        # local variables
        analyzer = QgsOverlayAnalyzer()
        area_idx = ToGrid.STAT_AREA_IDX
        cnt_idx = ToGrid.STAT_COUNT_IDX

        zone_names, zone_stat, zone_stat2, zone_totals = {}, {}, {}, {}

        # 1. find building count and total area for each zone
        # project geometry into mercator and get area in m2
        mercator_crs = QgsCoordinateReferenceSystem()
        mercator_crs.createFromEpsg(3395)
        mercator_transform = QgsCoordinateTransform(zone_layer.crs(), mercator_crs)

        try:
            # use zone geometry area
            self._create_zone_statistics(zone_layer, zone_field, count_field, zone_stat, zone_names)
        except Exception as err:
            raise OperatorError(str(err), self.__class__)

        # 2. create grids around extent of zone
        tmp_grid1 = "grid_" + get_unique_filename()
        tmp_grid1_file = self._tmp_dir + tmp_grid1 + ".shp"
        extent = zone_layer.extent()
        [x_min, y_min, x_max, y_max] = [extent.xMinimum(), extent.yMinimum(), extent.xMaximum(), extent.yMaximum()]
        tmp_grid_lyr1 = self._create_grid(
            tmp_grid1, tmp_grid1_file, x_min, y_min, x_max, y_max, DEFAULT_GRID_SIZE, DEFAULT_GRID_SIZE
        )

        # tally total building area if there is defined
        bldg_area_idx = layer_field_index(zone_layer, area_field)
        zone_area = {}
        zone_has_area = False
        if bldg_area_idx > 0:
            zone_has_area = True
            zone_gid_idx = layer_field_index(zone_layer, GID_FIELD_NAME)
            for _f in layer_features(zone_layer):
                gid = _f.attributeMap()[zone_gid_idx].toString()
                area = _f.attributeMap()[bldg_area_idx].toDouble()[0]
                if zone_area.has_key(gid):
                    zone_area[gid] = str(float(zone_area[gid])) + area
                else:
                    zone_area[gid] = area

        # 3. intersect grids and zones to obtain polygons with
        # - grid_id and zone_id
        # - ratio of grid covered by zone (polygon area / zone area)
        # apply ratio to zone building count to obtain count assigned to polygon
        tmp_join = "joined_%s" % get_unique_filename()
        tmp_join_file = "%s%s.shp" % (self._tmp_dir, tmp_join)
        try:
            # do intersection
            analyzer.intersection(tmp_grid_lyr1, zone_layer, tmp_join_file)
            tmp_join_layer = load_shapefile(tmp_join_file, tmp_join)
        except AssertionError as err:
            raise OperatorError(str(err), self.__class__)
        except Exception as err:
            raise OperatorError(str(err), self.__class__)

        # do tally
        zone_gid_idx = layer_field_index(tmp_join_layer, GID_FIELD_NAME)
        grid_gid_idx = layer_field_index(tmp_join_layer, "GRID_GID")
        bldg_cnt_idx = layer_field_index(tmp_join_layer, count_field)
        for _f in layer_features(tmp_join_layer):
            geom = _f.geometry()
            geom.transform(mercator_transform)
            area = geom.area()

            # generate all stats of interest
            zone_gid = _f.attributeMap()[zone_gid_idx].toString()
            grid_gid = _f.attributeMap()[grid_gid_idx].toString()
            stat = zone_stat[zone_gid]
            # calculate count/area as proportion of total zone area
            area_ratio = area / stat[area_idx]
            if bldg_cnt_idx > 0:
                bldg_cnt = _f.attributeMap()[bldg_cnt_idx].toDouble()[0] * area_ratio
            else:
                bldg_cnt = 0
            if zone_has_area:
                area = zone_area[zone_gid] * area_ratio
#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsCoordinateReferenceSystem [as 别名]
# 或者: from qgis.core.QgsCoordinateReferenceSystem import createFromEpsg [as 别名]
    def do_operation(self):
        """ perform footprint load operation """
        
        # input/output data checking already done during property set        
        # load and verify
        infile = self.inputs[0].value
        
        tmp_fp_layername = 'fp_%s' % get_unique_filename()
        tmp_fp_layer = load_shapefile(infile, tmp_fp_layername)
        if not tmp_fp_layer:
            raise OperatorError('Error loading footprint file' % (infile), self.__class__)

        if self._fp_ht_field is not None:
            ht_idx = layer_field_index(tmp_fp_layer, self._fp_ht_field)
        else:
            ht_idx = -1
        logAPICall.log('tmp_fp_layer.crs().epsg() %s ' % tmp_fp_layer.crs().epsg(),
                       logAPICall.DEBUG)
        if tmp_fp_layer.crs().epsg() != self._crs.epsg():
            transform = QgsCoordinateTransform(tmp_fp_layer.crs(), self._crs)
            transform_required = True
        else:
            transform_required = False
        
        mercator_crs = QgsCoordinateReferenceSystem()
        #mercator_crs.createFromProj4("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
        mercator_crs.createFromEpsg(3395)        
        mercator_transform = QgsCoordinateTransform(tmp_fp_layer.crs(), mercator_crs)
        
        # output grid
        fields = {
            0 : QgsField(GID_FIELD_NAME, QVariant.Int),
            1 : QgsField(LON_FIELD_NAME, QVariant.Double),
            2 : QgsField(LAT_FIELD_NAME, QVariant.Double),
            3 : QgsField(AREA_FIELD_NAME, QVariant.Double),
            4 : QgsField(HT_FIELD_NAME, QVariant.Int),
        }
        output_file = '%sfpc_%s.shp' % (self._tmp_dir, get_unique_filename())
        logAPICall.log('create outputfile %s ... ' % output_file, logAPICall.DEBUG)        
        try:
            writer = QgsVectorFileWriter(output_file, "utf-8", fields, QGis.WKBPoint, self._crs, "ESRI Shapefile")
            f = QgsFeature()
            gid = 0
            for _f in layer_features(tmp_fp_layer):
                # NOTE: geom.transform does projection in place to underlying
                #       C object, for some reason, multiple projection does not
                #       work correctly. following is a work-around
                 
                # 1. get geometry
                geom = _f.geometry()
                # 2. get original centroid point and project is required
                centroid  = geom.centroid().asPoint()
                if transform_required:
                    t_centroid = transform.transform(centroid)
                else:
                    t_centroid = centroid
                
                # 3. project into mercator and get area in m2
                geom.transform(mercator_transform)
                area = geom.area()
                
                # write to file
                gid += 1
                f.setGeometry(QgsGeometry.fromPoint(t_centroid))
                f.addAttribute(0, QVariant(gid))
                f.addAttribute(1, QVariant(t_centroid.x()))
                f.addAttribute(2, QVariant(t_centroid.y()))
                f.addAttribute(3, QVariant(area))
                if ht_idx != -1:
                    f.addAttribute(4, _f.attributeMap()[ht_idx])
                else:
                    f.addAttribute(4, QVariant(0))
                writer.addFeature(f)            
            del writer, f
        except Exception as err:
            remove_shapefile(output_file)
            raise OperatorError("error creating footprint centroids: %s" % err, self.__class__)

        fp_layer = load_shapefile(output_file, tmp_fp_layername)
        if not fp_layer:
            raise OperatorError('Error loading footprint centroid file' % (output_file), self.__class__)        
        
        # clean up
        del tmp_fp_layer
        
        # store data in output
        self.outputs[0].value = fp_layer
        self.outputs[1].value = output_file