Sunday, 8 July 2018

ogr - Changing shapefile longitude from -180 - 180 to 0 - 360?


I have a shapefile of world countries, with longitude between -180° and 180°. How can I change the shapefile (eg. with ogr2ogr) for having longitude between 0 and 360°.



Following OSGEO/Trac explanation about GenParms, I tried this without success (the result is still in -180 - 180:



ogr2ogr -wrapdateline -t_srs '+proj=latlong +datum=WGS84 +lon_wrap=-180 +over' /data/tmp/new_algo/world_new.shp /data/tmp/new_algo/world.shp



I need to have the new longitude defined as:



if Lon < 0: new_Lon = 360 + Lon else: new_Lon = Lon



Of course, the process must take care of the date line.



Answer




Using GDAL >= 1.10.0 compiled with SQLite and SpatiaLite:


ogr2ogr world_shifted.shp world.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,180,0) FROM world"

or:


ogr2ogr -s_srs EPSG:4326 -t_srs "+proj=longlat +ellps=WGS84 +pm=-180 +datum=WGS84 +no_defs" world_shifted.shp world.shp

Both commands produce a longitude offset of 180°, i.e. a prime meridian of -180° is considered. In fact:


>ogrinfo world_shifted.shp world_shifted | grep Extent
Extent: (0.000000, -90.000000) - (360.000000, 83.623596)


The difference between the two commands is that with a longitude offset (2nd try) data are simply reprojected using -180° as prime meridian, while shifting the coordinates geometries (1st try) are altered, even if the result is apparently the same.


EDIT
If there are parts in 0-180 that should not move, it's possible to adapt this working solution: https://gis.stackexchange.com/a/73164/22405


Clip the two parts:


ogr2ogr world_part1.shp world.shp -clipsrc -180 -90 0 90
ogr2ogr world_part2.shp world.shp -clipsrc 0 -90 180 90

Shift only the first part:


ogr2ogr world_part1_shifted.shp world_part1.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,360,0), CNTRY_NAME FROM world_part1"


Then, merge the second part and the first shifted:


ogr2ogr world_0_360_raw.shp world_part2.shp
ogr2ogr -update -append world_0_360_raw.shp world_part1_shifted.shp -nln world_0_360_raw

Finally, dissolve countries boundaries of world_0_360_raw.shp obtaining world_0_360.shp by country names. For instance:


ogr2ogr world_0_360.shp world_0_360_raw.shp -dialect sqlite -sql "SELECT ST_Union(Geometry), CNTRY_NAME FROM world_0_360_raw GROUP BY CNTRY_NAME"

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...