Thursday, 14 January 2016

GeoDjango: transform() not working?


I have a point: pnt = GEOSGeometry('POINT({0} {1})'.format(lon, lat)), and tried pnt.transform(900913), and it gives error:


GEOSException: Calling transform() with no SRID set is not supported

Then I tried transform() with my points in the database:


qs = T.objects.all().transform(900913)

and it works fine:


>>print qs[0].geom
>>POINT (-10887469.2136471699923277 3564159.8204717375338078)


Here's my model T:


class T(models.Model):
name = models.CharField(max_length=200)
geom = models.PointField(srid=4326, blank=True, null=True, default=None)
objects = models.GeoManager()

Any idea?



Answer



The reason why you get this error message is that the GEOSGeometry does not have a default SRID specified. When you instantiate an instance of a GEOSGeometry you have to explicitly specify the SRID of the geometry you are creating. To transform a geometry GEOS needs an "original" SRID and a "new" SRID to which it will transform to. Here is where you get the error, because there is no default, it does not know what the "original" SRID is.



This behavior is different for a GeoDjango Model like the one you defined. For a GeoDjango model there is a default srid, namely 4326. This is described in the GeoDjango Model reference. So for a specific instance of T or a queryset of T objects, GeoDjango assumes that the input geometry was in SRID 4326.


Here is an example how you can define a GEOSGeometry with an SRID specified. With the SRID also the transform function works fine.


(lon, lat) = ('1', '1')
pnt = GEOSGeometry('POINT({0} {1})'.format(lon, lat), srid=4326)
pnt.transform(900913)
pnt.wkt
>> u'POINT (111319.4907932723144768 111325.1428663848637370)'

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...