Python LineString.transform方法代码示例

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
 def clean_location(self):
     p = fromstr(self.cleaned_data['location'])
     line = LineString(p, self.competition.location.coordinates, srid=4326) # lon/lat (srid=4326)
     line.transform(53031) # transform to two-point equidistant projection (srid=53031)
     if line.length > self.competition.check_in_distance:
         raise forms.ValidationError(_("""You are not close enough to 
             enter the competition, or you GPS might not be turned on. In the latter 
             case, turn on your device's GPS and reload the page."""))
     return self.cleaned_data['location']

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
 def check_consistency(self, cleaned_data):
     if cleaned_data['city'].country != cleaned_data['country']:
         return (False, "The location's city and country do not match.")
     if cleaned_data['coordinates']:
         line = LineString(cleaned_data['coordinates'], cleaned_data['city'].coordinates, srid=4326)
         line.transform(53031)  # transform to two-point equidistant project
         # if coordinate is more than 500km away from city
         if line.length > 500000:
             return (False, "The location's coordinates are more than 500km away from its city.")
     return (True, "")

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
def sensors(point_list):
    route = LineString([[point['lng'], point['lat']] for point in point_list], srid=canonical_projection)
    route.transform(google_projection)
    sensor_map = []
    for sensor_ind, sensor in enumerate(sensors_transformed):
      if sensor['compare'].distance(route) < FUZZY_DIST:
          sensor_map.append({
              'loc':sensor['map'],
              'true':true_sensor_value[sensor_ind],
              'predicted':lsqr[sensor_ind]
          })
    return sensor_map

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
class Link:
    def __init__(self, geom, speed):
        self.line = LineString(geom['coordinates'], srid=4326)
        self.line.transform(32611)
        self.speed = speed
        self.cur_dist = 0
        self.total_dist = self.line.length

    def move(self, dt):
        distance = self.speed * dt
        overshot = distance + self.cur_dist - self.total_dist
        if overshot > 0:
            self.cur_dist = self.total_dist
            return float(overshot) / self.speed
        else:
            self.cur_dist += distance
            return 0

    def get_position(self):
        pt = self.line.interpolate(self.cur_dist)
        pt.transform(3857)
        return list(pt)

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
    def filter(self, queryset, containment="overlaps"):
        if not len(self):
            return queryset

        qs = queryset

        bbox = [None, None, None, None]
        srid = self.srid

        if srid is None:
            srid = 4326
        max_extent = crss.crs_bounds(srid)
        tolerance = crss.crs_tolerance(srid)

        for subset in self:
            if isinstance(subset, Slice):
                is_slice = True
                value = subset.value
            elif isinstance(subset, Trim):
                is_slice = False
                low = subset.low
                high = subset.high

            if subset.is_temporal:
                if is_slice:
                    qs = qs.filter(
                        begin_time__lte=value,
                        end_time__gte=value
                    )

                else:
                    if high is not None:
                        qs = qs.filter(
                            begin_time__lte=high
                        )
                    if low is not None:
                        qs = qs.filter(
                            end_time__gte=low
                        )
                
                    # check if the temporal bounds must be strictly contained
                    if containment == "contains":
                        if high is not None:
                            qs = qs.filter(
                                end_time__lte=high
                            )
                        if low is not None:
                            qs = qs.filter(
                                begin_time__gte=low
                            )

            else:
                if is_slice:
                    if subset.is_x:
                        line = LineString(
                            (value, max_extent[1]),
                            (value, max_extent[3])
                        )
                    else:
                        line = LineString(
                            (max_extent[0], value),
                            (max_extent[2], value)
                        )
                    line.srid = srid
                    if srid != 4326:
                        line.transform(4326)
                    qs = qs.filter(footprint__intersects=line)

                else:
                    if subset.is_x:
                        bbox[0] = subset.low
                        bbox[2] = subset.high
                    else:
                        bbox[1] = subset.low
                        bbox[3] = subset.high


        if bbox != [None, None, None, None]:
            bbox = map(
                lambda v: v[0] if v[0] is not None else v[1], 
                zip(bbox, max_extent)
            )

            bbox[0] -= tolerance; bbox[1] -= tolerance
            bbox[2] += tolerance; bbox[3] += tolerance

            logger.debug(
                "Applying BBox %s with containment '%s'." % (bbox, containment)
            )

            poly = Polygon.from_bbox(bbox)
            poly.srid = srid

            if srid != 4326:
                poly.transform(4326)
            if containment == "overlaps":
                qs = qs.filter(footprint__intersects=poly)
            elif containment == "contains":
                qs = qs.filter(footprint__within=poly)

