[GIS] GDAL to project from pixels to geo coordinates

coordinate systemgdal

I am using GDAL library in C++. On my dataset, I have obtained geo transform with GetGeoTransform to GT array (GCP is NULL).

Than I use this values for this equation (Source):

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

I also add 0.5 offset to Xpixel and Yline, to calculate coordinate in the "center" of pixel.

If I put [0,0], I got correct result in lat/lon. However, if I change Yline to something other than 0, my longitude is messed up – instead of 44.0, I got -349.2… (latitude is still correct). What am I doing wrong ?

My data informations –
origin, pixel size and top left are taken from GT array used in equation above; projection is GetProjectionRef

Projection is `GEOGCS["Coordinate System imported from GRIB file",DATUM["unknown",SPHEROID["Sphere",6371200,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
Origin = (0,0) - GT(2), GT(4)
Pixel Size = (0.5,-33.268) - GT(1), GT(5)
Top Left = (5.75,66.634) - GT(0), GT(3)

Best Answer

You should correct the size of your pixel in Y. -33 degree is definitely too much.

Related Solutions

[GIS] Converting projection

The shp2pgsql tool is used to load data into PostGIS. It is not meant for transforming projections. The -s flag of the command just tells shp2pgsql the SRID of your source data. It doesn't transform the projection at all. You can reproject it before you load it to PostGIS or after.

To transform the projection before you load it to PostGIS, you can use GIS tools like QGIS or GDAL.

To change the projection after loading the data into PostGIS, you can use the ST_Transform command.

Btw, you have to give the shp2pgsql tool the shape file's actual projection. According to Spatialreference.org, your shape file's SRID is 2154. You might want to import it properly this time, using -s 2154 instead of -s 4326 and then perform the ST_Transform on PostGIS.

Of the two options, I think it would be simpler to convert the file before loading it to PostGIS. That way, you won't have to mess around with the spatial database much.

On Identifying your shape file's projection

Stumbled upon prj2EPSG from after reading this thread. Pasted your definition on it and got the following results:

2154 - RGF93 / Lambert-93
3944 - RGF93 / CC44
3949 - RGF93 / CC49

You might want to investigate those three.

[GIS] Problem with Converting XY Coordinates to LatLon

You don't say, but I'm assuming the PROJCS string in your question is the string being read from your dataset. There are several errors in the string. Here's what I found on SpatialReference.org:

PROJCS["WGS 84 / UTM zone 41N",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",63],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","32641"],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]

Your central meridian is wrong, as is the false easting.

Related Question