Adding metadata to a geographic gridded data-set and projecting it

Question

I asked a question before that nobody answered...but I kind of figured it out myself. But I still remain with some things that I just don't seem to get right. Can anybody please find the time to help me? It really seems like something trivial... Thanks!

-

I want to transform a geographicly gridded latlon data-set into a polar stereographic projection with gdalwarp. For this I need to add the correct metadata to my gridded data-set, and I can't seem to get this right.

It is regularly spaced with 0.04 degrees from 0:360 longitude and -86:-60 Latitude (Antarctica). Only thing I have is the grid, so it comes without a header. I however have added the current header to the dataset:

ncols        9001
nrows        651
xllcorner    0
yllcorner    -86.000000000000
dx           0.04
dy           0.04

The first problem I get is that when I gdalinfo it I get:

Corner Coordinates:
Upper Left  (   0.0000000, -59.9600000) 
Lower Left  (   0.0000000, -86.0000000) 
Upper Right (     360.040,     -59.960) 
Lower Right (     360.040,     -86.000) 
Center      ( 180.0200000, -72.9800000) 

Which seems wrong (as it is 0.04 off). Possibly I have to use a gridsize of 9000/650 instead of 9001/651, but also this is confusing to me as I then omit a line. Do I set the header correctly? Other than defining the grid and the coordinates like this, what do I need to do with the projection data? Need I add -s_srs WGS84 (this is what I do) or something like -s_srs EPSG:3031 (South Polar Stereographic)? This gives me the following:

Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (0.000000000000000,-59.960000000000001)
Pixel Size = (0.040000000000000,-0.040000000000000)

I want to finally project it onto the following system:

PROJCS["unnamed",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",-71],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
 Origin = (-2800500.000000000000000,2800500.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
 INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-2800500.000, 2800500.000) ( 45d 0' 0.00"W, 54d40' 3.36"S)
Lower Left  (-2800500.000,-2800500.000) (135d 0' 0.00"W, 54d40' 3.36"S)
Upper Right ( 2800500.000, 2800500.000) ( 45d 0' 0.00"E, 54d40' 3.36"S)
Lower Right ( 2800500.000,-2800500.000) (135d 0' 0.00"E, 54d40' 3.36"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E, 90d 0' 0.00"S)

Now there are 2 things occuring with my results when I warp it like this:

gdalwarp -s_srs WGS84 -r cubic data1.txt data2.tif

(i use either an existing data-set (I want to remap so...) or set the -te and -tr parameters, and get the same results....)

1. I get a warped file which fails to do the entire file (there is probably something going wrong at the 0/360 longitude transition. As expected it omits the value that falls at -360.04
2. I get only half of the file!! Only the right side of the projection is done...

So my questions are:

1. How do i correctly set the correct projection parameters and how do I deal with polar stereographic data? Why can't I properly set the latitude and longitude range?
2. What in earth is going wrong when it only projects the right side?

Thanks for the help!

Answer 1

I can't answer all your questions, but I have noticed a couple of things that might help. First, your apparent off-by-0.04 is because by default, GDAL assumes the centre of a pixel represents the coordinate, not the corner. So by shifting your xllcorner and yllcorner by -0.02 will end up with the right coverage.

Secondly, most CRSs that I'm aware of have a range between -180 and +180 degrees longitude, whereas your bounds are 0 - 360. Try setting your xllcorner value to -180.02 (assuming you're using Greenwich as your prime meridian).