2

Yes, it's yet another transformation question.

I am working with ASTER DEM data, which is LatLong WGS84. I want to process the DEM map into an STL file for use in a CFD simulation. The script I am using to accomplish this requires metric coordinates, though.

Here's what my ASTER DEM mosaic looks like:

enter image description here

I have created a new location in GRASS and chosen UTM with the etrs89 datum and grs80 spheroid. (I don't know if this is the correct choice.) Following these instructions, I created a mosaic of some ASTER tiles using gdalwarp:

gdalwarp -t_srs "+proj=utm +zone=32 +datum=etrs89" -r bilinear -dstnodata -9999 -co "TILED=YES" AST_DEM.vrt AST_DEM_Rhein_Main.tif

and managed to successfully import them.

(EDIT: It would appear that proj4 does not support ETRS89 as a datum:

$ proj -ld
__datum_id__ __ellipse___ __definition/comments______________________________
       WGS84 WGS84        towgs84=0,0,0
      GGRS87 GRS80        towgs84=-199.87,74.79,246.62
                          Greek_Geodetic_Reference_System_1987
       NAD83 GRS80        towgs84=0,0,0
                          North_American_Datum_1983
       NAD27 clrk66       nadgrids=conus,ntv1_can.dat
                          North_American_Datum_1927
     potsdam bessel       towgs84=606.0,23.0,413.0
                          Potsdam Rauenberg 1950 DHDN
    carthage clark80      towgs84=-263.0,6.0,431.0
                          Carthage 1934 Tunisia
hermannskogel bessel       towgs84=653.0,-212.0,449.0
                          Hermannskogel
       ire65 modif_airy   towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15
                          Ireland 1965

Apparently proj doesn't fail with an error if it doesn't recognize a parameter, it just silently proceeds without warning the user, even if the result is inaccurate. Caveat emptor!)

I have a couple of problems with the result. First, when I imported the resulting mosaic, GRASS complained that it didn't recognized the datum:

r.in.gdal input=/home/sfbosch/Documents/Quelle/OpenFOAM/GIS/ASTER/Data/AST_DEM_Rhein_Main.tif output=AST_DEM_Rhein_Main --overwrite
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Projection of input dataset and current location appear to match
r.in.gdal complete. Raster map <AST_DEM_Rhein_Main> created.

Next, the mosaic obviously looks quite, well, warped:

UTM

I expected some change with the reprojection, but nothing that stark; could this be due to the non-existent (in proj4) datum? I am having trouble recognizing the topographical landmarks I had picked out for myself. Finally, I had been using OpenStreetMap data, and now it won't import because the projections don't match:

ERROR: Projection of dataset does not appear to match current location.

GRASS LOCATION PROJ_INFO is:
name: Universal Transverse Mercator
proj: utm
zone: 32
no_defs: defined
datum: etrs89
ellps: grs80
towgs84: 0.000,0.000,0.000

Import dataset PROJ_INFO is:
name: Lat/Lon
proj: ll
datum: wgs84
ellps: wgs84
no_defs: defined

You can use the -o flag to v.in.ogr to override this projection check.
Consider generating a new location with 'location' parameter from input data set.

I tried it with the override, but ended up with a layer that wasn't visible. I've since learned that I also need to reproject the vector layers with v.proj.

Really, all I want is to work with metric coordinates; it almost doesn't matter which system, the topography is the important part. Was there something wrong with my gdalwarp approach? Am I making this too complicated? What's the best way to accomplish what I want?

  • 1
    The error message seems to stem from "v.in.ogr" but you want to import with r.in.gdal... Please update your question. – markusN Sep 2 '13 at 20:13
  • @markusN I am reprojecting the raster data with r.in.gdal, but for the vector data from OpenStreetMap I attempted to do a plain import with v.in.ogr, which doesn't work. The vector map also needs to be reprojected. – Stephen Bosch Sep 2 '13 at 21:33
  • 1
    I think you should really use the proper EPSG code: spatialreference.org/ref/epsg/25832 – markusN Sep 5 '13 at 4:05
  • 1
    Stephen, just another small note: r.in.gdal does not re-project data. It simply imports data. The Location in which the data will be imported (well, the CRS upon which the Location is defined) should be identical to the CRS in which the data are referenced. To re-project data into another reference system, one needs to create & enter another new Location (based on the desired target CRS) and use r.proj for the reprojection. – Nikos Alexandris Oct 13 '13 at 11:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.