[gdal-dev] ogr2ogr: issues going from polar2geographic

[Even Rouault]

Le samedi 09 juillet 2016 01:00:46, Even Rouault a écrit :
> 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_4
> > km .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
> 

Update: I've implemented recently in trunk the above trick, so ogr2ogr will 
now work correctly out of the box for polar projection (and projections 
spanning the antimeridian) to EPSG:4326

Even

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