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:
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 -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
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:
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?