GRASS GIS 7 just got better: When reprojecting vector data, now automated vertex densification is applied. This reduces the reprojection error for long lines (or polygon boundaries). The needed improvement has been kindly added in v.proj by Markus Metz.

Example

As an (extreme?) example, we generate a box in LatLong/WGS84 (EPSG: 4326) which is of 10 degree side length (see below for screenshot and at bottom for SHAPE file download of this “box” map):

[neteler@oboe ~]$ grass70 ~/grassdata/latlong/grass7/
# for the ease of generating the box, set computational region:
g.region n=60 s=40 w=0 e=30 res=10 -p
projection: 3 (Latitude-Longitude)
zone:       0
datum:      wgs84
ellipsoid:  wgs84
north:      60N
south:      40N
west:       0
east:       30E
nsres:      10
ewres:      10
rows:       2
cols:       3
cells:      6
# generate the box according to current computational region:
v.in.region box_latlong_10deg
exit

Next we start GRASS GIS in a metric projection, here the EU LAEA:

# EPSG 3035, metric EU LAEA:
grass70 ~/grassdata/europe_laea/user1/
GRASS 7.0.0svn (europe_laea): >

Now we first reproject the map the “traditional way” (no vertex densification as in most GIS, here enforced by smax=0):

v.proj box_latlong_10deg out=box_latlong_10deg_no_densification
location=latlong mapset=grass7 smax=0

Then we do a second reprojection with new automated vertex densification (here we use the default values for smax which is a 10km vertex distance in the reprojected map by default):

v.proj box_latlong_10deg out=box_latlong_10deg_yes_densification
location=latlong mapset=grass7

Eventually we can compare both reprojected maps:

g.region vect=box_latlong_10deg_no_densification

# compare:
d.mon wx0
d.vect box_latlong_10deg_no_densification color=red
d.vect box_latlong_10deg_yes_densification color=green fill_color=none

Comparison of the reprojection of a 10 degree large LatLong box to the metric EU LAEA (EPSG 3035): before in red and new in green. The grid is based on WGS84 at 5 degree spacing.

The result shows how nicely the projection is now performed in GRASS GIS 7: the red line output is old, the green color line is the new output (its area filling is not shown).

Consider to benchmark this with other GIS… the reprojected map should not become a simple trapezoid.

Sample dataset download

Download of box_latlong_10deg.shp for own tests (1kB).

Get vector map extent

You can easily grep the map extent of a vector map (bounding box):

ogrinfo /cdrom/ITALY_GC.SHP ITALY_GC | grep Extent
Extent: (9.299460, 43.787362) – (13.911350, 47.070188)

Merge of two SHAPE files

Merge of two SHAPE files ‘file1.shp’ and ‘file2.shp’ into a new file ‘file_merged.shp’ is performed like this:

ogr2ogr file_merged.shp file1.shp
ogr2ogr -update -append file_merged.shp file2.shp -nln file_merged file2

The second command is opening file_merged.shp in update mode, and trying to find existing layers and append the features being copied. The -nln option sets the name of the layer to be copied to.

Vector map reprojection

We reproject from the source projection (as defined in .prj file) to WGS84/LL:

ogr2ogr vmap0rd_ll.shp -t_srs “EPSG:4326” vmap0rd.shp

If the .prj file is missing, you can use the ‘epsg_tr.py’ utility to create it if you know the EPSG code:

epsg_tr.py -wkt 4326 > cities.prj

Reproject to current GRASS location projection:

ogr2ogr -t_srs “`g.proj -wf`” polbnda_italy_GB_ovest.shp polbnda_italy_LL.shp

Cut out a piece of a vector map
Use spatial query extents: -spat xmin ymin xmax ymax (W S E N)

ogr2ogr ARC_BZ.shp -spat 10 45 13 47 ARC.shp

Get VMAP0 metadata info:

ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr | grep bnd
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘roadl@trans(*)_line’

MAP0: Extract spatial subregion, reproject from NAD83 to WGS84

# coordinate order: W S E N
ogr2ogr -spat 19.95035 -26.94755 29.42989 -17.72624 -t_srs ‘EPSG:4326’ \
polbnda_botswana.shp gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr \
‘polbnda@bnd(*)_area’

