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

Tue Jul 16 07:46:10 PDT 2013

On Tue, Jul 16, 2013 at 2:10 PM, Nikos Alexandris
<nik at nikosalexandris.net> wrote:
> 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.

It should be possible.

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

Yes, or ogrinfo.
>
>
>> 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...

The projection as recognized by GRASS is indeed different from the one
reported by GDAL.

What says g.proj georef=your_geofile -p?

What reports r.in.gdal when trying to import in a location that
definitively does not match the SRS?

Both tests should be done both with the NITF file and with a
subdataset of the NITF file.

Markus M

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