Map Reprojection

Question: How to change map/image projections, datums, etc in GRASS GIS?

Answer: For quality reasons, GRASS GIS handles one projection per location. You can reproject between locations as well as reproject during data import if needed.

Rationale

Reprojection on the fly can easily introduce artifacts. Especially,

• vector maps: lines and polygons need to have a sufficiently high amount of vertices which is not guaranteed, hence v.split be used;
• raster maps:
  • categorical maps need to be reprojected with nearest-neighbor method, while
  • floating point data might be better reprojected with bilinear or cubic convolution methods)

Since many GIS users tend to blindly go ahead without selecting proper methods, the simple approach was adopted: put data which are in different projections into separate locations. That's it.

Reprojecting data in GRASS GIS

Variant 1:

1. first create a location from the dataset itself (e.g. using the GRASS Location Wizard)
2. second: use r.import or v.import

Ready.

Variant 2:

1. first create a location in the source projection/geodetic datum (using the GRASS Location Wizard; note that you can generate the location from the dataset itself)
2. import the map/image into the location,
3. create a second destination location in the projection/geodetic datum you want to reproject the map/image into,
4. Within the destination location session, use r.proj or v.proj (depending on whether the map/image is raster or vector) to reproject the map/image from the import location to the destination location.

Ready.

Variant 3:
Another (similar) strategy frequently needed when maps need to be imported using, e.g. v.in.ogr or r.in.gdal, is to tell these programs to create a new location (argument "location"). If the imported data contains proper PROJ_INFO files, they can be reprojected into the current mapset, using the the projection in the current location/mapset.

Raster data: GDAL alternative

Export as a GeoTIFF, use gdalwarp to change the projection and then import the new GeoTIFF file into GRASS GIS. If you have installed gdal, you almost certainly have gdalwarp as well. Indications for gdalwarp are found at http://www.gdal.org/gdalwarp.html. The options for gdalwarp may be a bit confusing for newbies. In the following example of projecting a GeoTIFF based on SRTM elevation data to UTM 37N, -t_srs is the output file projection, AfricaHornElev.tif is the input file and AfricaHornElev37n.tif is the output file.

# use bilinear resampling (-r; use bilinear for floating point maps)
# use 2GB of RAM as cache
# target resolution (-tr): ensure 100m pixel resolution
gdalwarp --config GDAL_CACHEMAX 2000 -t_srs '+proj=utm +zone=37 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' \
         -r bilinear -tr 100 100 AfricaHornElev.tif AfricaHornElev37n.tif

Extracting spatial subset (subregion)

If you just want to reproject a part of the image, you may extract a spatial subset first

# projwin order: W N E S
gdal_translate -of GTiff -projwin 636861 5152686 745617 5054047.5 \
     p192r28_5t19920809_nn1.tif test1_utm.tif

Vector data: have enough vertices in your vector data

When reprojecting vector data, it is important to have enough vertices in the data. Starting from GRASS GIS 7 the vertices are automatically densified.

Example: Box (polygon) created with v.in.region:

g.region n=60 s=40 w=0 e=30 res=10 -p
v.in.region box_latlong_10deg

4 corner rectangle, LatLong (no further vertices)

Update 12/2014: vertex densification is now automatically done in v.proj of GRASS GIS 7! Just use v.proj and it will be right.

v.split -n box_latlong_10deg output=box_latlong_10deg_vertices_200km length=200 units=kilometers

Split rectangle, LatLong

As obvious from the results, enough vertices must be present in the map to avoid (severe) distortion. Note that this applies likewise to "ogr2ogr" and QGIS:

# switch to a different location.
# Properly reproject in GRASS GIS (when having used v.split):
grass70 -c box_latlong_10deg_EU_LAEA.shp ~/grassdata/eulaea
v.proj box_latlong_10deg location=latlong_wgs84 mapset=neteler
v.proj box_latlong_10deg_vertices_200km location=latlong_wgs84 mapset=neteler

Results in EU LAEA (EPSG: 3035)

# Test in OGR/QGIS
v.out.ogr box_latlong_10deg dsn=box_latlong_10deg.shp
ogr2ogr -t_srs epsg:3035 box_latlong_10deg_EU_LAEA.shp box_latlong_10deg.shp
qgis box_latlong_10deg_EU_LAEA.shp
# ... oops... (result not shown here).

See also

• Location and Mapsets
• WxGUI
• Projections

Other publications and sites:

• CRS Explorer - PROJ codes
• Interactive Album of Map Projections 2.0
• Compare Map Projections: View Map Projections