Map Reprojection: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(minor cleanup)
(cleanup with Reprojection content)
Line 4: Line 4:




'''Rationale''': Reprojection on the fly can easily introduce artifacts. Especially,
=== 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 {{cmd|v.split}} be used;
* vector maps: lines and polygons need to have a sufficiently high amount of vertices which is not guaranteed, hence {{cmd|v.split}} be used;
* raster maps:
* raster maps:
Line 12: Line 14:




'''Reprojecting data in GRASS GIS:'''
=== Reprojecting data in GRASS GIS ===


# first create an location in the source projection/geodetic datum (using the [[GRASS Location Wizard]] d you can generate the location from the dataset itself)
# first create an location in the source projection/geodetic datum (using the [[GRASS Location Wizard]] d you can generate the location from the dataset itself)
Line 23: Line 25:
Another (similar) strategy frequently needed when maps need to be imported using, e.g. {{cmd|v.in.ogr}} or {{cmd|r.in.gdal}}, is to tell these programs to create a new location (argument "location"). If the imported data contains proper <code>PROJ_INFO</code> files, they can be reprojected into the current mapset, using the the projection in the current location/mapset.
Another (similar) strategy frequently needed when maps need to be imported using, e.g. {{cmd|v.in.ogr}} or {{cmd|r.in.gdal}}, is to tell these programs to create a new location (argument "location"). If the imported data contains proper <code>PROJ_INFO</code> files, they can be reprojected into the current mapset, using the the projection in the current location/mapset.


=== Raster data: GDAL alternative ===


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


<source lang="bash">
<source lang="bash">
gdalwarp -t_srs '+proj=utm +zone=37 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' -r bilinear AfricaHornElev.tif AfricaHornElev37n.tif
gdalwarp -t_srs '+proj=utm +zone=37 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' -r bilinear AfricaHornElev.tif AfricaHornElev37n.tif
</source>
</source>
=== Vector data: have enough vertices in your vector data ===
When reprojecting vector data, it is important to have enough vertices in the data.
Example: Box (polygon) created with {{cmd|v.in.region}}:
g.region n=60 s=40 w=0 e=30 res=10 -p
v.in.region box_latlong_10deg
[[Image:V proj input unsplit boundaries.png|300px|center|thumb|4 corner rectangle, LatLong (no further vertices)]]<br>
v.split -n box_latlong_10deg output=box_latlong_10deg_vertices_200km length=200 units=kilometers
[[File:V proj input split boundaries.png|300px|center|thumb|Split rectangle, LatLong]]<br>
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:
# 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
# ... ops.
# 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
[[File:V proj result.png|300px|center|thumb|Results in EU LAEA (EPSG: 3035)]]


== See also ==
== See also ==

Revision as of 13:12, 16 May 2014

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

Answer: For quality reasons, GRASS GIS handles one projection per location.


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

  1. first create an location in the source projection/geodetic datum (using the GRASS Location Wizard d 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.

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.

gdalwarp -t_srs '+proj=utm +zone=37 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' -r bilinear AfricaHornElev.tif AfricaHornElev37n.tif

Vector data: have enough vertices in your vector data

When reprojecting vector data, it is important to have enough vertices in the data.

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)


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:

# 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
# ... ops.
# 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)

See also