[GRASSLIST:1544] Re: r.in.gdal and a GeoTIFF

Thu Oct 23 05:43:28 EDT 2003

On Wed, Oct 22, 2003 at 09:36:49AM -0500, Ed Davison wrote:
> On Wed, 2003-10-22 at 03:02, Markus Neteler wrote:
> > On Tue, Oct 21, 2003 at 10:46:42AM -0500, bfdi533 wrote:
> > > I am trying to load the following geotiff into GRASS with r.in.gdal:
> > > 
> ...
> > > GRASS:~/gis/data/nationalatlas/shdrlf > r.in.gdal input=shdrlfi020l.tif
> > > output=shdrlf target=tx
> > > ERROR: Projection of dataset does not appear to match current location.
> > >                                                                                                                                                                                    
> > >        cellhd.proj = 0 (unreferenced)
> > >                                                                                                                                                                                    
> > >        You can use the -o flag to r.in.gdal to override this check.
> > > 
> > > So, I tried with the -o flag and now I get:
> > > 
> > > GRASS:~/gis/data/nationalatlas/shdrlf > r.in.gdal input=shdrlfi020l.tif
> > > output=shdrlf -o
> > > WARNING: G_set_window(): Illegal latitude for North
> > 
> > Do you try to import into a Lat/Long location? What does
> > g.region -dp
> > 
> > print?
> > 
> 
> GRASS:~/gis/data/nationalatlas/shdrlf > g.region -dp
> projection: 3 (Latitude-Longitude)
> zone:       0
> datum:      wgs84
> ellipsoid:  wgs84
> north:      37N
> south:      25N
> west:       107W
> east:       93W
> nsres:      0:00:30
> ewres:      0:00:30
> rows:       1440
> cols:       1680
> 
> Yup, it looks like it.  I guess the question is, then, how do I get this into my lat/long location?

Searching around I found this file.
http://edcftp.cr.usgs.gov/pub/data/nationalatlas/shdrlfi020l.tar.gz

A good idea is always to use 'gdalinfo' to get information about Projection,
coordinate system, datum etc.:

gdalinfo shdrlfi020l.tif
Driver: GTiff/GeoTIFF
Size is 10366, 7273
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",45],
    PARAMETER["longitude_of_center",-100],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-6086629.000000,4488761.000000)
Pixel Size = (1000.000000,-1000.000000)
Metadata:
  TIFFTAG_IMAGEDESCRIPTION=Color shaded relief of North America
  TIFFTAG_SOFTWARE=USGS DRG production software
  TIFFTAG_DATETIME=2003:02:06 13:19:34
  TIFFTAG_XRESOLUTION=254
  TIFFTAG_YRESOLUTION=254
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Corner Coordinates:
Upper Left  (-6086629.000, 4488761.000) (156d 3'31.10"E, 37d45'55.65"N)
Lower Left  (-6086629.000,-2784239.000) (154d27'44.36"W,  3d 2'12.16"N)
Upper Right ( 4279371.000, 4488761.000) (  4d 0'23.08"W, 53d59'22.66"N)
Lower Right ( 4279371.000,-2784239.000) ( 61d 9'32.19"W, 11d18'12.56"N)
Center      ( -903629.000,  852261.000) (113d12'19.80"W, 51d59'8.84"N)
Band 1 Block=10366x1 Type=Byte, ColorInterp=Red
Band 2 Block=10366x1 Type=Byte, ColorInterp=Green
Band 3 Block=10366x1 Type=Byte, ColorInterp=Blue

-> This is "Lambert_Azimuthal_Equal_Area".

Method 1:
 You have to define a location accordingly to the parameters reported above.
 Then you can import into this location, exit, change to your Lat/Long
 location and run 'r.proj' there to pull the LAEA projected map into the
 Lat/Long location.

Method 2:
You run 'gdalwarp' (provided you have the full GDAL installed, at time
the CVS version or you wait for the next official release) and reproject
the GeoTIFF file to another projection, stored as new GeoTIFF file:

#get EPSG code from /usr/local/share/proj/epsg
gdalwarp -t_srs '+init=epsg:4326'  shdrlfi020l.tif shdrlfi020l_LL_wgs84.tif

#verify:
gdalinfo shdrlfi020l_LL_wgs84.tif
Driver: GTiff/GeoTIFF
Size is 20653, 4802
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-178.618743,86.014709)
Pixel Size = (0.017281,-0.017281)
Corner Coordinates:
Upper Left  (-178.6187426,  86.0147085) (178d37'7.47"W, 86d 0'52.95"N)
Lower Left  (-178.6187426,   3.0332766) (178d37'7.47"W,  3d 1'59.80"N)
Upper Right ( 178.2774494,  86.0147085) (178d16'38.82"E, 86d 0'52.95"N)
Lower Right ( 178.2774494,   3.0332766) (178d16'38.82"E,  3d 1'59.80"N)
Center      (  -0.1706466,  44.5239925) (  0d10'14.33"W, 44d31'26.37"N)
Band 1 Block=20653x1 Type=Byte, ColorInterp=Red
Band 2 Block=20653x1 Type=Byte, ColorInterp=Green
Band 3 Block=20653x1 Type=Byte, ColorInterp=Blue

This was just an example which you may adapt to your needs.

Import: I have made a script "create_location.sh" which creates a
GRASS location from a GDAL supported raster map:
Fetch from:
  http://mpa.itc.it/grasstutor/examples.phtml

Cheers

 Markus