aeqd projection and lat_0 and lon_0 meaning

[regodon]

Hello,

I don't know if it's a problem related to gdal or proj.4.

It seems that lat_0 and lon_0 are being interpreted as top-left corner pixel instead of center pixel when using proj=aeqd.

Example:

<VRTDataset rasterXSize="480" rasterYSize="480">
  <SRS>+proj=aeqd +lat_0=39.18 +lon_0=-0.25 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m</SRS>
    <GeoTransform>0, 1000, 0, 0, 0, -1000</GeoTransform>
    <VRTRasterBand dataType="Byte" band="1" subClass="VRTRawRasterBand">
        <SourceFilename relativeToVRT="1">data.raw</SourceFilename>
    </VRTRasterBand>
</VRTDataset>

This is a raster centered on coordinates (-0.25, 39.18) with a pixel width of 1000 meters.

gdalinfo output for this vrt is:

Driver: VRT/Virtual Raster
Files: data.vrt
       data.raw
Size is 480, 480
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6371000,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Azimuthal_Equidistant"],
    PARAMETER["latitude_of_center",39.18],
    PARAMETER["longitude_of_center",-0.25],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) (  0d15'0.00"W, 39d10'48.00"N)
Lower Left  (       0.000, -480000.000) (  0d15'0.00"W, 34d51'47.72"N)
Upper Right (  480000.000,       0.000) (  5d18'42.59"E, 39d 2'51.56"N)
Lower Right (  480000.000, -480000.000) (  5d 0'2.46"E, 34d44'31.82"N)
Center      (  240000.000, -240000.000) (  2d27'6.99"E, 36d59'23.90"N)
Band 1 Block=480x1 Type=Byte, ColorInterp=Undefined

As you can see, "Center" appears as (2d27'6.99"E, 36d59'23.90"N) instead of (-0.25, 39.18).

Is it a geotransform problem? A Proj.4 problem?

Thanks.

[regodon]

Hello again,

With this GeoTransform:

<GeoTransform>-240000, 1000, 0, 240000, 0, -1000</GeoTransform>

i get this gdalinfo output:

Driver: GTiff/GeoTIFF
Files: radar.tiff
Size is 480, 480
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6371000,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Azimuthal_Equidistant"],
    PARAMETER["latitude_of_center",39.18],
    PARAMETER["longitude_of_center",-0.25],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-240000.000000000000000,240000.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -240000.000,  240000.000) (  3d 7'22.77"W, 41d18'13.17"N)
Lower Left  ( -240000.000, -240000.000) (  2d57'6.99"W, 36d59'23.90"N)
Upper Right (  240000.000,  240000.000) (  2d37'22.77"E, 41d18'13.17"N)
Lower Right (  240000.000, -240000.000) (  2d27'6.99"E, 36d59'23.90"N)
Center      (   0.0000000,   0.0000000) (  0d15'0.00"W, 39d10'48.00"N)
Band 1 Block=480x17 Type=Byte, ColorInterp=Gray

Layer is displayed right now. I don't know if there is a bug or just a lack of documentation.

Please close this ticket if you find that nothing needs to be done.

Thanks.

[Even Rouault]

This is not a bug (if you're not sure if a behaviour is a bug or not, use the gdal-dev mailing list first). In your first case, your geotransform <GeoTransform>0, 1000, 0, 0, 0, -1000</GeoTransform> with a raster of size 480x480 leads to a center at coordinates (0+240*1000,0-240*1000), which is not (0,0) so it is logical you don't get the (lon_0, lat_0) values. In the second case, the center is at coordinates (-240000+240*1000,2400000-240*1000)=(0,0)