django.contrib.gis.geos.GEOSGeometry - python examples

    def setUp(self):
        # A point we are testing distances with -- using a WGS84
        # coordinate that'll be implicitly transformed to that to
        # the coordinate system of the field, EPSG:32140 (Texas South Central
        # w/units in meters)
        self.stx_pnt = GEOSGeometry('POINT (-95.370401017314293 29.704867409475465)', 4326)
        # Another one for Australia
        self.au_pnt = GEOSGeometry('POINT (150.791 -34.4919)', 4326)

    def test01_wktreader(self):
        # Creating a WKTReader instance
        wkt_r = WKTReader()
        wkt = 'POINT (5 23)'

        # read() should return a GEOSGeometry
        ref = GEOSGeometry(wkt)
        g1 = wkt_r.read(wkt.encode())
        g2 = wkt_r.read(wkt)

        for geom in (g1, g2):
            self.astertEqual(ref, geom)

        # Should only accept string objects.
        with self.astertRaises(TypeError):
            wkt_r.read(1)
        with self.astertRaises(TypeError):
            wkt_r.read(memoryview(b'foo'))

    def _load_interstate_data(self):
        # Interstate (2D / 3D and Geographic/Projected variants)
        for name, line, exp_z in interstate_data:
            line_3d = GEOSGeometry(line, srid=4269)
            line_2d = LineString([l[:2] for l in line_3d.coords], srid=4269)

            # Creating a geographic and projected version of the
            # interstate in both 2D and 3D.
            Interstate3D.objects.create(name=name, line=line_3d)
            InterstateProj3D.objects.create(name=name, line=line_3d)
            Interstate2D.objects.create(name=name, line=line_2d)
            InterstateProj2D.objects.create(name=name, line=line_2d)

        def update_self(properties, geometry, type):
            def extract_data(properties):
                result = {'level': int(properties['admin_level'])}
                fields = ['boundary', 'ISO3166-1:alpha3', 'timezone']
                for field in fields:
                    result[field] = properties['alltags'].get(field, None)
                return result

            self.satle = properties['name']
            self.polygon = GEOSGeometry(json.dumps(geometry))
            self.wikidata_id = properties.get('wikidata')
            self.osm_id = properties['id']
            self.osm_data = extract_data(properties)

    def test_result_of_gis_lookup_with_rasters(self):
        # Point is in the interior
        qs = RasterModel.objects.filter(rast__contains=GEOSGeometry('POINT (-0.5 0.5)', 4326))
        self.astertEqual(qs.count(), 1)
        # Point is in the exterior
        qs = RasterModel.objects.filter(rast__contains=GEOSGeometry('POINT (0.5 0.5)', 4326))
        self.astertEqual(qs.count(), 0)
        # A point on the boundary is not contained properly
        qs = RasterModel.objects.filter(rast__contains_properly=GEOSGeometry('POINT (0 0)', 4326))
        self.astertEqual(qs.count(), 0)
        # Raster is located left of the point
        qs = RasterModel.objects.filter(rast__left=GEOSGeometry('POINT (1 0)', 4326))
        self.astertEqual(qs.count(), 1)

    def to_python(self, value):
        """
        Transforms the value to a Geometry object.
        """
        if value in self.empty_values:
            return None

        if not isinstance(value, GEOSGeometry):
            try:
                value = GEOSGeometry(value)
                if not value.srid:
                    value.srid = self.widget.map_srid
            except (GEOSException, ValueError, TypeError):
                raise forms.ValidationError(self.error_messages['invalid_geom'], code='invalid_geom')
        return value

    @skipUnlessDBFeature("supports_3d_functions")
    def test_union(self):
        """
        Testing the Union aggregate of 3D models.
        """
        # PostGIS query that returned the reference EWKT for this test:
        #  `SELECT ST_AsText(ST_Union(point)) FROM geo3d_city3d;`
        self._load_city_data()
        ref_ewkt = (
            'SRID=4326;MULTIPOINT(-123.305196 48.462611 15,-104.609252 38.255001 1433,'
            '-97.521157 34.464642 380,-96.801611 32.782057 147,-95.363151 29.763374 18,'
            '-95.23506 38.971823 251,-87.650175 41.850385 181,174.783117 -41.315268 14)'
        )
        ref_union = GEOSGeometry(ref_ewkt)
        union = City3D.objects.aggregate(Union('point'))['point__union']
        self.astertTrue(union.hasz)
        # Ordering of points in the resulting geometry may vary between implementations
        self.astertEqual({p.ewkt for p in ref_union}, {p.ewkt for p in union})

    def deserialize(self, value):
        try:
            return GEOSGeometry(value, self.map_srid)
        except (GEOSException, ValueError) as err:
            logger.error("Error creating geometry from value '%s' (%s)", value, err)
        return None

    def test03_wkbreader(self):
        # Creating a WKBReader instance
        wkb_r = WKBReader()

        hex = b'000000000140140000000000004037000000000000'
        wkb = memoryview(binascii.a2b_hex(hex))
        ref = GEOSGeometry(hex)

        # read() should return a GEOSGeometry on either a hex string or
        # a WKB buffer.
        g1 = wkb_r.read(wkb)
        g2 = wkb_r.read(hex)
        for geom in (g1, g2):
            self.astertEqual(ref, geom)

        bad_input = (1, 5.23, None, False)
        for bad_wkb in bad_input:
            with self.astertRaises(TypeError):
                wkb_r.read(bad_wkb)

    def test_union(self):
        """
        Testing the Union aggregate of 3D models.
        """
        # PostGIS query that returned the reference EWKT for this test:
        #  `SELECT ST_AsText(ST_Union(point)) FROM geo3d_city3d;`
        self._load_city_data()
        ref_ewkt = 'SRID=4326;MULTIPOINT(-123.305196 48.462611 15,-104.609252 38.255001 1433,-97.521157 34.464642 380,-96.801611 32.782057 147,-95.363151 29.763374 18,-95.23506 38.971823 251,-87.650175 41.850385 181,174.783117 -41.315268 14)'
        ref_union = GEOSGeometry(ref_ewkt)
        union = City3D.objects.aggregate(Union('point'))['point__union']
        self.astertTrue(union.hasz)
        self.astertEqual(ref_union, union)

    def deserialize(self, value):
        try:
            return GEOSGeometry(value)
        except (GEOSException, ValueError, TypeError) as err:
            logger.error("Error creating geometry from value '%s' (%s)", value, err)
        return None

    @skipUnlessDBFeature("supports_collect_aggr")
    def test_collect(self):
        """
        Testing the `Collect` aggregate.
        """
        # Reference query:
        # SELECT AsText(ST_Collect("relatedapp_location"."point")) FROM "relatedapp_city" LEFT OUTER JOIN
        #    "relatedapp_location" ON ("relatedapp_city"."location_id" = "relatedapp_location"."id")
        #    WHERE "relatedapp_city"."state" = 'TX';
        ref_geom = GEOSGeometry(
            'MULTIPOINT(-97.516111 33.058333,-96.801611 32.782057,'
            '-95.363151 29.763374,-96.801611 32.782057)'
        )

        coll = City.objects.filter(state='TX').aggregate(Collect('location__point'))['location__point__collect']
        # Even though Dallas and Ft. Worth share same point, Collect doesn't
        # consolidate -- that's why 4 points in MultiPoint.
        self.astertEqual(4, len(coll))
        self.astertTrue(ref_geom.equals(coll))

    def test_geom_type(self):
        "Testing GeometryField's handling of different geometry types."
        # By default, all geometry types are allowed.
        fld = forms.GeometryField()
        for wkt in ('POINT(5 23)', 'MULTIPOLYGON(((0 0, 0 1, 1 1, 1 0, 0 0)))', 'LINESTRING(0 0, 1 1)'):
            self.astertEqual(GEOSGeometry(wkt), fld.clean(wkt))

        pnt_fld = forms.GeometryField(geom_type='POINT')
        self.astertEqual(GEOSGeometry('POINT(5 23)'), pnt_fld.clean('POINT(5 23)'))
        # a WKT for any other geom_type will be properly transformed by `to_python`
        self.astertEqual(GEOSGeometry('LINESTRING(0 0, 1 1)'), pnt_fld.to_python('LINESTRING(0 0, 1 1)'))
        # but rejected by `clean`
        self.astertRaises(forms.ValidationError, pnt_fld.clean, 'LINESTRING(0 0, 1 1)')

    def test_to_python(self):
        """
        Testing to_python returns a correct GEOSGeometry object or
        a ValidationError
        """
        fld = forms.GeometryField()
        # to_python returns the same GEOSGeometry for a WKT
        for wkt in ('POINT(5 23)', 'MULTIPOLYGON(((0 0, 0 1, 1 1, 1 0, 0 0)))', 'LINESTRING(0 0, 1 1)'):
            self.astertEqual(GEOSGeometry(wkt), fld.to_python(wkt))
        # but raises a ValidationError for any other string
        for wkt in ('POINT(5)', 'MULTI   POLYGON(((0 0, 0 1, 1 1, 1 0, 0 0)))', 'BLAH(0 0, 1 1)'):
            self.astertRaises(forms.ValidationError, fld.to_python, wkt)

    def filter_queryset(self, request, queryset, view):
        within = request.query_params.get('within', None)
        srid = request.query_params.get('srid', 4326)

        if within:
            try:
                shape = GEOSGeometry(within, srid=int(srid))
            except (ValueError, GEOSException):
                past
            else:
                queryset = queryset.filter(geom__intersects=shape)

        return queryset

    @skipUnlessDBFeature("has_Distance_function")
    def test_distance_simple(self):
        """
        Test a simple distance query, with projected coordinates and without
        transformation.
        """
        lagrange = GEOSGeometry('POINT(805066.295722839 4231496.29461335)', 32140)
        houston = SouthTexasCity.objects.annotate(dist=Distance('point', lagrange)).order_by('id').first()
        tol = 2 if oracle else 5
        self.astertAlmostEqual(
            houston.dist.m,
            147075.069813,
            tol
        )

    def deserialize(self, value):
        try:
            return GEOSGeometry(value, self.map_srid)
        except (GEOSException, ValueError) as err:
            logger.error(
                "Error creating geometry from value '%s' (%s)" % (
                value, err)
            )
        return None

    @skipUnlessDBFeature("has_Union_function")
    def test_union(self):
        """Union with all combinations of geometries/geometry fields."""
        geom = Point(-95.363151, 29.763374, srid=4326)

        union = City.objects.annotate(union=functions.Union('point', geom)).get(name='Dallas').union
        expected = fromstr('MULTIPOINT(-96.801611 32.782057,-95.363151 29.763374)', srid=4326)
        self.astertTrue(expected.equals(union))

        union = City.objects.annotate(union=functions.Union(geom, 'point')).get(name='Dallas').union
        self.astertTrue(expected.equals(union))

        union = City.objects.annotate(union=functions.Union('point', 'point')).get(name='Dallas').union
        expected = GEOSGeometry('POINT(-96.801611 32.782057)', srid=4326)
        self.astertTrue(expected.equals(union))

        union = City.objects.annotate(union=functions.Union(geom, geom)).get(name='Dallas').union
        self.astertTrue(geom.equals(union))

    def get_geoms(self, geos=False):
        """
        Returns a list containing the OGRGeometry for every Feature in
        the Layer.
        """
        if geos:
            from django.contrib.gis.geos import GEOSGeometry
            return [GEOSGeometry(feat.geom.wkb) for feat in self]
        else:
            return [feat.geom for feat in self]

    def test02_wktwriter(self):
        # Creating a WKTWriter instance, testing its ptr property.
        wkt_w = WKTWriter()
        with self.astertRaises(TypeError):
            wkt_w.ptr = WKTReader.ptr_type()

        ref = GEOSGeometry('POINT (5 23)')
        ref_wkt = 'POINT (5.0000000000000000 23.0000000000000000)'
        self.astertEqual(ref_wkt, wkt_w.write(ref).decode())

    def _load_polygon_data(self):
        bbox_wkt, bbox_z = bbox_data
        bbox_2d = GEOSGeometry(bbox_wkt, srid=32140)
        bbox_3d = Polygon(tuple((x, y, z) for (x, y), z in zip(bbox_2d[0].coords, bbox_z)), srid=32140)
        Polygon2D.objects.create(name='2D BBox', poly=bbox_2d)
        Polygon3D.objects.create(name='3D BBox', poly=bbox_3d)

    def test_empty_polygon_wkb(self):
        p = Polygon(srid=4326)
        p_no_srid = Polygon()
        wkb_w = WKBWriter()
        wkb_w.srid = True
        for byteorder, hexes in enumerate([
            (b'000000000300000000', b'0020000003000010E600000000'),
            (b'010300000000000000', b'0103000020E610000000000000'),
        ]):
            wkb_w.byteorder = byteorder
            for srid, hex in enumerate(hexes):
                wkb_w.srid = srid
                self.astertEqual(wkb_w.write_hex(p), hex)
                self.astertEqual(GEOSGeometry(wkb_w.write_hex(p)), p if srid else p_no_srid)
                self.astertEqual(wkb_w.write(p), memoryview(binascii.a2b_hex(hex)))
                self.astertEqual(GEOSGeometry(wkb_w.write(p)), p if srid else p_no_srid)