#.........这里部分代码省略.........

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
def render_static(request, height=None, width=None, format='png',
                  background='satellite', bounds=None, center=None, render_srid=3857):

# width and height
    width = int(width)
    height = int(height)
    if width > settings.MAX_IMAGE_DIMENSION or \
        height > settings.MAX_IMAGE_DIMENSION or \
        width <= 1 or height <= 1:
        logging.debug("Invalid size")
        return HttpResponseBadRequest("Invalid image size, both dimensions must be in range %i-%i" % (1, settings.MAX_IMAGE_DIMENSION))

# image format
    if format not in IMAGE_FORMATS:
        logging.error("unknown image format %s" % format)
        return HttpResponseBadRequest("Unknown image format, available formats: " + ", ".join(IMAGE_FORMATS))

    if format.startswith('png'):
        mimetype = 'image/png'
    elif format.startswith('jpeg'):
        mimetype = 'image/jpeg'

# bounds
    bounds_box = None
    if bounds:
        bounds_components = bounds.split(',')
        if len(bounds_components) != 4:
            return HttpResponseBadRequest("Invalid bounds, must be 4 , separated numbers")
        bounds_components = [float(f) for f in bounds_components]

        if not (-180 < bounds_components[0] < 180) or not (-180 < bounds_components[2] < 180):
            logging.error("x out of range %f or %f" % (bounds_components[0], bounds_components[2]))
            return HttpResponseBadRequest("x out of range %f or %f" % (bounds_components[0], bounds_components[2]))
        if not (-90 < bounds_components[1] < 90) or not (-90 < bounds_components[3] < 90):
            logging.error("y out of range %f or %f" % (bounds_components[1], bounds_components[3]))
            return HttpResponseBadRequest("y out of range %f or %f" % (bounds_components[1], bounds_components[3]))

        ll = Point(bounds_components[0], bounds_components[1], srid=4326)
        ll.transform(render_srid)

        ur = Point(bounds_components[2], bounds_components[3], srid=4326)
        ur.transform(render_srid)
        bounds_box = mapnik.Box2d(ll.x, ll.y, ur.x, ur.y)
    elif center:
        center_components = center.split(',')
        if len(center_components) != 3:
            return HttpResponseBadRequest()
        lon = float(center_components[0])
        lat = float(center_components[1])
        zoom = int(center_components[2])
        # todo calc bounds from center and zoom

# baselayer
    if background not in settings.BASE_LAYERS and background != 'none':
        return HttpResponseNotFound("Background not found")

# GeoJSON post data
    if request.method == "POST" and len(request.body):
        input_data = json.loads(request.body)
    else:
        input_data = None

    if not bounds and not center and not input_data:
        return HttpResponseBadRequest("Bounds, center, or post data is required.")

# initialize map
    m = mapnik.Map(width, height)
    m.srs = '+init=epsg:' + str(render_srid)

# add a tile source as a background
    if background != "none":
        background_file = settings.BASE_LAYERS[background]
        background_style = mapnik.Style()
        background_rule = mapnik.Rule()
        background_rule.symbols.append(mapnik.RasterSymbolizer())
        background_style.rules.append(background_rule)
        m.append_style('background style', background_style)
        tile_layer = mapnik.Layer('background')
        tile_layer.srs = '+init=epsg:' + str(render_srid)
        tile_layer.datasource = mapnik.Gdal(base=settings.BASE_LAYER_DIR, file=background_file)
        tile_layer.styles.append('background style')
        m.layers.append(tile_layer)

# add features from geojson
    if input_data and input_data['type'] == "Feature":
        features = [input_data]
    elif input_data and input_data['type'] == "FeatureCollection":
        if 'features' not in input_data:
            return HttpResponseBadRequest()
        features = input_data['features']
    else:
        features = []

    logging.debug("Adding %d features to map" % len(features))

    geometries = []
    point_features = []
    fid = 0
    for feature in features:
        if 'geometry' not in feature:
