Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

In QGIS do I have a EPSG:4326 shapefile that shows the coordinates "-7148300 -3688500" and I need them in the format -64º -31º for further processing in a geodjango app.

However using ogrinfo I get this information from the shapefile:

Feature Count: 22
Extent: (-64.279736, -31.509785) - (-64.092643, -31.332019)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
id: Integer (10.0)

How can I convert the values? I found useful both solutions:

1) Getting the coordinates of a given feature id, within Geodjango shell, so I have a reference and work from there using a custom conversion function.

2) Know the mathematical equivalence between the two pairs of coordinates.

A snippet or idea in geodjango would be nice, but any hint will be appreciated!

Thank you in advance.

share|improve this question
    
I'm voting to delete my own answer, because I don't know enough about QGIS to figure out what may be the answer. The shapefile description looks fine. I think the QGIS "project" might be set to a different coordinate reference system and/or the units are set to meters. –  mkennedy Mar 22 '13 at 23:51
    
Thank you. Indeed, the units are set to meters. This is the CRS definition (EPSG:3857): +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs –  Martin0x777 Mar 22 '13 at 23:56
add comment

4 Answers 4

Your QGIS project CRS is set to Google Mercator, which is needed if you want Google or Openstreetmap background by using Openlayers plugin.

The layer CRS can (and should be) different from that, in your case WGS 84 in degrees.

If you want the coordinate display in degrees, but have a Google/OSM background: make a screen copy by File -> Save Picture as, add that picture to the canvas, delete the openlayers layer, Rightclick on the shapefile layer -> Set Layer CRS for project

The picture may look distorted the more away you are from the aequator, but thats what lat/lon coordinates are.

share|improve this answer
    
Thank you Andre. I think you refer to the "Set Project CRS from layer", which is a nice solution for me. However I still need to develop some lookup functionality via a python routine or shell script, so I'll keep digging. Thank you again! –  Martin0x777 Mar 25 '13 at 22:59
add comment

This is the updated solution I have found so far:

First populating the PostGIS database

import os
from django.contrib.gis.utils import LayerMapping
from models import UJ

mj = '/home/martin/data_dir/mapa/tramos.shp'
# Auto-generated `LayerMapping` dictionary for uj model
uj_mapping = {
    'id' : 'id',
    'geom' : 'MULTILINESTRING',
}

def run(verbose=True):
    lm = LayerMapping(UJ, mj, uj_mapping,
                      transform=False, encoding='utf-8')

    lm.save(strict=True, verbose=verbose)

Then Using cs2cs from the comment of @isolier I get the converted coordinates:

cs2cs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs Then I enter "-7148300 -3688500 " without quotes and get: 64d12'51.377"W 31d25'41.939"S 0.000

Then Using PostGIS

from mapa.models import UJ
from django.contrib.gis.geos import *
from django.contrib.gis.measure import D
pbt = fromstr('POINT(-64.202117 -31.476626)', srid=4326)
qs = UJ.objects.filter(geom__distance_lte = (pbt, D(km=3)))
len(qs)
for uj in qs.distance(pbt): print(uj.id, uj.distance)

I still look for a way to automate this, but I'm closer.

Thank you everyone.

share|improve this answer
add comment

Check out the PROJ.4 command line tool cs2cs. This is certainly not the Geodjango shell solution that you would prefer, but may be best suited for what you are after. Good luck!

share|improve this answer
1  
nice one thanks :) This is the command that I issue to manually enter EPSG:900113 coordinates and get EPSG:4326 coordinates: cs2cs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs Then I enter "-7148300 -3688500 " without quotes and get: 64d12'51.377"W 31d25'41.939"S 0.000 –  Martin0x777 Mar 25 '13 at 22:39
add comment
up vote 0 down vote accepted

Another approach: The Geodjango snippet :)

>>> from django.contrib.gis.gdal import SpatialReference, CoordTransform
>>> ct = CoordTransform(SpatialReference(3857), SpatialReference(4326))
>>> pbt = fromstr('POINT(-7146947 -3694804)', srid=3857)
>>> pbt.transform(ct)
>>> print 'x: %s; y: %s; srid: %s' % (pbt.x, pbt.y, pbt.srid)
x: -64.2021172489; y: -31.4766257292; srid: 4326
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.