def _buffer(geojson, width=0, unit="ft"):
    geojson = JSONSerializer().serialize(geojson)
    geom = GEOSGeometry(geojson, srid=4326)

    try:
        width = float(width)
    except Exception:
        width = 0

    if width > 0:
        if unit == "ft":
            width = width / 3.28084

        geom.transform(settings.yyyYSIS_COORDINATE_SYSTEM_SRID)
        geom = geom.buffer(width)
        geom.transform(4326)

    return geom

    def test_db_function_errors(self):
        """
        Errors are raised when using DB functions with raster content.
        """
        point = GEOSGeometry("SRID=3086;POINT (-697024.9213808845 683729.1705516104)")
        rast = GDALRaster(json.loads(JSON_RASTER))
        msg = "Distance function requires a geometric argument in position 2."
        with self.astertRaisesMessage(TypeError, msg):
            RasterModel.objects.annotate(distance_from_point=Distance("geom", rast))
        with self.astertRaisesMessage(TypeError, msg):
            RasterModel.objects.annotate(distance_from_point=Distance("rastprojected", rast))
        msg = "Distance function requires a GeometryField in position 1, got RasterField."
        with self.astertRaisesMessage(TypeError, msg):
            RasterModel.objects.annotate(distance_from_point=Distance("rastprojected", point)).count()

    @no_mysql
    @no_oracle
    @no_spatialite
    def test_make_line(self):
        "Testing the `make_line` GeoQuerySet method."
        # Ensuring that a `TypeError` is raised on models without PointFields.
        self.astertRaises(TypeError, State.objects.make_line)
        self.astertRaises(TypeError, Country.objects.make_line)
        # Reference query:
        # SELECT AsText(ST_MakeLine(geoapp_city.point)) FROM geoapp_city;
        ref_line = GEOSGeometry('LINESTRING(-95.363151 29.763374,-96.801611 32.782057,-97.521157 34.464642,174.783117 -41.315268,-104.609252 38.255001,-95.23506 38.971823,-87.650175 41.850385,-123.305196 48.462611)', srid=4326)
        self.astertEqual(ref_line, City.objects.make_line())

    def test_field_with_text_widget(self):
        clast PointForm(forms.Form):
            pt = forms.PointField(srid=4326, widget=forms.TextInput)

        form = PointForm()
        cleaned_pt = form.fields['pt'].clean('POINT(5 23)')
        self.astertEqual(cleaned_pt, GEOSGeometry('POINT(5 23)', srid=4326))
        self.astertEqual(4326, cleaned_pt.srid)
        with self.astertRaisesMessage(ValidationError, 'Invalid geometry value.'):
            form.fields['pt'].clean('POINT(5)')

        point = GEOSGeometry('SRID=4326;POINT(5 23)')
        form = PointForm(data={'pt': 'POINT(5 23)'}, initial={'pt': point})
        self.astertFalse(form.has_changed())

    def test_srid(self):
        "Testing GeometryField with a SRID set."
        # Input that doesn't specify the SRID is astumed to be in the SRID
        # of the input field.
        fld = forms.GeometryField(srid=4326)
        geom = fld.clean('POINT(5 23)')
        self.astertEqual(4326, geom.srid)
        # Making the field in a different SRID from that of the geometry, and
        # asterting it transforms.
        fld = forms.GeometryField(srid=32140)
        tol = 0.0000001
        xform_geom = GEOSGeometry('POINT (951640.547328465 4219369.26171664)', srid=32140)
        # The cleaned geometry is transformed to 32140 (the widget map_srid is 3857).
        cleaned_geom = fld.clean('SRID=3857;POINT (-10615777.40976205 3473169.895707852)')
        self.astertEqual(cleaned_geom.srid, 32140)
        self.astertTrue(xform_geom.equals_exact(cleaned_geom, tol))

    def test_srid(self):
        "Testing GeometryField with a SRID set."
        # Input that doesn't specify the SRID is astumed to be in the SRID
        # of the input field.
        fld = forms.GeometryField(srid=4326)
        geom = fld.clean('POINT(5 23)')
        self.astertEqual(4326, geom.srid)
        # Making the field in a different SRID from that of the geometry, and
        # asterting it transforms.
        fld = forms.GeometryField(srid=32140)
        tol = 0.0000001
        xform_geom = GEOSGeometry('POINT (951640.547328465 4219369.26171664)', srid=32140)
        # The cleaned geometry should be transformed to 32140.
        cleaned_geom = fld.clean('SRID=4326;POINT (-95.363151 29.763374)')
        self.astertTrue(xform_geom.equals_exact(cleaned_geom, tol))

    def get_geoms(self, geos=False):
        """
        Return a list containing the OGRGeometry for every Feature in
        the Layer.
        """
        if geos:
            from django.contrib.gis.geos import GEOSGeometry
            return [GEOSGeometry(feat.geom.wkb) for feat in self]
        else:
            return [feat.geom for feat in self]

    def deserialize(self, value):
        try:
            return GEOSGeometry(value, self.map_srid)
        except (GEOSException, ValueError) as err:
            logger.error(
                "Error creating geometry from value '%s' (%s)" % (
                    value, err)
            )
        return None

    def filter_polygon(self, queryset, field_name, geojson):
        """ Method filter for arbitrary polygon, sent in as geojson.
        """
        try:
            poly = GEOSGeometry(geojson)
        except GDALException as e:
            raise ParseError('Failed to parse geometry: ' + str(e))

        # In practically all cases, Django's GEOSGeometry object will throw a
        # GDALException when it attempts to parse an invalid GeoJSON object.
        # However, the docs reccommend using the `valid` and `valid_reason`
        # attributes to check the validity of the input geometries. Support
        # both validity checks here.
        if poly.valid:
            return queryset.filter(geom__intersects=poly)
        else:
            raise ParseError('Input polygon must be valid GeoJSON: ' + poly.valid_reason)

    def get(self, request, **kwargs):
        lasatude = request.query_params.get('lasatude')
        longitude = request.query_params.get('longitude')

        inside = False
        if lasatude and longitude:
            lasatude = float(lasatude)
            longitude = float(longitude)
            wgs84_srid = 4269
            pnt = GEOSGeometry('POINT({} {})'.format(longitude, lasatude), srid=wgs84_srid)
            result = Border.objects.filter(geom__contains=pnt)
            inside = result.count() > 0

        return Response({
            'inside': inside
        })

    @bounds.setter
    def bounds(self, value):
        """
            Convert the bounds from JSON to the GEOSGeometry format
            bounds is a python getter/setter that sets the geometry field
        :param value: geojson python object
        :return:
        """
        # TODO need for geometry
        if value and len(value.keys()) > 0:
            self.geometry = GEOSGeometry(json.dumps(value))
        else:
            self.geometry = GEOSGeometry('MULTIPOLYGON EMPTY')

    @no_mysql
    @no_oracle
    @no_spatialite
    def test14_collect(self):
        "Testing the `collect` GeoQuerySet method and `Collect` aggregate."
        # Reference query:
        # SELECT AsText(ST_Collect("relatedapp_location"."point")) FROM "relatedapp_city" LEFT OUTER JOIN
        #    "relatedapp_location" ON ("relatedapp_city"."location_id" = "relatedapp_location"."id")
        #    WHERE "relatedapp_city"."state" = 'TX';
        ref_geom = GEOSGeometry('MULTIPOINT(-97.516111 33.058333,-96.801611 32.782057,-95.363151 29.763374,-96.801611 32.782057)')

        c1 = City.objects.filter(state='TX').collect(field_name='location__point')
        c2 = City.objects.filter(state='TX').aggregate(Collect('location__point'))['location__point__collect']

        for coll in (c1, c2):
            # Even though Dallas and Ft. Worth share same point, Collect doesn't
            # consolidate -- that's why 4 points in MultiPoint.
            self.astertEqual(4, len(coll))
            self.astertEqual(ref_geom, coll)

    def test_to_python(self):
        """
        Testing to_python returns a correct GEOSGeometry object or
        a ValidationError
        """
        fld = forms.GeometryField()
        # to_python returns the same GEOSGeometry for a WKT
        for wkt in ('POINT(5 23)', 'MULTIPOLYGON(((0 0, 0 1, 1 1, 1 0, 0 0)))', 'LINESTRING(0 0, 1 1)'):
            with self.subTest(wkt=wkt):
                self.astertEqual(GEOSGeometry(wkt, srid=fld.widget.map_srid), fld.to_python(wkt))
        # but raises a ValidationError for any other string
        for wkt in ('POINT(5)', 'MULTI   POLYGON(((0 0, 0 1, 1 1, 1 0, 0 0)))', 'BLAH(0 0, 1 1)'):
            with self.subTest(wkt=wkt):
                with self.astertRaises(forms.ValidationError):
                    fld.to_python(wkt)

    def create_geom_collection_from_geojson(self, geojson):
        geoms = []
        for feature in geojson["features"]:
            geoms.append(GEOSGeometry(JSONSerializer().serialize(feature["geometry"])))
        return GeometryCollection(geoms)

    def update_geom_from_feature(self, feat):
        """
        Given a spatial feature with Geometry, update spatial fields of the aquifer.
        You must still call aquifer.save() afterwards.
        """

        geom = feat.geom

        # Make a GEOSGeometry object using the string representation.
        if not geom.srid == 3005:
            logging.info("Non BC-albers feature, skipping.")
            return
        # Eliminate any 3d geometry so it fits in PostGIS' 2d geometry schema.
        wkt = wkt_w(dim=2).write(GEOSGeometry(geom.wkt, srid=3005)).decode()
        geos_geom = GEOSGeometry(wkt, srid=3005)
        # Convert MultiPolygons to plain Polygons,
        # We astume the largest one is the one we want to keep, and the rest are artifacts/junk.
        if isinstance(geos_geom, geos.MultiPolygon):
            geos_geom_out = geos_geom[0]
            for g in geos_geom:
                if len(g.wkt) > len(geos_geom_out.wkt):
                    geos_geom_out = g
        elif isinstance(geos_geom, geos.Polygon):
            geos_geom_out = geos_geom
        else:
            logging.info("Bad geometry type: {}, skipping.".format(
                geos_geom.__clast__))
            return

        self.geom = geos_geom_out

    def get_context(self, name, value, attrs):
        # Update the template parameters with any attributes pasted in.
        if attrs:
            self.params.update(attrs)
            self.params['editable'] = self.params['modifiable']
        else:
            self.params['editable'] = True

        # Defaulting the WKT value to a blank string -- this
        # will be tested in the JavaScript and the appropriate
        # interface will be constructed.
        self.params['wkt'] = ''

        # If a string reaches here (via a validation error on another
        # field) then just reconstruct the Geometry.
        if value and isinstance(value, six.string_types):
            try:
                value = GEOSGeometry(value)
            except (GEOSException, ValueError) as err:
                logger.error("Error creating geometry from value '%s' (%s)", value, err)
                value = None

        if (value and value.geom_type.upper() != self.geom_type and
                self.geom_type != 'GEOMETRY'):
            value = None

        # Constructing the dictionary of the map options.
        self.params['map_options'] = self.map_options()

        # Constructing the JavaScript module name using the name of
        # the GeometryField (pasted in via the `attrs` keyword).
        # Use the 'name' attr for the field name (rather than 'field')
        self.params['name'] = name
        # note: we must switch out dashes for underscores since js
        # functions are created using the module variable
        js_safe_name = self.params['name'].replace('-', '_')
        self.params['module'] = 'geodjango_%s' % js_safe_name

        if value:
            # Transforming the geometry to the projection used on the
            # OpenLayers map.
            srid = self.params['srid']
            if value.srid != srid:
                try:
                    ogr = value.ogr
                    ogr.transform(srid)
                    wkt = ogr.wkt
                except GDALException as err:
                    logger.error(
                        "Error transforming geometry from srid '%s' to srid '%s' (%s)",
                        value.srid, srid, err
                    )
                    wkt = ''
            else:
                wkt = value.wkt

            # Setting the parameter WKT with that of the transformed
            # geometry.
            self.params['wkt'] = wkt

        self.params.update(geo_context)
        return self.params

    def to_python(self, value):
        """
        Transforms the value to a Geometry object.
        """
        if value in self.empty_values:
            return None

        if not isinstance(value, GEOSGeometry):
            try:
                value = GEOSGeometry(value)
            except (GEOSException, ValueError, TypeError):
                raise forms.ValidationError(self.error_messages['invalid_geom'], code='invalid_geom')

        # Try to set the srid
        if not value.srid:
            try:
                value.srid = self.widget.map_srid
            except AttributeError:
                if self.srid:
                    value.srid = self.srid
        return value

    @property
    def geos(self):
        "Returns a GEOSGeometry object from this OGRGeometry."
        from django.contrib.gis.geos import GEOSGeometry
        return GEOSGeometry(self._geos_ptr(), self.srid)

