If you have data from the poles, avoid EPSG:3857, because that is undefined at the poles. Reprojection might fail, and the rest of the data might get lost.
Try EPSG:4326 instead. To get the full picture, include the target extent (for the Arctic region):
gdalwarp -t_srs EPSG:4326 -te -180 -90 180 90 northpsg.20141027 output.tif
and you will get your ice (under Natural Earth vector data):
EDIT
For the southern hemisphere it is a bit more work. Gdalinfo does not like the stored projection information, read as
+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-260 +k=-90 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs
which makes indeed little sense to me. So i looked into the information wgrib reports:
polar stereo: Lat1 -36.866000 Long1 -220.194000 Orient -260.000000
south pole (345 x 355) Dx 25400 Dy 25400 scan 64 mode 0
Lat1 and Long1 are the WGS84 coordinates of one corner cell (mid cell!), resolution is the same as for the arctic, but size is different.
So we have to override the projection information with gdal_translate before reprojecting.
The Orient = -260 equals a center meridian of 100°E. To get the corner coordinate, I converted the given latlong to polar:
gdaltransform -s_srs EPSG:4326 -t_srs "+proj=stere +lat_0=-90 +lat_ts=-60 +lon_0=100 +k=1 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs" <ll1south.txt >ll1-southpol.txt
with 139.806 -36.8660 as input and 3805876.67386847 4566982.49310513 as output. These have to be increased by half pixel resolution 12700.
The final reprojection is:
gdal_translate -a_srs "+proj=stere +lat_0=-90 +lat_ts=-60 +lon_0=100 +k=1 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs" -a_ullr 3818607 -4437313 -4944393 4579687 -of VRT southpsg.20141027 temp.vrt
gdalwarp -overwrite -t_srs EPSG:4326 -te -180 -90 180 90 temp.vrt outputsouth.tif
As an intermediate step, I added the .vrt to a QGIS project with the polar stereo projection. For a first guess I took +-nx * Dx/2 and +-ny * Dy/2 as extent and played with those values until it was fitting well.
You see that the reference point is in the upper right corner. This picture helps to calculate the other corner coordinates exactly by subtracting nx * Dx and ny * Dy.
This is the final picture in EPSG:4326:
You find the reference point as a little yellow dot to the south of Australia.
The projection methods commonly used in Australia are Albers equal area, Lambert conformal conic and transverse mercator.
Since the first two need one or two latitudes as parameter, I assume that transverse mercator will fit. UTM uses the same method.
You may take EPSG:3113 GDA94/BCSG02 as an example:
+proj=tmerc +lat_0=-28 +lon_0=153 +k=0.9999900000000001 +x_0=50000 +y_0=100000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
and exchange the parameters with those you have got. Make sure to use decimal degrees for the central meridian. lat_0 seems to be missing, so you might try the equator with +lat_0=0.
Basically, the units could be something else than meters, but this seems to be unusual in Australia.
From the extended information you give, the projection string should be:
+proj=tmerc +lat_0=0 +lon_0=120.95 +k=0.9999 +x_0=446904.02 +y_0=2879827.84 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Your Origin point leads me to a point east of the town named Newman, West Australia.
and the full command line:
gdalwarp -s_srs "+proj=tmerc +lat_0=0 +lon_0=120.95 +k=0.9999 +x_0=4 46904.02 +y_0=2879827.84 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -t_srs EPSG:28350 im.tif out.tif
Best Answer
According to spatilreference.org, your proj4 string is not correct. -a_srs defines the projection but does not modify the coordinates, so your output vrt file is likely to be incorrect (based on a sphere and not an ellipsoid, longitude of origin is 100 and intersection at parallel -70 instead of -60). Below is the proj4 string for EPSG:3412
+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs
If you are sure that your input data is in NSIDC Sea Ice Polar Stereographic South, you don't need the first line of code and you could write :
EDIT: 1) in the proj4 string, you can read the parameters of the projection:
lat_0 is the latitude of origin, that is the summit of the cone in case of stereographic projection
lat_ts is the latitude of the parallel where the cone intersects the ellipsoid
lon_0 is the longitude of the origin of the projection
k is a scale factor
x_0 and y_0 are the shifts of the origin of the cartesian coordinate system (it is used to avoid negative coordinates with some projection)
a and b are the size (in m) of the semi-major axes of the ellipsoid (when they are the same, you have a sphere)
2) Your output pixel size will be in degrees and you can set it using -tr. An average pixel size could be your original pixel size in meter divided by 110000 (this is roughly the size in meter of one degree along a meridian). However, stereographic projection and WGS84 coordinates are very different. Therefore there will therefore be a large amount of resampling due to distortions, and you might have to optimize your pixel size for the location where your work is focused.