Search code examples
pythongdal

Python (GDAL): projection conversion from WGS84 to NZTM2000 is not correct


I am newbie in Python and GDAL. I was converting the following file into NZTM2000 coordinate converter. I imported the file as csv and then used Lat/Long to convert it into an XY output. However, the result is not correct to me. I don't understand the output unit/scale. I was expecting something of 6 digits.

Script:

"reads a CSV file and converts its coordinates"
import os
import csv
from osgeo import gdal
from osgeo import osr
from osgeo import ogr

# Coordinate Reference System (CRS)
SourceEPSG = 4326  
TargetEPSG = 2193

source = osr.SpatialReference()
source.ImportFromEPSG(SourceEPSG)

target = osr.SpatialReference()
target.ImportFromEPSG(TargetEPSG)


# Input file details
fullpath = os.path.abspath("\ew0001\NavigationsCMP")

def CRSTransform(Lat, Long):
    transform = osr.CoordinateTransformation(source, target)
    point = ogr.CreateGeometryFromWkt("POINT ("+Lat+" "+Long+")")
    point.Transform(transform)
    print point.GetX(), "   ", point.GetY()


print "Reading CSV"
inCSV = csv.DictReader(open(fullpath+".csv"))
for row in inCSV:
    lat = row['Lat']
    long = row['Long']
    CRSTransform(lat, long)

Input:

Lat Long    CMP Year    Month   Day Hours   Mins    Sec XXX Line    Vintage
-44.419134  172.243651  264 2000    1   23  6   11  10  180 10  EW0001
-44.419176  172.243706  265 2000    1   23  6   11  12  681 10  EW0001
-44.419214  172.243759  266 2000    1   23  6   11  15  181 10  EW0001
-44.419259  172.24382   267 2000    1   23  6   11  17  711 10  EW0001

Output:

-44.419134     172.243651
-44.419176     172.243706
-44.419214     172.243759
-44.419259     172.24382

Solution

  • Original answer for GDAL <3

    You need to swap the axis order, so that X = Longitude, Y = Latitude.

    def CRSTransform(Lat, Long):
        transform = osr.CoordinateTransformation(source, target)
        point = ogr.Geometry(ogr.wkbPoint)
        point.SetPoint_2D(0, float(Long), float(Lat))
        point.Transform(transform)
        print point.GetX(), "   ", point.GetY()
    
    CRSTransform(-44.419134, 172.243651)  # 1539788.868     5081294.99354
    

    Updated answer for GDAL 3+

    Axis ordering has changed from GDAL 3.0 (e.g. see GDAL's OSR API tutorial). The default order is "authority compliant", .e.g. Lat/Long. To get "traditional axis orders" of X/Y use:

    source.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    target.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    

    And if you wanted to keep Lat/Long for source, it would be the default or specified:

    source.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
    

    Be aware that EPSG:2193 has default authority Northing/Easting (or Y/X).