OGR and SQL

Sample ‘where’ statements (use -sql for PostgreSQL driver):

# -where ‘fac_id in (195,196)’
# -where ‘fac_id = 195’
ogrinfo -ro -where ‘fac_id in (195,196)’ \
gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’

VMAP0 examples

Find out the Countries VMAP0 coding:

ogdi_info -u gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr \
-l ‘polbnda@bnd(*)_area’ -f area | grep Botswana

or read the VMAP0 Military specs, page 75

Extract Botswana, reproject on the fly from NAD83 to WGS84, store as SHAPE:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘BC'” polbnda_botswana.shp \
gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’

Extract Germany, reproject on the fly from NAD83 to WGS84, store as SHAPE:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘GM'” polbnda_germany.shp
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘polbnda@bnd(*)_area’
ogrinfo -summary polbnda_germany.shp polbnda_germany | grep Extent
# Extent: (5.865639, 47.275776) – (15.039889, 55.055637)
# W S E N

VMAP0 Contour lines for Germany:

ogr2ogr -t_srs “EPSG:4326” -spat 5.865639 47.275776 15.039889 55.055637 \
contour_lines.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘contourl@elev(*)_line’

VMAP0 elevation spots (points) for Germany:

ogr2ogr -t_srs “EPSG:4326” -spat 5.865639 47.275776 15.039889 55.055637 \
elevation_spots.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘elevp@elev(*)_point’

VMAP0 lakes of Trentino province in Italy:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘IT'” \
-spat 10.340029 45.261888 10.98727 45.98993 \
lakes_italy.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘inwatera@hydro(*)_area’

Connect OGR and PostgreSQL/PostGIS

ogrinfo PG:’host=grass.itc.it user=postgres dbname=ogc_simple’
ogr2ogr out.shape PG:’host=grass.itc.it user=postgres dbname=ogc_simple’ lake_geom

GRASS 6 and OGR

Convert GRASS 6 vector map to SHAPE (needs GDAL-OGR-GRASS plugin):

# -nln is “new layer name” for the result:
ogr2ogr archsites.shp grassdata/spearfish60/PERMANENT/vector/archsites/head 1 \
-nln archsites

Using WKT files with ogr2ogr

The definition is in ESRI WKT format. If you save it to a text file called out.wkt you can do the following in a translation to reproject input latlong points to this coordinate system:

ogr2ogr -s_srs WGS84 -t_srs ESRI::out.wkt out_dir indatasource

Most comand line options for GDAL/OGR tools that accept a coordinate system will allow you to give the name of a file containing WKT. And if you prefix the filename with ESRI:: the library will interprete the WKT as being ESRI WKT and convert to “standard” format accordingly. The -s_srs switch is assigning a source coordinate system to your input data (in case it didn’t have this properly defined already), and the -t_srs is defining a target coordinate system to reproject to.

TIGER files in OGR

# linear features:
ogr2ogr tiger_lines.shp tgr46081.rt1 CompleteChain

# area features:
export PYTHONPATH=/usr/local/lib/python2.5/site-packages
tigerpoly.py tgr46081.rt1 tiger_area.shp

OGR CSV driver: easily indicate column types

You can now write a little csv help file to indicate the columns types to OGR. It works as follows. Suppose you have a foobar.csv file that looks like this:

“ID”,”X”,”Y”,”AREA”,”NAME”
“1”,”1023.5″,”243.56″,”675″,”FOOBAR”
…

Now write a foobar.csvt file like this one:

“Integer”,”Real”,”Real”,”Integer”,”String”

The driver will then use the types you specified for the csv columns. The types recognized are Integer, Real and String, DateTime, and Date.

Convert KML to CSV (WKT)

First find layers:

ogrinfo -so myfile.kml

Then convert KML to CSV:

ogr2ogr -f CSV out.csv myfile.kml -sql “select *,OGR_GEOM_WKT from myfilelayer”
cat out.csv

Or use the cool online converter: https://geoconverter.hsr.ch

GRASS GIS 7: Vector data reprojection with automated vertex densification