Good morning,<br><br>Here is a code snippet I'm using in gdal-python to transform a point from GDM2000 (geodetic) to BRSO (projected):<br><br><pre><font size="2"><i><span style="font-family: verdana,sans-serif;">"""</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> transform a pair of coordinates from the input CRS</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> to the output GEOGCS and output PROJCS</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> all codes are EPSG integers only</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> (used internally)</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> """</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> in_CRS = osr.SpatialReference()</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> out_GEOGCS = osr.SpatialReference()</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> out_PROJCS = osr.SpatialReference()</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> in_CRS.ImportFromEPSG(inCRS)</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> out_GEOGCS.ImportFromEPSG(outGEOGCS)</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> out_PROJCS.ImportFromEPSG(outPROJCS)</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> transform = osr.CoordinateTransformation(in_CRS, out_GEOGCS)</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> coords = transform.TransformPoint(self.lon_value,self.lat_value)</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> self.latitude = coords[1]</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> self.longitude = coords[0] </span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> proj_transform = osr.CoordinateTransformation(out_GEOGCS, out_PROJCS)</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> proj_coords = proj_transform.TransformPoint(float(self.longitude), float(self.latitude))</span><br style="font-family: verdana,sans-serif;"><span style="font-family: verdana,sans-serif;"> self.easting = proj_coords[0]</span><br style="font-family: verdana,sans-serif;">
<span style="font-family: verdana,sans-serif;"> self.northing = proj_coords[1]</span></i></font><br></pre>I'm using EPSG codes 4742 (GDM2000) and 3376 (BRSO) with this test point:<br>ϕ = 6° 52′ 42.45022″ N<br>λ = 116° 50′ 47.58835″ E<br>
<br>This should give E/N of:<br>793704.631 mE 762081.047 mN<br><br>but I am getting 204217.0mE 318529.9mN. I have not made any changes to the GDAL files gcs.csv or pcs.csv<br><br>Can anyone see what mistake I am making please?<br>
<br>Thanks, Roger<br><br clear="all"><br>