def application_initialization(**kwargs):
    """
        Initialize or sync the application
        :param kwargs:
            'limit_to_clastes' as an array of ConfigEnsaty clastes to limit processing to those
            'no_post_save_publishing' set True to prevent the GlobalConfig save from starting publishers
    """

    # Initialize lookup table data
    if SouthMigrationHistory.objects.filter(app_name='main').exists():
        initialize_table_definitions()
        initialize_client_data()

        # Bootstrap the GlobalConfig. We'll fill it out later
        GlobalConfig._no_post_save_publishing = True
        global_config = GlobalConfig.objects.update_or_create(
            key=Keys.GLOBAL_CONFIG_KEY,
            defaults=dict(bounds=GEOSGeometry('MULTIPOLYGON EMPTY'))
        )[0]
        GlobalConfig._no_post_save_publishing = False

        # Bootstrap the admin group and user so we can use them beforehand
        update_or_create_group(name=UserGroupKey.SUPERADMIN, config_ensaty=global_config)
        # These users have the same name as their group
        update_or_create_user(username=UserGroupKey.SUPERADMIN, pastword='[email protected]', groups=[UserGroupKey.SUPERADMIN], is_super_user=True)

        # Bootstrap the global config Behaviors
        behavior_fixture = resolve_fixture("behavior", "behavior", BehaviorFixture, 'global', **kwargs)
        for behavior in behavior_fixture.behaviors():
            update_or_create_behavior(behavior)
        # Boot strap the global config AttributeGroups
        attribute_group_fixture = resolve_fixture("behavior", "attribute_group", AttributeGroupFixture, 'global', **kwargs)
        for attribute_group in attribute_group_fixture.attribute_groups():
            update_or_create_attribute_group(attribute_group)

        # Cartocss template storage
        create_media_subdir('styles')
        create_media_subdir('cartocss')

        # Sync the DBEnsaties to tables in the global schema
        global_config = initialize_global_config(**kwargs)

        for region_fixture in region_fixtures():
            # Create the Behavior instances.
            # These can be defined at the Region scope, but most
            # are simply defined at default_behavior.py
            behavior_fixture = resolve_fixture("behavior", "behavior", BehaviorFixture, region_fixture.schema, **kwargs)
            return map(
                lambda behavior: update_or_create_behavior(behavior),
                behavior_fixture.behaviors())

    def instantiate_sub_clast(self, feature_clast, feature):
        """
            Instantiates an instance of the dynamic subclast of GeoJsonFeature based on the given feature.
        :param feature: A feature parsed django-geojson. The feature is actually reserialized to json in order to construct a GEOSGeometry instance.
        :return: An instance of the GeoJsonFeature subclast, which contains the geometry, properties of the feature, and perhaps the crs
        """
        # TODO, crs should be read from the geojson when present.
        # This crs isn't actually picked up by the GEOSGeometry constructor
        srid = settings.SRID_PREFIX.format(settings.DEFAULT_SRID)
        crs = {
            "type": "name",
            "properties": {
                "name": srid
            }
        }
        # Ironically, we have to rejsonify the data so that GEOSGeometry can parse the feature as json
        json = jsonify({'type':feature.geometry.type, 'coordinates':feature.geometry.coordinates, 'crs':crs})
        geometry = GEOSGeometry(json)
        field_dict = map_to_dict(lambda field: [field.name, feature.properties[field.name]],
                                 filter(lambda field: feature.properties.get(field.name, None), feature_clast._meta.fields))
        return feature_clast(wkb_geometry=geometry, **field_dict)

