From: Zheng-Xiang Li Date: 8. August 2013 14:42:51 MESZ To: "Georg Feulner (PIK)" Subject: Re: Map data for Li et al. (2008) Hi George, Please see attached slightly updated maps. I did not attach the 580 Ma one as I do not think it is as reliable as that of our 2008 paper. The related paper is also attached here for you. Cheers, Li > I have come across your 2008 paper in Precambrian Research, and I > wonder whether you could provide me with electronic data for the > paleogeography maps which could then be converted into topography > maps to be used in our climate mode. I would be particularly > interested in the maps in the range 720-600 Ma, and it would be > great if you could provide me with any data for those. Apparenty, the maps were exported from GPlates in GMT xy format. The general mapping tools GMT can read and display them. Export into NetCDF or some other usable-for-format was not possible. Approach: Write a small program that converts the .xy files into ESRI SHAPE files. The .xy files apparently consist of a set polygons, separated by metadata lines starting with '>' . Lines without leading '>' consist of lon/lat coordinate pairs, giving polygon vertices. Each sequence of coordinate-pair-lines is taken to be polygon of its own. If the last verte does not match the first one, i.e. the polygon appears to open, the first vertex is replicated at the end, to make it closed. The .xy file has some polygons which cross the date boundary -180 to 180, which makes very ugly artifacts in GRASS. Thus, the western hemisphere is moved from -180,0 to 180,360 . if (x < 0) x += 360; Most probably, large polygons are split up in the .xy files, and it might be possible to collect and assemble pieces with matching end/start points into larger polygons. However, this was not implemented. gcc -g -Wall -o xy2shp -lshp xy2shp.c for i in *.xy; do ./xy2shp $i `basename $i .xy`; done This creates for each input file foo.xy three output files foo.shp, foo.shx, foo.dbf, which together form the SHAPE file. Those can now be imported into GRASS. Force grass to ask for a new LOCATION and MAPSET: rm ~/.grassrc6 grass Select a grassdata directory, e.g. /home/petri/wk/grassdata Create a new LOCATION, named e.g. Neoproterozoic-0-360 For the projection, select ``arbitrary XY projection'', without any further specifications. When asked if you want to select a region, do so. In the dialogue, select -90 to 90 lat, 0 West to 360 East lon, 1 deg resolution in each direction. Then import the above generated shape file: v.in.ogr dsn=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/r_0_620.shp layer=r_0_620 output=r_0_620 Convert it to raster map: v.to.rast --verbose input=r_0_620@PERMANENT column=dummy output=rast_620 And write it out as NetCDF: r.out.gdal --v input=rast_620 \ format=netCDF \ type=Int16 \ output=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/reconstructed_0_620.00.nc # r.out.gdal names the variables always BAND1. # There are hints in the documentation that it should be possible to # set option for writing metadata (metaopt=NAME=VALUE), but I could not get that to work. v.in.ogr --verbose dsn=reconstructed_0_620.shp layer=reconstructed_0_620 output=reconstructed_0_620 v.in.ogr --verbose dsn=reconstructed_635.shp layer=reconstructed_635 output=reconstructed_635 v.in.ogr --verbose dsn=reconstructed_720.shp layer=reconstructed_720 output=reconstructed_720 v.to.rast --verbose input=reconstructed_0_620@PERMANENT column=dummy output=rast_620 v.to.rast --verbose input=reconstructed_635@PERMANENT column=dummy output=rast_635 v.to.rast --verbose input=reconstructed_720@PERMANENT column=dummy output=rast_720 r.out.gdal --v input=rast_620 format=netCDF type=Int16 output=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/reconstructed_0_620.00.nc r.out.gdal --v input=rast_635 format=netCDF type=Int16 output=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/reconstructed_635.00Ma.nc r.out.gdal --v input=rast_720 format=netCDF type=Int16 output=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/reconstructed_720.00Ma.nc Seufz. Das reconstructed_720.00Ma.xy enthaelt sowohl Teile, die ueber die Datumsgrenze gehen, als auch welche, die ueber den Null-Meridian gehen. Da muessen wir eine ganz eigene Rotation finden: Neue Option fuer xy2shp: ./xy2shp -10 reconstructed_720.00Ma.xy reconstructed_720 shifts all longitudes by 10 deg Eastwards. Then in grass, as above: v.in.ogr --verbose dsn=reconstructed_720.shp layer=reconstructed_720 output=reconstructed_720 v.to.rast --verbose input=reconstructed_720@PERMANENT column=dummy output=rast_720 r.out.gdal --v input=rast_720 format=netCDF type=Int16 output=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/reconstructed_720.00Ma-shifted10.nc Now the .nc File needs lon-rotation Westwars by 10 degrees. Do that with ncap2: ncap2 -O -v -c -s '*tmp=Band1;Band1(:,0:349)=Band1(:,10:359);Band1(:,350:359)=tmp(:,0:9);' \ reconstructed_720.00Ma-shifted10.nc reconstructed_720.00Ma.nc #Aber: in den Files steht unsinniges Zeugs. # #> name ISLANDS AROUND BAFFIN ISLAND, NORTH AMERICA #> oldPlatesHeader OldPlatesHeader 9922 67 ISLANDS AROUND BAFFIN ISLAND, #NORTH AMERICA 107 1750 -999 CS 1 107 1 # -143.608643 -19.306530 # -143.472009 -19.392650 # -143.608643 -19.306530 # #Das ist einfach ein geschlossenes Polygon aus zwei Punkten. Sprich: eine Linie. Von solchen Teiln wimmelt es. # #> name Caledonides of Grenland #> oldPlatesHeader OldPlatesHeader 9010 910 Caledonides of Grenland 910 #3999 -999 DS 5 910 1 # -139.050656 -33.716610 # -138.487207 -33.788819 #... # -129.322931 -23.422659 # -129.816236 -23.587934 # #Das ist ein einzelnes offenes Polygon, das aber offenbar die kleinen Inseln um Groenland herum darstellen soll. #Um eine tektonische Platte zu modellieren mag das OK sein, aber als Land-See-Maske taugt das nicht. # Auch Polygonzuege mit einem einzigen Vertex kommen vor. # Vermutlich kann man die Teilen anhand der #> reconstructionPlateId 108 # irgendwie wieder zusammensetzen. # O.a. Id 108 z.B. kommt 11 mal in reconstructed_0_620.00.xy vor. Aber # einige Teile davon sind in sich schon abgeschlossen, andere haben # gemeinsame Anfangspunkte, einmal finde ich den Endpunkt eines Teils # als Anfangspunkt eines anderen Teils. # Da finde ich noch kein Sortierprinzip drin. # Du solltest ihn mal fragen, ob er die Daten nicht auch als # Shape-format schicken kann? Lt. Online-Manual sollte GPlates das # exportieren koennen. ----- From: Zheng-Xiang Li Date: 15. August 2013 09:49:43 MESZ To: "Georg Feulner (PIK)" Subject: Re: Map data for Li et al. (2008) Hi George, Hope you will have better luck with this format _ let me know. attached: Li et al 2013_720 Ma shape filkes.zip Those files contain a database column named ``TIME''. GRASS does not like that, it must be renamed upon import, thus the ``cnames='' list. use dbfinfo filename.dbf to get a list of the columns in the file. v.in.ogr dsn=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/Li_et_al_2013_720_Ma_shapefiles/reconstructed_720.00Ma/reconstructed_720.00Ma_polyline.shp output=rec_reconstructed_720_00Ma_poly min_area=0.0001 snap=-1 cnames=cat,ANCHOR,TIME_,FILE1,FILE2,FILE3,PLATE_ID,TYPE,GPGIM_TYPE,FROMAGE,TOAGE,NAME,DESCR,FEATURE_ID,PLATEID2,RECON_METH,L_PLATE,R_PLATE Number of nodes: 481 Number of primitives: 345 Number of points: 0 Number of lines: 345 Number of boundaries: 0 Number of centroids: 0 Number of areas: 0 Number of isles: 0 v.to.rast --v input=rec_reconstructed_720_00Ma_poly column=cat output=rec_rast_720_poly v.in.ogr dsn=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/Li_et_al_2013_720_Ma_shapefiles/WA_lines/reconstructed_720.00Ma.shp output=WA_reconstructed_720_00Ma min_area=0.0001 snap=-1 cnames=cat,ANCHOR,TIME_,FILE1,FILE2,FILE3,CODE,NAME,MAXAGE,MINAGE,COMMENT,REFERENCE,GEON,SOURCETHM,LENGTH Number of nodes: 8 Number of primitives: 5 Number of points: 0 Number of lines: 5 Number of boundaries: 0 Number of centroids: 0 Number of areas: 0 Number of isles: 0 v.to.rast --v input=WA_reconstructed_720_00Ma column=cat output=WA_rast_720 v.in.ogr dsn=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/Li_et_al_2013_720_Ma_shapefiles/schina_lin/reconstructed_720.00Ma.shp output=shin_reconstructed_720_00Ma min_area=0.0001 snap=-1 cnames=cat,ANCHOR,TIME_,FILE1,FILE2,FILE3,CODE,NAME,MAXAGE,MINAGE,COMMENT,REFERENCE,GEON,SOURCETHM,LENGTH,PLATEID,FROMAGE,TOAGE,DESCR,TYPE Number of nodes: 135 Number of primitives: 94 Number of points: 0 Number of lines: 94 Number of boundaries: 0 Number of centroids: 0 Number of areas: 0 Number of isles: 0 v.to.rast --v input=shin_reconstructed_720_00Ma column=cat output=shin_rast_720 v.in.ogr dsn=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/Li_et_al_2013_720_Ma_shapefiles/li1\ data\ without\ poles/reconstructed_720.00Ma/reconstructed_720.00Ma_polyline.shp output=li_R_reconstructed_720_00Ma_poly min_area=0.0001 snap=-1 cnames=cat,ANCHOR,TIME_,FILE1,FILE2,FILE3,PLATE_ID,TYPE,GPGIM_TYPE,FROMAGE,TOAGE,NAME,DESCR,FEATURE_ID,PLATEID2,RECON_METH,L_PLATE,R_PLATE Number of nodes: 338 Number of primitives: 246 Number of points: 0 Number of lines: 246 Number of boundaries: 0 Number of centroids: 0 Number of areas: 0 Number of isles: 0 v.to.rast --v input=li_R_reconstructed_720_00Ma_poly column=cat output=li_R_rast_720_poly The files in the subdirectories WA_lines/ and schina_lin/ seem to be auxiliary tectonic lines. Subdirectory reconstructed_720.00Ma/ seems to be merge of subdirectories WA_lines/, schina_lin/, and li1\ data\ without\ poles/reconstructed_720.00Ma/ Thus the files in subdirectory li1\ data\ without\ poles/reconstructed_720.00Ma/ seem most apropirate for our purposes. Those shape files apparently contain (open) polyline definitions, not (closed) polygon definitions. Thus, GRASS is not able to determine any inside, outside or area of the shapes, and converting to raster only produces outlines, but without any further meaning for GRASS. So, we can try to fill in the land areas by hand. r.out.gdal --v input=li_R_rast_720_poly format=netCDF type=Int16 output=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/reconstructed_720.00-outline.nc r.out.gdal --v input=li_R_rast_720_poly format=PNG output=/iplex/01/climber3/petri/mom5.0.2/exp/Neoproterozoic/reconstructed_720.00-outline.png ---------------- 19-Aug-2013 Georg used GIMP to fill the outlines in the .png version. Stefan used gdal to convert it to NetCDF -------------- Wed May 31 12:30:37 CEST 2017 For visualisation, Stefan Fuchs wants to have the continental outlines of the 720Ma map in resolution 2400x1800 Lets try that with shapefiles and GRASS. grass In the dialog box, select ``New'' GRASS Location. Name it Neoproterozoic-0-360-highres ``Choose method for creating a new location'' -> Create a generic Cartesian coordinate system (XY) ``Do you want to set the default region extents and resolution now?'' -> Yes In the dialogue, select -90 to 90 lat, 0 West to 360 East lon, For E-W resolution we get 360/2400 = 0.15 For N-S resolution we get 180/1800 = 0.10 Next dialog box asks ``Do you want to create a new mapset?'' proposal is ``petri''. I have no idea what this is good for. The default mapset ``PERMANENT'' should be OK for us. Thus klick ``cancel'' here. Then ``Start GRASS session''. First attempt of using v.in.ogr complained: ERROR: Projection of dataset does not appear to match current location. Import dataset PROJ_INFO is: name: GCS_WGS_1984 datum: wgs84 ellps: wgs84 proj: ll no_defs: defined In case of no significant differences in the projection definitions, use the -o flag to ignore them and use current location definition. Consider generating a new location with 'location' parameter from input data set. We just use -o to override the projection information: v.in.ogr --verbose -o input=/p/projects/climber3/petri/POEM/exp/Neoproterozoic/Li_et_al_2013_720_Ma_shapefiles/li1\ data\ without\ poles/reconstructed_720.00Ma/reconstructed_720.00Ma_polyline.shp output=li_R_reconstructed_720_00Ma_poly min_area=0.0001 snap=-1 columns=cat\,ANCHOR\,TIME_\,FILE1\,FILE2\,FILE3\,PLATE_ID\,TYPE\,GPGIM_TYPE\,FROMAGE\,TOAGE\,NAME\,DESCR\,FEATURE_ID\,PLATEID2\,RECON_METH\,L_PLATE\,R_PLATE Over-riding projection check Check if OGR layer contains polygons... 100% Using native format Column name