4

I simply want to project a polar stereographic raster (EPSG:3411) to Mercator (EPSG:3857) or WGS84 (EPSG:4326). The data I is from NOAA.

$gdalwarp -s_srs EPSG:3411 -t_srs EPSG:3857 northpsg.20141027 output.tif

The resulting raster is croped near the height of 75°N - which is exactly the opposite behavior of projecting into polar stereographic that zetah has

sceenshot

On the other hand the projection into WGS84 results in a output with a height of 0.

Do you have any suggestions how this could work?

I already contacted the GDAL community but didn't get any responses yet.

2

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): enter image description here


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.

enter image description here

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:

enter image description here

You find the reference point as a little yellow dot to the south of Australia.

  • Thank you so much, that helps me out! This complexity is unbelievable! I tried something like this with gdal_edit.py before. Unfortunately the transfer from the wgrib parameters to the new ullr-parameter is not really clear for me. How did you calculated them? In my case I am not able to project a different file with the following parameters $ wgrib polar stereo: Lat1 -41.501000 Long1 -135.000000 Orient 0.000000 north pole (790 x 830) Dx 10000 Dy 10000 scan 64 mode 0 $proj4-string +proj=stere +lat_0=90 +lon_0=0 +k=0.9330127018922193 +x_0=0 +y_0=0 +a=6371229 +b=6371229 +units=km +no_defs – dnltsk Oct 29 '14 at 15:46
  • I added a few more steps from my workaround. It is a bit trial and error. Looking at the map in polar stereo helps to see if it fits to Natural Earth vectors or not. Note that your new projection has units=km. – AndreJ Oct 29 '14 at 16:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.