[GRASS-user] [gdal-dev] Working with NTF files

[Nikos Alexandris]

Markus Metz:

> >> This is mixing different concepts, i.e. import and georeferencing. You
> >> can use i.rectify after importing raster data. RPC is not supported by
> >> GRASS, thus the gdalwarp approach is most probably the preferred
> >> solution. In good GRASS tradition, you could also figure out a
> >> reasonable target grid geometry and set this with the corresponding
> >> option to gdalwarp.

Nikos Alexandris:

> > Ok -- On/Off: would it make sense to be able to say
> > grass -c Some_NITF_Container /grassdb/Location_Based_On_NITF

Markus Metz wrote:
 
> You probably mean
> grass -c Some_GDAL/OGR_recognized_Container
> /grassdb/Location_Based_On_GDAL/OGR_recognized_dataset

Absolutely.
 
> which is already possible.

Is it (always)?  Please, see below whenever time permits.

 
> > and, if need be, select the SDS that contains/is the SRS, i.e.
> > grass -c Some_NITF_Container /grassdb/Location_Based_On_NITF sds=2
 
> No, for e.g. hdf and NITF you can specify the sds directly, no need to
> have a sds= option.

 
> In general, grass -c <geofile> would use the SRS of the geodataset,

which is?  The one that gdalinfo reports on top, right?


> not the SRS of the embedded GCPs or any other geotransform
> information. Transformation/georeferencing would be done after import
> with i.rectify or before import with gdalwarp.

Got it. I, think, I have a problem then...

Working with

# GDAL
gdal-config --version

1.10.0

# which?
which gdal-config

/usr/bin/gdal-config


# formats:
gdalinfo --formats | grep NI

  NITF (rw+vs): National Imagery Transmission Format

# grass?
g.version -bg

version=7.0.svn
date=2013
revision=57074
build_date=2013-06-12

 ./configure  --with-cxx --with-includes=/usr/include/ --with-libs=/usr/lib64/ 
--with-proj --with-proj-includes=/usr/include/ --with-proj-libs=/usr/lib64/ --
with-proj-share=/usr/share/proj/ --with-geos --with-geos=/usr/bin/geos-config 
--with-gdal=/usr/bin/gdal-config --with-x --with-motif --with-cairo --with-
opengl-libs=/usr/include/GL --without-ffmpeg --with-python=yes --with-
python=/usr/bin/python2.7-config --with-wxwidgets --with-freetype=yes --with-
freetype-includes=/usr/include/freetype2/ --with-odbc=yes --with-sqlite=yes --
with-mysql=yes --with-mysql-includes=/usr/include/mysql --with-mysql-
libs=/usr/lib/mysql --with-postgres=yes --with-postgresql=yes --with-postgres-
includes=/usr/include/postgresql --with-opencl --with-openmp --with-pthread --
with-lapack --with-fftw --with-readline --with-regex --with-nls --with-jpeg --
with-tiff --with-png --with-netcdf --without-opendwg --enable-largefile=yes


Consider the following example

# below, obviously not epsg=4326, but UTM 37S
gdalinfo -proj4 -nogcp -nomd -norat -noct -nofl 
15JUN11IK0101000po_697515_pan_0000000.ntf

Driver: NITF/National Imagery Transmission Format
Files: 15JUN11IK0101000po_697515_pan_0000000.ntf
Size is 13768, 44528
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",39],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000]]
PROJ.4 string is:
'+proj=utm +zone=37 +south +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m 
+no_defs '
Origin = (545471.000000000000000,9533235.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left  (  545471.000, 9533235.000) ( 39d24'35.06"E,  4d13'22.02"S)
Lower Left  (  545471.000, 9488707.000) ( 39d24'35.85"E,  4d37'32.19"S)
Upper Right (  559239.000, 9533235.000) ( 39d32' 1.67"E,  4d13'21.75"S)
Lower Right (  559239.000, 9488707.000) ( 39d32' 2.71"E,  4d37'31.89"S)
Center      (  552355.000, 9510971.000) ( 39d28'18.81"E,  4d25'26.98"S)
Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Undefined
  Min=0.000 Max=2047.000 
  Minimum=0.000, Maximum=2047.000, Mean=392.957, StdDev=204.156

# use it
grass7 -c 15JUN11IK0101000po_697515_pan_0000000.tif /geo/grassdb/Location_Pan

# where am I in?
g.gisenv

MAPSET=PERMANENT
GISDBASE=/geo/grassdb
LOCATION_NAME=Location_Pan
GUI=text

# where?
g.proj -p

-PROJ_INFO-------------------------------------------------
name       : Lat/Lon
proj       : ll
datum      : wgs84
ellps      : wgs84
no_defs    : defined
-PROJ_UNITS------------------------------------------------
unit       : degree
units      : degrees
meters     : 1.0


This is not correct, is it?


Then, reporting on _all_ multi-spectral bands like

gdalinfo -proj4 -nogcp -nomd -norat -noct -nofl po_697515_red_0000000.ntf | 
grep proj

returns

'+proj=utm +zone=37 +south +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m 
+no_defs '

Using each of the MS bands to create a Location, gives

g.proj -p

XY location (unprojected)

?

Best, Nikos