#.........这里部分代码省略.........

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
class Path(AddPropertyMixin, MapEntityMixin, AltimetryMixin,
           TimeStampedModelMixin, StructureRelated):
    geom = models.LineStringField(srid=settings.SRID, spatial_index=False)
    geom_cadastre = models.LineStringField(null=True, srid=settings.SRID, spatial_index=False,
                                           editable=False)
    valid = models.BooleanField(db_column='valide', default=True, verbose_name=_(u"Validity"),
                                help_text=_(u"Approved by manager"))
    visible = models.BooleanField(db_column='visible', default=True, verbose_name=_(u"Visible"),
                                  help_text=_(u"Shown in lists and maps"))
    name = models.CharField(null=True, blank=True, max_length=20, db_column='nom', verbose_name=_(u"Name"),
                            help_text=_(u"Official name"))
    comments = models.TextField(null=True, blank=True, db_column='remarques', verbose_name=_(u"Comments"),
                                help_text=_(u"Remarks"))

    departure = models.CharField(null=True, blank=True, default="", max_length=250, db_column='depart', verbose_name=_(u"Departure"),
                                 help_text=_(u"Departure place"))
    arrival = models.CharField(null=True, blank=True, default="", max_length=250, db_column='arrivee', verbose_name=_(u"Arrival"),
                               help_text=_(u"Arrival place"))

    comfort = models.ForeignKey('Comfort',
                                null=True, blank=True, related_name='paths',
                                verbose_name=_("Comfort"), db_column='confort')
    source = models.ForeignKey('PathSource',
                               null=True, blank=True, related_name='paths',
                               verbose_name=_("Source"), db_column='source')
    stake = models.ForeignKey('Stake',
                              null=True, blank=True, related_name='paths',
                              verbose_name=_("Maintenance stake"), db_column='enjeu')
    usages = models.ManyToManyField('Usage',
                                    blank=True, related_name="paths",
                                    verbose_name=_(u"Usages"), db_table="l_r_troncon_usage")
    networks = models.ManyToManyField('Network',
                                      blank=True, related_name="paths",
                                      verbose_name=_(u"Networks"), db_table="l_r_troncon_reseau")
    eid = models.CharField(verbose_name=_(u"External id"), max_length=1024, blank=True, null=True, db_column='id_externe')
    draft = models.BooleanField(db_column='brouillon', default=False, verbose_name=_(u"Draft"), db_index=True)

    objects = PathManager()
    include_invisible = PathInvisibleManager()

    is_reversed = False

    @property
    def length_2d(self):
        if self.geom:
            return self.geom.length
        else:
            return None

    @classproperty
    def length_2d_verbose_name(cls):
        return _(u"2D Length")

    @property
    def length_2d_display(self):
        return round(self.length_2d, 1)

    def kml(self):
        """ Exports path into KML format, add geometry as linestring """
        kml = simplekml.Kml()
        geom3d = self.geom_3d.transform(4326, clone=True)  # KML uses WGS84
        line = kml.newlinestring(name=self.name,
                                 description=plain_text(self.comments),
                                 coords=geom3d.coords)
        line.style.linestyle.color = simplekml.Color.red  # Red
        line.style.linestyle.width = 4  # pixels
        return kml.kml()

    def __unicode__(self):
        return self.name or _('path %d') % self.pk

    class Meta:
        db_table = 'l_t_troncon'
        verbose_name = _(u"Path")
        verbose_name_plural = _(u"Paths")
        permissions = MapEntityMixin._meta.permissions + [("add_draft_path", "Can add draft Path"),
                                                          ("change_draft_path", "Can change draft Path"),
                                                          ("delete_draft_path", "Can delete draft Path"),
                                                          ]

    @classmethod
    def closest(cls, point, exclude=None):
        """
        Returns the closest path of the point.
        Will fail if no path in database.
        """
        # TODO: move to custom manager
        if point.srid != settings.SRID:
            point = point.transform(settings.SRID, clone=True)
        qs = cls.objects.exclude(draft=True)
        if exclude:
            qs = qs.exclude(pk=exclude.pk)
        return qs.exclude(visible=False).distance(point).order_by('distance')[0]

    def is_overlap(self):
        return not PathHelper.disjoint(self.geom, self.pk)

    def reverse(self):
        """
        Reverse the geometry.
#.........这里部分代码省略.........

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
def find_shared_segments(R, S):
  """Return a list of LineStrings where two routes overlap

  Keyword arguments:
  route_a -- LineString route
  route_b -- LineString route
  This method should be commutative.
  
  Algorithmic idea:
  For a given set of routes that are not identical:
    - Go through each route (R and S), and make a list, P(*) for each that contains points not on the other
      route.
    - Starting with P(R), iterate through elements p_i = P(R)_i:
      - Find p_i in original polyline, R.
      - Find the vertex of the item before p_i in R, v_s. v_s should be on the other route S.
      - Starting from this element in R, move forward through route until you find the first vertex
        that is on S, call this v_f. Remove elements not on S from P(R), including p_i
      - Starting from v_s, find closest vertex on S, call it g. This could be before or after where v_s
        intersects the route.
      - Move along S in both directions from g, taking ordered pairs where order is index in S, and identify
        the first pair where the first element is on R, but the second element is not.
        - If the first element of this tuple is after g, then that element is end of a shared leg.
        - Otherwise, v_s is the end of that shared leg.
      - Starting from v_f, find closest vertex on S, call it g, this could be before or after where v_f
        intersects the route.
      - Move along S in both directions from g, taking ordered pairs where order is index in S, and identify
        the first pair where the first element is not on R, but the second element is.
        - If the second element of this tuple is before g, then that element is the start of a shared leg.
        - Otherwise, v_f is the start of that shared leg.
    - At this point, we have a list of starting points and ending points of shared legs, combining them
      in to a polyline from the two routes seems tractable. TODO: figure this out when brain is less fried.
  """
  R_points_not_on_S = []
  S_points_not_on_R = []
  S_trans = S.clone()
  S_trans.set_srid(4326)
  S_trans.transform(900913)
  for pt in R:
      place = Point(pt)
      place.set_srid(4326)
      place.transform(900913)
      if place.distance(S_trans) > FUZZY_DIST:
          R_points_not_on_S.append(pt)
  #TODO: refactor these into single function call
  S_trans = S.clone()
  S_trans.set_srid(4326)
  S_trans.transform(900913)
  R_trans = R.clone()
  R_trans.set_srid(4326)
  R_trans.transform(900913)
  for pt in S:
      place = Point(pt, srid=4326)
      place.transform(900913)
      if place.distance(R_trans) > FUZZY_DIST:
          S_points_not_on_R.append(pt)
  # we know they start at the same point
  shared_leg_list = []
  shared_leg = []
  shared_leg_start_index = 0
  n = len(R_points_not_on_S)
  while n > 0:
      p_i = R_points_not_on_S[0]
      j = R.index(p_i)
      v_s = R[j - 1]
      f = j
      v_f = p_i
      while v_f in R_points_not_on_S:
          idx = R_points_not_on_S.index(v_f)
          del R_points_not_on_S[idx]
          n = n - 1
          f = f + 1
          v_f = R[f]
      # We know v_f is fuzzy-on S,  so we can iterate through pairs of S and find the segment it is fuzzy-on
      before_index = 0
      for i, start_vertex in enumerate(S):
          if i == len(S):
              break
          end_vertex = S[i + 1]
          line = LineString([start_vertex, end_vertex])
          line.set_srid(4326)
          line.transform(900913)
          pt = Point(v_s)
          pt.set_srid(4326)
          pt.transform(900913)
          if pt.distance(line) < FUZZY_DIST:
              before_index = i
              break
      # At this point, we know shared_leg_start_index..before_index is certainly on the path.
      shared_leg = S[shared_leg_start_index:(before_index+1)]
      shared_leg.append(v_s)
      after_index = before_index + 1
      # Either v_s is the end of the previous shared segment, or after_index or something following that,
      # so go until you find a point on S not on R, starting at after_index
      pt = Point(S[after_index], srid=4326)
      pt.transform(900913)
      while pt.distance(R_trans) < FUZZY_DIST:
          shared_leg.append(S[after_index])
          after_index = after_index + 1
      # should check that shared_leg is not just first element. In fact, TODO: go back an check what happens
      # if the first element is the only shared element.
#.........这里部分代码省略.........

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
    def matches(self, eo_object, containment="overlaps"):
        """ Check if the given :class:`EOObject
        <eoxserver.resources.coverages.models.EOObject>` matches the given
        subsets.

        :param eo_object: the ``EOObject`` to match
        :param containment: either "overlaps" or "contains"
        :returns: a boolean value indicating if the object is contained in the
                  given subsets
        """
        if not len(self):
            return True

        bbox = [None, None, None, None]
        srid = self.srid
        if srid is None:
            srid = 4326
        max_extent = crss.crs_bounds(srid)
        tolerance = crss.crs_tolerance(srid)

        # check if time intervals are configured as "open" or "closed"
        config = get_eoxserver_config()
        reader = SubsetConfigReader(config)
        # note that the operator is inverted from filter() above as the
        # filters use an inclusive search whereas here it's exclusive
        if reader.time_interval_interpretation == "closed":
            gt_op = operator.gt
            lt_op = operator.lt
        else:
            gt_op = operator.ge
            lt_op = operator.le

        footprint = eo_object.footprint
        begin_time = eo_object.begin_time
        end_time = eo_object.end_time

        for subset in self:
            if isinstance(subset, Slice):
                is_slice = True
                value = subset.value
            elif isinstance(subset, Trim):
                is_slice = False
                low = subset.low
                high = subset.high
                # we need the value in case low == high
                value = low

            if subset.is_temporal:
                if is_slice or (low == high and containment == "overlaps"):
                    if begin_time > value or end_time < value:
                        return False
                elif low == high:
                    if begin_time < value or end_time > value:
                        return False
                else:
                    # check if the temporal bounds must be strictly contained
                    if containment == "contains":
                        if high is not None:
                            if gt_op(end_time, high):
                                return False
                        if low is not None:
                            if lt_op(begin_time, low):
                                return False
                    # or just overlapping
                    else:
                        if high is not None:
                            if gt_op(begin_time, high):
                                return False
                        if low is not None:
                            if lt_op(end_time, low):
                                return False

            else:
                if is_slice:
                    if subset.is_x:
                        line = LineString(
                            (value, max_extent[1]),
                            (value, max_extent[3])
                        )
                    else:
                        line = LineString(
                            (max_extent[0], value),
                            (max_extent[2], value)
                        )
                    line.srid = srid
                    if srid != 4326:
                        line.transform(4326)

                    if not line.intersects(footprint):
                        return False

                else:
                    if subset.is_x:
                        bbox[0] = subset.low
                        bbox[2] = subset.high
                    else:
                        bbox[1] = subset.low
                        bbox[3] = subset.high

        if bbox != [None, None, None, None]:
#.........这里部分代码省略.........

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
    def filter(self, queryset, containment="overlaps"):
        """ Filter a :class:`Django QuerySet <django.db.models.query.QuerySet>`
        of objects inheriting from :class:`EOObject
        <eoxserver.resources.coverages.models.EOObject>`.

        :param queryset: the ``QuerySet`` to filter
        :param containment: either "overlaps" or "contains"
        :returns: a ``QuerySet`` with additional filters applied
        """
        if not len(self):
            return queryset

        qs = queryset

        bbox = [None, None, None, None]
        srid = self.srid

        if srid is None:
            srid = 4326
        max_extent = crss.crs_bounds(srid)
        tolerance = crss.crs_tolerance(srid)

        # check if time intervals are configured as "open" or "closed"
        config = get_eoxserver_config()
        reader = SubsetConfigReader(config)
        if reader.time_interval_interpretation == "closed":
            gt_op = "__gte"
            lt_op = "__lte"
        else:
            gt_op = "__gt"
            lt_op = "__lt"

        for subset in self:
            if isinstance(subset, Slice):
                is_slice = True
                value = subset.value
            elif isinstance(subset, Trim):
                is_slice = False
                low = subset.low
                high = subset.high
                # we need the value in case low == high
                value = low

            if subset.is_temporal:
                if is_slice or (high == low and containment == "overlaps"):
                    qs = qs.filter(
                        begin_time__lte=value,
                        end_time__gte=value
                    )

                elif high == low:
                    qs = qs.filter(
                        begin_time__gte=value,
                        end_time__lte=value
                    )

                else:
                    # check if the temporal bounds must be strictly contained
                    if containment == "contains":
                        if high is not None:
                            qs = qs.filter(**{
                                "end_time" + lt_op: high
                            })
                        if low is not None:
                            qs = qs.filter(**{
                                "begin_time" + gt_op: low
                            })
                    # or just overlapping
                    else:
                        if high is not None:
                            qs = qs.filter(**{
                                "begin_time" + lt_op: high
                            })
                        if low is not None:
                            qs = qs.filter(**{
                                "end_time" + gt_op: low
                            })

            else:
                if is_slice:
                    if subset.is_x:
                        line = LineString(
                            (value, max_extent[1]),
                            (value, max_extent[3])
                        )
                    else:
                        line = LineString(
                            (max_extent[0], value),
                            (max_extent[2], value)
                        )
                    line.srid = srid
                    if srid != 4326:
                        line.transform(4326)
                    qs = qs.filter(footprint__intersects=line)

                else:
                    if subset.is_x:
                        bbox[0] = subset.low
                        bbox[2] = subset.high
                    else:
#.........这里部分代码省略.........

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
    def matches(self, eo_object, containment="overlaps"):
        if not len(self):
            return True

        bbox = [None, None, None, None]
        srid = self.srid
        if srid is None:
            srid = 4326
        max_extent = crss.crs_bounds(srid)
        tolerance = crss.crs_tolerance(srid)

        footprint = eo_object.footprint
        begin_time = eo_object.begin_time
        end_time = eo_object.end_time

        for subset in self:
            if isinstance(subset, Slice):
                is_slice = True
                value = subset.value
            elif isinstance(subset, Trim):
                is_slice = False
                low = subset.low
                high = subset.high

            if subset.is_temporal:
                if is_slice:
                    if begin_time > value or end_time < value:
                        return False
                elif low is None and high is not None:
                    if begin_time > high:
                        return False
                elif low is not None and high is None:
                    if end_time < low:
                        return False
                else:
                    if begin_time > high or end_time < low:
                        return False

            else:
                if is_slice:
                    if subset.is_x:
                        line = LineString(
                            (value, max_extent[1]),
                            (value, max_extent[3])
                        )
                    else:
                        line = LineString(
                            (max_extent[0], value),
                            (max_extent[2], value)
                        )
                    line.srid = srid
                    if srid != 4326:
                        line.transform(4326)

                    if not line.intersects(footprint):
                        return False

                else:
                    if subset.is_x:
                        bbox[0] = subset.low
                        bbox[2] = subset.high
                    else:
                        bbox[1] = subset.low
                        bbox[3] = subset.high

        if bbox != [None, None, None, None]:
            bbox = map(
                lambda v: v[0] if v[0] is not None else v[1],
                zip(bbox, max_extent)
            )

            bbox[0] -= tolerance
            bbox[1] -= tolerance
            bbox[2] += tolerance
            bbox[3] += tolerance

            logger.debug(
                "Applying BBox %s with containment '%s'." % (bbox, containment)
            )

            poly = Polygon.from_bbox(bbox)
            poly.srid = srid

            if srid != 4326:
                poly.transform(4326)
            if containment == "overlaps":
                if not footprint.intersects(poly):
                    return False
            elif containment == "contains":
                if not footprint.within(poly):
                    return False
        return True

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
 def line(self):
     geom = LineString([(x, y) for x, y, _, _ in self.gcps])
     geom.srid = self._srid
     if self._srid != 4326:
         geom.transform(4326)
     return geom

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
def save_into_database(rib_path, rmb_path, putdict, sewerdict, rmberrors):
    # Get sewerage name, try to create sewerage
    # If it exists, return with an error
    sewerage_name = os.path.basename(rmb_path)[:-4]  # Minus ".RMB"

    if models.Sewerage.objects.filter(name=sewerage_name).exists():
        rmberrors.append(Error(
                line_number=0,
                message=("Er bestaat al een stelsel met de naam {name}. "
                 "Verwijder het op de archiefpagina, of gebruik "
                 "een andere naam.").format(name=sewerage_name)))
        return

    # Files are copied only at the end
    sewerage = models.Sewerage.objects.create(
        name=sewerage_name,
        rib=None,  # Filled in later
        rmb=None,
        active=True)

    # Save the puts, keep a dictionary
    saved_puts = dict()
    for put_id, putinfo in putdict.items():
        saved_puts[put_id] = models.Manhole.objects.create(
            sewerage=sewerage,
            code=put_id,
            sink=int(putinfo['is_sink']),
            ground_level=putinfo['surface_level'],
            the_geom=Point(*putinfo['coordinate']))

    # Save the sewers, use the dictionary
    saved_sewers = dict()
    for sewer_id, sewerinfo in sewerdict.items():
        manhole1 = saved_puts[sewerinfo['manhole_code_1']]
        manhole2 = saved_puts[sewerinfo['manhole_code_2']]
        sewer_line_rd = LineString(manhole1.the_geom, manhole2.the_geom)
        sewer_line_rd.set_srid(4326)
        sewer_line_rd.transform(RD)

        saved_sewers[sewer_id] = models.Sewer.objects.create(
            sewerage=sewerage,
            code=sewer_id,
            quality=models.Sewer.QUALITY_UNKNOWN,
            diameter=sewerinfo['diameter'],
            manhole1=manhole1,
            manhole2=manhole2,
            bob1=sewerinfo['bob_1'],
            bob2=sewerinfo['bob_2'],
            the_geom=LineString(manhole1.the_geom, manhole2.the_geom),
            the_geom_length=sewer_line_rd.length)

    # Save the measurements
    sewer_measurements_dict = dict()
    for sewer_id, sewerinfo in sewerdict.items():
        measurements = sewerinfo['measurements']
        sewer = saved_sewers[sewer_id]

        if measurements:
            sewer_measurements = [
                # Create the SewerMeasurement objects, but don't save
                # them yet!
                models.SewerMeasurement(
                    sewer=sewer,
                    dist=m['dist'],
                    virtual=False,
                    water_level=None,
                    flooded_pct=None,
                    bob=m['bob'],
                    obb=m['bob'] + sewerinfo['diameter'],
                    the_geom=Point(*m['coordinate']))
                for m in measurements]

            # Quality
            sewer.judge_quality(sewer_measurements)
            sewer.save()

            # BOB correction ("sawtooth" phenomenon)
            correct_bob_values(sewer, sewer_measurements)

            # Create two virtual sewer measurements for the start and
            # end of the sewer
            virtual_start = models.SewerMeasurement(
                sewer=sewer, dist=0, virtual=True, water_level=None,
                flooded_pct=None, bob=sewer.bob1,
                obb=sewer.bob1 + sewerinfo['diameter'],
                the_geom=sewer.manhole1.the_geom)
            virtual_end = models.SewerMeasurement(
                sewer=sewer, dist=sewer.the_geom_length,
                virtual=True, water_level=None,
                flooded_pct=None, bob=sewer.bob2,
                obb=sewer.bob2 + sewerinfo['diameter'],
                the_geom=sewer.manhole2.the_geom)

            # Note: we MUST add those two virtual points only after
            # doing the sawtooth correction, otherwise the sawtooth
            # correction will think that everything is fine already
            # since the first and end points would be equal to the
            # bobs of the sewer...
            sewer_measurements = (
                [virtual_start] + sewer_measurements + [virtual_end])
#.........这里部分代码省略.........

# 需要导入模块: from django.contrib.gis.geos import LineString [as 别名]
# 或者: from django.contrib.gis.geos.LineString import transform [as 别名]
    def handle(self, *args, **options):
        
        # Parse arguments.
        if len(args) < 2 or len(args) > 3:
            raise CommandError("Missing arguments. Try running with --help.")

        shapefile_base = args[0]
        agency_name = args[1]
        srid = int(args[2]) if len(args) > 2 else None

        # Verify the agency name.
        try:
            agency = Agency.objects.get(agency_name=agency_name)
        except Agency.DoesNotExist:
            self.stderr.write("No agency found for '%s'.\n" % agency_name)
            if Agency.objects.all():
                self.stderr.write("Agency name must be one of:\n")
                for agency in Agency.objects.all():
                    self.stderr.write("- %s\n" % agency.agency_name)
            return
    
        # Read in the shapefile.
        self.stdout.write("Reading shapefile from '%s'.\n" % shapefile_base)
        sf = shapefile.Reader(shapefile_base)

        # Verify that there is a NAME field in this shapefile.
        # Ignore first field ('DeletionField' added by shapefile library).
        name_index = [f[0] for f in sf.fields[1:]].index('NAME')

        if not srid:
            # Try opening a .prj file to get the shapefile projection.
            prj_file = open(shapefile_base + '.prj')
            prj_text = prj_file.read()

            # Connect to API to get SRID from projection, if not provided.
            # Inspired by http://gis.stackexchange.com/a/7784
            self.stdout.write("Querying prj2epsg to find SRID.\n")
            url = 'http://prj2epsg.org/search.json'
            query = urlencode(
                {'exact': True,
                 'error': True,
                 'mode': 'wkt',
                 'terms': prj_text})
            response = urlopen(url, query)
            data = json.loads(response.read())
            srid = int(data['codes'][0]['code'])
            self.stdout.write("Using SRID %d.\n" % srid)

        for shaperecord in sf.shapeRecords():
            shape = shaperecord.shape
            record = shaperecord.record

            # Get name of the record.
            name = record[name_index]
            self.stdout.write("Finding patterns for shape '%s'.\n" % name)

            # Verify that this shape is a PolyLine (type 3).
            if 3 != shape.shapeType:
                pass

            # Create the LineString with the appropriate SRID
            ls = LineString([list(pt) for pt in shape.points], srid=srid)

            # Transform into a common SRID (4326).
            ls.transform(4326)

            # Find the routes and patterns that are covered by this shape.
            # Consider all patterns where both the start and end stops are
            # within a certain distance to the shape.
            distance_threshold = 1609

            cursor = connection.cursor()

            # Find matching pattern ID, shape ID, and first and last stop locations.
            query = 'SELECT a.pattern_id, gp.shape_id, first_stop_location, last_stop_location FROM (SELECT pattern_id, gs.stop_id AS first_stop_id, location AS first_stop_location, stop_name AS first_stop_name FROM gtfs_patternstop gps INNER JOIN gtfs_stop gs ON gps.stop_id = gs.id WHERE gps.order = 0) a INNER JOIN (SELECT gps.pattern_id, gs.stop_id AS last_stop_id, location AS last_stop_location, stop_name AS last_stop_name FROM gtfs_patternstop gps INNER JOIN (SELECT pattern_id, MAX(gps.order) max_order FROM gtfs_patternstop gps GROUP BY pattern_id) gps2 ON gps.pattern_id = gps2.pattern_id AND gps.order = gps2.max_order INNER JOIN gtfs_stop gs ON gps.stop_id = gs.id) b ON a.pattern_id = b.pattern_id INNER JOIN gtfs_pattern gp ON a.pattern_id = gp.id INNER JOIN gtfs_route gr ON gp.route_id = gr.id WHERE gr.agency_id = %s AND ST_Distance_Sphere(first_stop_location, ST_Line_Interpolate_Point(ST_GeomFromEWKB(%s), ST_Line_Locate_Point(ST_GeomFromEWKB(%s), first_stop_location))) < %s AND ST_Distance_Sphere(last_stop_location, ST_Line_Interpolate_Point(ST_GeomFromEWKB(%s), ST_Line_Locate_Point(ST_GeomFromEWKB(%s), last_stop_location))) < %s AND (gr.route_short_name ILIKE %s OR gr.route_long_name ILIKE %s) AND gp.shape_id IS NOT NULL'
            args = [agency.id, ls.ewkb, ls.ewkb, distance_threshold, ls.ewkb, ls.ewkb, distance_threshold, name, name]

            cursor.execute(query, args)
            data = cursor.fetchall()

            for row in data:

                pattern_id = row[0]
                shape_id = row[1]
                first_stop_location = row[2]
                last_stop_location = row[3]
            
                # Cut the shape between the first and last stops.
                query = 'SELECT ST_Line_Substring(ST_GeomFromEWKB(%s), LEAST(ST_Line_Locate_Point(%s, %s), ST_Line_Locate_Point(%s, %s)), GREATEST(ST_Line_Locate_Point(%s, %s), ST_Line_Locate_Point(%s, %s)))'
                args = [ls.ewkb, ls.ewkb, first_stop_location, ls.ewkb, last_stop_location, ls.ewkb, first_stop_location, ls.ewkb, last_stop_location]
                cursor.execute(query, args)
                substring = cursor.fetchall()[0]

                # Update the database with the new partial shape.
                query = 'UPDATE gtfs_shape SET line=%s WHERE id=%s'
                args = [substring, shape_id]

                cursor.execute(query, args)
                transaction.commit_unless_managed()