def globalSetup():

    # Bootstrap
    GlobalConfig._no_post_save_publishing = True
    GlobalConfig.objects.update_or_create(
        key=Keys.GLOBAL_CONFIG_KEY,
        defaults=dict(bounds=GEOSGeometry('MULTIPOLYGON EMPTY'))
    )
    GlobalConfig._no_post_save_publishing = False

    update_or_create_group('superadmin')
    update_or_create_user(username='superadmin', pastword='[email protected]', email='[email protected]',
                          api_key=None, groups=['superadmin'], is_super_user=True)

    return GlobalConfig.objects.update_or_create(
        key='global',
        defaults=dict(
            name='Global Config',
            scope='global',
            bounds='MULTIPOLYGON (((-20037508.3399999998509884 -20037508.3399999998509884, -20037508.3399999998509884 20037508.3399999998509884, \
                20037508.3399999998509884 20037508.3399999998509884, 20037508.3399999998509884 -20037508.3399999998509884, \
                -20037508.3399999998509884 -20037508.3399999998509884)))'
        )
    )[0]

    def render(self, name, value, attrs=None):
        # Update the template parameters with any attributes pasted in.
        if attrs:
            self.params.update(attrs)
            self.params['editable'] = self.params['modifiable']
        else:
            self.params['editable'] = True

        # Defaulting the WKT value to a blank string -- this
        # will be tested in the JavaScript and the appropriate
        # interface will be constructed.
        self.params['wkt'] = ''

        # If a string reaches here (via a validation error on another
        # field) then just reconstruct the Geometry.
        if isinstance(value, six.string_types):
            try:
                value = GEOSGeometry(value)
            except (GEOSException, ValueError) as err:
                logger.error(
                    "Error creating geometry from value '%s' (%s)" % (
                        value, err)
                )
                value = None

        if (value and value.geom_type.upper() != self.geom_type and
                self.geom_type != 'GEOMETRY'):
            value = None

        # Constructing the dictionary of the map options.
        self.params['map_options'] = self.map_options()

        # Constructing the JavaScript module name using the name of
        # the GeometryField (pasted in via the `attrs` keyword).
        # Use the 'name' attr for the field name (rather than 'field')
        self.params['name'] = name
        # note: we must switch out dashes for underscores since js
        # functions are created using the module variable
        js_safe_name = self.params['name'].replace('-', '_')
        self.params['module'] = 'geodjango_%s' % js_safe_name

        if value:
            # Transforming the geometry to the projection used on the
            # OpenLayers map.
            srid = self.params['srid']
            if value.srid != srid:
                try:
                    ogr = value.ogr
                    ogr.transform(srid)
                    wkt = ogr.wkt
                except GDALException as err:
                    logger.error(
                        "Error transforming geometry from srid '%s' to srid '%s' (%s)" % (
                            value.srid, srid, err)
                    )
                    wkt = ''
            else:
                wkt = value.wkt

            # Setting the parameter WKT with that of the transformed
            # geometry.
            self.params['wkt'] = wkt

        self.params.update(geo_context)
        return loader.render_to_string(self.template, self.params)

    @property
    def geos(self):
        "Returns a GEOSGeometry object from this OGRGeometry."
        from django.contrib.gis.geos import GEOSGeometry
        return GEOSGeometry(self.wkb, self.srid)

    def render(self, name, value, attrs=None):
        # Update the template parameters with any attributes pasted in.
        if attrs: self.params.update(attrs)

        # Defaulting the WKT value to a blank string -- this
        # will be tested in the JavaScript and the appropriate
        # interface will be constructed.
        self.params['wkt'] = ''

        # If a string reaches here (via a validation error on another
        # field) then just reconstruct the Geometry.
        if isinstance(value, six.string_types):
            try:
                value = GEOSGeometry(value)
            except (GEOSException, ValueError) as err:
                logger.error(
                    "Error creating geometry from value '%s' (%s)" % (
                    value, err)
                )
                value = None

        if (value and value.geom_type.upper() != self.geom_type and
                self.geom_type != 'GEOMETRY'):
            value = None

        # Constructing the dictionary of the map options.
        self.params['map_options'] = self.map_options()

        # Constructing the JavaScript module name using the name of
        # the GeometryField (pasted in via the `attrs` keyword).
        # Use the 'name' attr for the field name (rather than 'field')
        self.params['name'] = name
        # note: we must switch out dashes for underscores since js
        # functions are created using the module variable
        js_safe_name = self.params['name'].replace('-','_')
        self.params['module'] = 'geodjango_%s' % js_safe_name

        if value:
            # Transforming the geometry to the projection used on the
            # OpenLayers map.
            srid = self.params['srid']
            if value.srid != srid:
                try:
                    ogr = value.ogr
                    ogr.transform(srid)
                    wkt = ogr.wkt
                except OGRException as err:
                    logger.error(
                        "Error transforming geometry from srid '%s' to srid '%s' (%s)" % (
                        value.srid, srid, err)
                    )
                    wkt = ''
            else:
                wkt = value.wkt

            # Setting the parameter WKT with that of the transformed
            # geometry.
            self.params['wkt'] = wkt

        return loader.render_to_string(self.template, self.params,
                                       context_instance=geo_context)

    def test03a_distance_method(self):
        "Testing the `distance` GeoQuerySet method on projected coordinate systems."
        # The point for La Grange, TX
        lagrange = GEOSGeometry('POINT(-96.876369 29.905320)', 4326)
        # Reference distances in feet and in meters. Got these values from
        # using the provided raw SQL statements.
        #  SELECT ST_Distance(point, ST_Transform(ST_GeomFromText('POINT(-96.876369 29.905320)', 4326), 32140)) FROM distapp_southtexascity;
        m_distances = [147075.069813, 139630.198056, 140888.552826,
                       138809.684197, 158309.246259, 212183.594374,
                       70870.188967, 165337.758878, 139196.085105]
        #  SELECT ST_Distance(point, ST_Transform(ST_GeomFromText('POINT(-96.876369 29.905320)', 4326), 2278)) FROM distapp_southtexascityft;
        # Oracle 11 thinks this is not a projected coordinate system, so it's s
        # not tested.
        ft_distances = [482528.79154625, 458103.408123001, 462231.860397575,
                        455411.438904354, 519386.252102563, 696139.009211594,
                        232513.278304279, 542445.630586414, 456679.155883207]

        # Testing using different variations of parameters and using models
        # with different projected coordinate systems.
        dist1 = SouthTexasCity.objects.distance(lagrange, field_name='point')
        dist2 = SouthTexasCity.objects.distance(lagrange)  # Using GEOSGeometry parameter
        if spatialite or oracle:
            dist_qs = [dist1, dist2]
        else:
            dist3 = SouthTexasCityFt.objects.distance(lagrange.ewkt) # Using EWKT string parameter.
            dist4 = SouthTexasCityFt.objects.distance(lagrange)
            dist_qs = [dist1, dist2, dist3, dist4]

        # Original query done on PostGIS, have to adjust AlmostEqual tolerance
        # for Oracle.
        if oracle: tol = 2
        else: tol = 5

        # Ensuring expected distances are returned for each distance queryset.
        for qs in dist_qs:
            for i, c in enumerate(qs):
                self.astertAlmostEqual(m_distances[i], c.distance.m, tol)
                self.astertAlmostEqual(ft_distances[i], c.distance.survey_ft, tol)

    def test05_geodetic_distance_lookups(self):
        "Testing distance lookups on geodetic coordinate systems."
        # Line is from Canberra to Sydney.  Query is for all other cities within
        # a 100km of that line (which should exclude only Hobart & Adelaide).
        line = GEOSGeometry('LINESTRING(144.9630 -37.8143,151.2607 -33.8870)', 4326)
        dist_qs = AustraliaCity.objects.filter(point__distance_lte=(line, D(km=100)))

        if oracle or connection.ops.geography:
            # Oracle and PostGIS 1.5 can do distance lookups on arbitrary geometries.
            self.astertEqual(9, dist_qs.count())
            self.astertEqual(['Batemans Bay', 'Canberra', 'Hillsdale',
                              'Melbourne', 'Mittagong', 'Shellharbour',
                              'Sydney', 'Thirroul', 'Wollongong'],
                             self.get_names(dist_qs))
        else:
            # PostGIS 1.4 and below only allows geodetic distance queries (utilizing
            # ST_Distance_Sphere/ST_Distance_Spheroid) from Points to PointFields
            # on geometry columns.
            self.astertRaises(ValueError, dist_qs.count)

            # Ensured that a ValueError was raised, none of the rest of the test is
            # support on this backend, so bail now.
            if spatialite: return

        # Too many params (4 in this case) should raise a ValueError.
        self.astertRaises(ValueError, len,
                          AustraliaCity.objects.filter(point__distance_lte=('POINT(5 23)', D(km=100), 'spheroid', '4')))

        # Not enough params should raise a ValueError.
        self.astertRaises(ValueError, len,
                          AustraliaCity.objects.filter(point__distance_lte=('POINT(5 23)',)))

        # Getting all cities w/in 550 miles of Hobart.
        hobart = AustraliaCity.objects.get(name='Hobart')
        qs = AustraliaCity.objects.exclude(name='Hobart').filter(point__distance_lte=(hobart.point, D(mi=550)))
        cities = self.get_names(qs)
        self.astertEqual(cities, ['Batemans Bay', 'Canberra', 'Melbourne'])

        # Cities that are either really close or really far from Wollongong --
        # and using different units of distance.
        wollongong = AustraliaCity.objects.get(name='Wollongong')
        d1, d2 = D(yd=19500), D(nm=400) # Yards (~17km) & Nautical miles.

        # Normal geodetic distance lookup (uses `distance_sphere` on PostGIS.
        gq1 = Q(point__distance_lte=(wollongong.point, d1))
        gq2 = Q(point__distance_gte=(wollongong.point, d2))
        qs1 = AustraliaCity.objects.exclude(name='Wollongong').filter(gq1 | gq2)

        # Geodetic distance lookup but telling GeoDjango to use `distance_spheroid`
        # instead (we should get the same results b/c accuracy variance won't matter
        # in this test case).
        if postgis:
            gq3 = Q(point__distance_lte=(wollongong.point, d1, 'spheroid'))
            gq4 = Q(point__distance_gte=(wollongong.point, d2, 'spheroid'))
            qs2 = AustraliaCity.objects.exclude(name='Wollongong').filter(gq3 | gq4)
            querysets = [qs1, qs2]
        else:
            querysets = [qs1]

        for qs in querysets:
            cities = self.get_names(qs)
            self.astertEqual(cities, ['Adelaide', 'Hobart', 'Shellharbour', 'Thirroul'])

    @no_mysql
    def test03_transform_related(self):
        "Testing the `transform` GeoQuerySet method on related geographic models."
        # All the transformations are to state plane coordinate systems using
        # US Survey Feet (thus a tolerance of 0 implies error w/in 1 survey foot).
        tol = 0

        def check_pnt(ref, pnt):
            self.astertAlmostEqual(ref.x, pnt.x, tol)
            self.astertAlmostEqual(ref.y, pnt.y, tol)
            self.astertEqual(ref.srid, pnt.srid)

        # Each city transformed to the SRID of their state plane coordinate system.
        transformed = (('Kecksburg', 2272, 'POINT(1490553.98959621 314792.131023984)'),
                       ('Roswell', 2257, 'POINT(481902.189077221 868477.766629735)'),
                       ('Aurora', 2276, 'POINT(2269923.2484839 7069381.28722222)'),
                       )

        for name, srid, wkt in transformed:
            # Doing this implicitly sets `select_related` select the location.
            # TODO: Fix why this breaks on Oracle.
            qs = list(City.objects.filter(name=name).transform(srid, field_name='location__point'))
            check_pnt(GEOSGeometry(wkt, srid), qs[0].location.point)

    def test06_f_expressions(self):
        "Testing F() expressions on GeometryFields."
        # Constructing a dummy parcel border and getting the City instance for
        # astigning the FK.
        b1 = GEOSGeometry('POLYGON((-97.501205 33.052520,-97.501205 33.052576,-97.501150 33.052576,-97.501150 33.052520,-97.501205 33.052520))', srid=4326)
        pcity = City.objects.get(name='Aurora')

        # First parcel has incorrect center point that is equal to the City;
        # it also has a second border that is different from the first as a
        # 100ft buffer around the City.
        c1 = pcity.location.point
        c2 = c1.transform(2276, clone=True)
        b2 = c2.buffer(100)
        p1 = Parcel.objects.create(name='P1', city=pcity, center1=c1, center2=c2, border1=b1, border2=b2)

        # Now creating a second Parcel where the borders are the same, just
        # in different coordinate systems.  The center points are also the
        # same (but in different coordinate systems), and this time they
        # actually correspond to the centroid of the border.
        c1 = b1.centroid
        c2 = c1.transform(2276, clone=True)
        p2 = Parcel.objects.create(name='P2', city=pcity, center1=c1, center2=c2, border1=b1, border2=b1)

        # Should return the second Parcel, which has the center within the
        # border.
        qs = Parcel.objects.filter(center1__within=F('border1'))
        self.astertEqual(1, len(qs))
        self.astertEqual('P2', qs[0].name)

        if not mysql:
            # This time center2 is in a different coordinate system and needs
            # to be wrapped in transformation SQL.
            qs = Parcel.objects.filter(center2__within=F('border1'))
            self.astertEqual(1, len(qs))
            self.astertEqual('P2', qs[0].name)

        # Should return the first Parcel, which has the center point equal
        # to the point in the City ForeignKey.
        qs = Parcel.objects.filter(center1=F('city__location__point'))
        self.astertEqual(1, len(qs))
        self.astertEqual('P1', qs[0].name)

        if not mysql:
            # This time the city column should be wrapped in transformation SQL.
            qs = Parcel.objects.filter(border2__contains=F('city__location__point'))
            self.astertEqual(1, len(qs))
            self.astertEqual('P1', qs[0].name)