[gdal-dev] ogr2ogr: issues going from polar2geographic

[Even Rouault]

Le vendredi 08 juillet 2016 21:05:05, Robb Wright a écrit :
> I'm having problems going from a Stereographic North Pole to 4326.
> Contiguous polys across the 180 line in polar are splitting and wrapping
> across the 4326 space, instead of just splitting at the 180.  The data
> doesn't go south of ~50 degrees latitude.
> 
> ogr2ogr --config CHECK_WITH_INVERT_PROJ=YES "ESRI Shapefile"
> output_wgs84.shp masie_ice_r00_v01_2016187_4km.shp -s_srs "+proj=stere
> +lat_0=90 +lat_ts=60 +lon_0=-80 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m
> +no_defs" -t_srs EPSG:4326
> 
> I've tried it with/without wrapdateline, with/without
> CHECK_WITH_INVERT_PROJ but all output looks the same.
> The proj4 came from gdalsrsinfo on the source file.
> 
> Source shapefile:
> ftp://sidads.colorado.edu/DATASETS/NOAA/G02186/latest/4km/maise_ice_r00_4km
> .zip
> 
> Using GDAL 2.1.0, released 2016/04/25 from the built GISInternals site.
> 
> If anyone has thoughts on how to go from a polar to anything else, it
> would be appreciated.

Robb,

ogr2ogr -wrapdateline heuristics aren't ready to deal with such complex cases
as polygons in polar projection that contain the pole itself and/or that cross the
merdian +/- 180 deg (here we have both)

I've finally found a method in 4 steps :

1) create needle.csv with the following content :

id,WKT
1,"POLYGON((179.99999 0,179.99999 89.99999,-179.99999 89.99999,-179.99999 0,179.99999 0))"

2) project it to the polar system :

$ ogr2ogr needle_polar.shp needle.csv -s_srs EPSG:4326 -t_srs "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -overwrite

3) erase the needle from the original shapefile :

$ python ogr_layer_algebra.py Erase -input_ds masie_ice_r00_v01_2016189_4km.shp -method_ds needle_polar.shp  -output_ds polar_without_pole.shp -overwrite

where ogr_layer_algebra.py comes from
https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ogr_layer_algebra.py

This step makes sure that no polygons include a point at longitude +/- 180 deg, as well as
making sure that the ~ 90 latitude appears in the coordinates of the polygon that contains the pole.

4) Finally project it to EPSG:4326

$ ogr2ogr result.shp polar_without_pole.shp -t_srs EPSG:4326 -overwrite


Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com