[GRASS-user] Working with NTF files

GRASS GRASS user
,
1

Hi list & apologies for cross-posting.

In short,

I am trying to handle ntf files (NITF) as easy as possible -- working with
some WorldView 01/02 and QuickBird imagery.

1) The question is: how do you correctly import in GRASS-GIS QuickBird &
WorldView imagery from N(I)TF containers?

2) A second, of less importance, question, is: how important are the warnings
like: "Warning 1: Unable to save auxilary information in
/vsisubfile/3884_471349721,10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf.aux.xml." ?

My working solution (seems to be),

is to gdal-warp the SUBDATASET that contains the raster spectral bands _and_
by using the -rpc switch (mentioned in some past post in GDAL-dev's ML [19]).
Otherwise, the resulting warped image is heavily skewed and not ready for
analysis. E.g.:

--%<---
gdalwarp -rpc NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003_NorthUp.tif
--->%--

This post is, also, sort of a BUMP related to my previous latest message in
the GDAL-dev-ML's thread entitled "Can't read NITF images" :-! I think that
N(I)TF files, no matter whatsoever the status of the related drivers are, as
well as the "feelings" towards favouring it or not, deserves some more
examples, a page in GRASS-Wiki perhaps, etc.

Now, the long story is that,

I 've been through the list (grass-user, gdal-dev), manuals (grass importing
related stuff & gdal's driver's documentation), well known and easy to find
pages in GRASS-Wiki, GDAL tutorials, etc. I haven't found any concrete
examples on working with NITF images and GFOSS. [**See various links in the
end]

I would appreciate any hints (while I am actively searching How To do this
best) to save time in the quest to find an optimal solution which can be
scripted to handle various images from the same source. Unfortunately, there
is no access in GeoTiffs (so far/yet)!

Set aside the internal compression "obstacle" (JPEG2000) which one can
overcome using the OpenJPEG driver (or, necessarily some proprietary (?) bound
solution [0, 1, 2, 3]), NITF files that contain multiple SUBDATASETS may look
like:

# code ###################################################
gdalinfo -nomd -norat -noct -nofl 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf

Driver: NITF/National Imagery Transmission Format
Files: 10APR13WV020600013APR10075059-P1BS-500060446050_05_P003.ntf
Size is 35840, 27648
Coordinate System is `'
GCP Projection =
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"]]
GCP[ 0]: Id=UpperLeft, Info=
          (0.5,0.5) -> (43.8088888888889,-23.3544444444444,0)
GCP[ 1]: Id=UpperRight, Info=
          (35839.5,0.5) -> (43.6219444444444,-23.3344444444444,0)
GCP[ 2]: Id=LowerRight, Info=
          (35839.5,27647.5) -> (43.6225,-23.1958333333333,0)
GCP[ 3]: Id=LowerLeft, Info=
          (0.5,27647.5) -> (43.8083333333333,-23.2152777777778,0)
Subdatasets:
  SUBDATASET_1_NAME=NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
  SUBDATASET_1_DESC=Image 1 of 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
  SUBDATASET_2_NAME=NITF_IM:1:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
  SUBDATASET_2_DESC=Image 2 of 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,27648.0)
Upper Right (35840.0, 0.0)
Lower Right (35840.0,27648.0)
Center (17920.0,13824.0)
Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Gray
  Overviews: 17920x13824, 8960x6912, 4480x3456, 2240x1728, 1120x864
  Overviews: arbitrary
###################################################

The 1st SUBDATASET is reported as:

# code ###################################################
gdalinfo -sd 1 -nomd -norat -noct -nofl 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf

[ -- cut -- some Warnings -- ]

Driver: NITF/National Imagery Transmission Format
Files: 10APR13WV020600013APR10075059-P1BS-500060446050_05_P003.ntf
Size is 35840, 27648
Coordinate System is `'
GCP Projection =
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"]]
GCP[ 0]: Id=UpperLeft, Info=
          (0.5,0.5) -> (43.8088888888889,-23.3544444444444,0)
GCP[ 1]: Id=UpperRight, Info=
          (35839.5,0.5) -> (43.6219444444444,-23.3344444444444,0)
GCP[ 2]: Id=LowerRight, Info=
          (35839.5,27647.5) -> (43.6225,-23.1958333333333,0)
GCP[ 3]: Id=LowerLeft, Info=
          (0.5,27647.5) -> (43.8083333333333,-23.2152777777778,0)
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,27648.0)
Upper Right (35840.0, 0.0)
Lower Right (35840.0,27648.0)
Center (17920.0,13824.0)
Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Gray
  Overviews: 17920x13824, 8960x6912, 4480x3456, 2240x1728, 1120x864
  Overviews: arbitrary

[ -- cut -- some Warnings -- ]
###################################################

The 2nd SUBDATASET is reported as:

# code ###################################################
gdalinfo -sd 2 -nomd -norat -noct -nofl 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf

[ -- cut -- some Warnings -- ]

Driver: NITF/National Imagery Transmission Format
Files: 10APR13WV020600013APR10075059-P1BS-500060446050_05_P003.ntf
Size is 97, 78
Coordinate System is:
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"]]
GeoTransform =
  43.8096587433111, -0.001954571759259283, -1.803751803795437e-06
  -23.35525135093495, 0.0002068865740741814, 0.001823593073593144
Corner Coordinates:
Upper Left ( 43.8096587, -23.3552514) ( 43d48'34.77"E, 23d21'18.90"S)
Lower Left ( 43.8095181, -23.2130111) ( 43d48'34.26"E, 23d12'46.84"S)
Upper Right ( 43.6200653, -23.3351834) ( 43d37'12.24"E, 23d20' 6.66"S)
Lower Right ( 43.6199246, -23.1929431) ( 43d37'11.73"E, 23d11'34.60"S)
Center ( 43.7147917, -23.2740972) ( 43d42'53.25"E, 23d16'26.75"S)
Band 1 Block=97x1 Type=Byte, ColorInterp=Undefined
###################################################

Likewise, reporting on a Multi-Spectral NTF container, the 2nd SUBDATASET is
the Spatial Reference System information.

Unfortunately, though expected, r.in.gdal complains about the missing
projection info for the 1st SUBDATASET.

# code ###################################################
r.in.gdal in=NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf out=m_sd_1
ERROR: Projection of dataset does not appear to match current location.

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

       Import dataset PROJ_INFO is:
       cellhd.proj = 0 (unreferenced/unknown)

       You can use the -o flag to r.in.gdal to override this check and use
       the location definition for the dataset.
       Consider generating a new location from the input dataset using the
       'location' parameter.

[ -- cut -- some Warnings -- ]
###################################################

Adding the "-o" flag, gets to the next error:

# code ###################################################
r.in.gdal in=NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf out=m_sd_1 -o

WARNING: Over-riding projection check
ERROR: Illegal latitude for North

[ -- cut -- some Warnings -- ]
###################################################

Adding the "-l" switch, takes too much time "waiting" in 0% for the import
process. Now, the GCP's are given in the 2nd SUBDATASET(s). Is there any way
to avoid usnig "r.region" and make the long story short, just by using some
gdal-* magic?

Using gdalwarp to go NorthUp:

#################################################### code #
gdalwarp -of NITF -co "ICORDS=G" NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003_NorthUp.ntf
###################################################

allows for imoprting the image via r.in.gdal. However, this image is
obviously and heavily skewed/shifted. Using r.region doesn't make any
sense...

since the image boundaries coincide!

#################################################### code #
gdalinfo -nomd -noct -norat -nofl -nogcp 10APR13WV020600013APR10075059-
M1BS-500060446050_05_P003_NorthUp.tif

Corner Coordinates:
Upper Left ( 43.6218649, -23.1956837) ( 43d37'18.71"E, 23d11'44.46"S)
Lower Left ( 43.6218649, -23.3543246) ( 43d37'18.71"E, 23d21'15.57"S)
Upper Right ( 43.8088360, -23.1956837) ( 43d48'31.81"E, 23d11'44.46"S)
Lower Right ( 43.8088360, -23.3543246) ( 43d48'31.81"E, 23d21'15.57"S)
###################################################

and minor differences when gdal-warping in to a NITF file

#################################################### code #
Corner Coordinates:
Upper Left ( 43.6219339, -23.1955450) ( 43d37'18.96"E, 23d11'43.96"S)
Lower Left ( 43.6219339, -23.3544550) ( 43d37'18.96"E, 23d21'16.04"S)
Upper Right ( 43.8088994, -23.1955450) ( 43d48'32.04"E, 23d11'43.96"S)
Lower Right ( 43.8088994, -23.3544550) ( 43d48'32.04"E, 23d21'16.04"S)
###################################################

After importing, one band reports

#################################################### code #
r.info -g M1BS_from_TIF.1

north=-23.1956836847222
south=-23.3543246272222
east=43.8088360363889
west=43.6218648544444
###################################################

The corner coordinates are practically identical. Hence, it doesn't look like
boundary (re-)definition/correction is required for these rasters.

So, how do you go about your NITF images?

The only "nice-working" attempt so far, is when I extacted a specific fragment
out of the Multi-Spectral container, for example, using the "-te" switch while
warping, i.e. (using values, of course, in place of West, South, North and
East):

gdalwarp -s_srs 'epsg:4326' -te West South North East
03APR13WV010500013APR03073141-P1BS-500060446050_03_P003.ntf test.tif

These small fragments coincide perfectly with some vectorised area of intere
st.

And, finally, and fortunately, adding the "-rpc" magic, seems to do the trick.
Importing an rpc-warped image in GRASS, looks, after some cross-checks, nice
:slight_smile:

Any other solutions?

Thanks, Nikos

---
[0]
<http://www.jpeg.org/faq.phtml?action=show_answer&question_id=q3f042a68b1081&gt;
[1] <http://trac.osgeo.org/gdal/wiki/ECW&gt;
[2] <http://trac.osgeo.org/gdal/wiki/JasPer&gt;
[3] <http://trac.osgeo.org/gdal/wiki/JP2KAK&gt;
[4] <http://trac.osgeo.org/gdal/wiki/MrSID&gt;

NITF and GDAL

[5] <http://www.gdal.org/frmt_nitf.html&gt;
[6] <http://www.gdal.org/frmt_nitf_advanced.html&gt;

NITF Related posts (even a bit)

[7] <http://lists.maptools.org/pipermail/fwtools/2008-March/001193.html&gt;,
2008, refers also to license stuff.
[8] <http://osgeo-org.1560.x6.nabble.com/gdal-dev-Can-t-read-NITF-image-td5064590.html&gt;
[9] <http://lists.osgeo.org/pipermail/gdal-dev/2012-April/032431.html&gt;
[10] <http://web.archiveorange.com/archive/v/H2CPnHaFxXHV9S1zxqmE&gt;
[11] <http://lists.osgeo.org/pipermail/grass-user/2000-May/003502.html&gt;, very
old!
[12] <http://lists.osgeo.org/pipermail/grass-user/2006-March/033353.html&gt;, old
[13] <http://lists.osgeo.org/pipermail/grass-user/2006-March/033300.html&gt;, old
[14] <http://osgeo-org.1560.x6.nabble.com/gdal-dev-NITF-with-MultiDataSet-td5050649.html&gt;, "Try building a vrt with the NITF datasets."
[15] <http://www.mail-archive.com/gdal-dev@lists.osgeo.org/msg02376.html&gt;,
refers to "-co".
[16] <http://marc.info/?l=gdal-dev&m=121323053222339&w=2&gt;
[17] <http://osdir.com/ml/gdal-development-gis-osgeo/2009-01/msg00044.html&gt;,
about compression/compressed output.

Advanced coding discussions

[18] <http://osgeo-org.1560.x6.nabble.com/NITF-image-metadata-domain-field-in-GDAL-SUBDATASETS-td3762940.html&gt;
[19] <http://osgeo-org.1560.x6.nabble.com/gdal-dev-How-to-use-RPC-Metadata-from-NITF-Files-in-a-MapServer-TileIndex-td4981951.html&gt;, "The images also
have rpc metadata and using the -rpc option in gdalwarp reprojects them quite
well."

2

On Sun, Jul 14, 2013 at 2:53 PM, Nikos Alexandris
<nik@nikosalexandris.net>wrote:

Hi list & apologies for cross-posting.

In short,

I am trying to handle ntf files (NITF) as easy as possible -- working with
some WorldView 01/02 and QuickBird imagery.

1) The question is: how do you correctly import in GRASS-GIS QuickBird &
WorldView imagery from N(I)TF containers?

Nikos,

Your working solution with gdalwarp seems reasonable. Is there a problem
with it? Perhaps you were hoping to do the RPC based warp within GRASS?

2) A second, of less importance, question, is: how important are the
warnings
like: "Warning 1: Unable to save auxilary information in
/vsisubfile/3884_471349721,10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf.aux.xml." ?

This is not important. I presume you are getting this with JPEG2000
compressed
image streams in NITF? If you file a bug I might be able to clear this up
but there
are layers of complexity that make it challening.

My working solution (seems to be),

is to gdal-warp the SUBDATASET that contains the raster spectral bands
_and_
by using the -rpc switch (mentioned in some past post in GDAL-dev's ML
[19]).
Otherwise, the resulting warped image is heavily skewed and not ready for
analysis. E.g.:

--%<---
gdalwarp -rpc NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003_NorthUp.tif
--->%--

This post is, also, sort of a BUMP related to my previous latest message in
the GDAL-dev-ML's thread entitled "Can't read NITF images" :-! I think
that
N(I)TF files, no matter whatsoever the status of the related drivers are,
as
well as the "feelings" towards favouring it or not, deserves some more
examples, a page in GRASS-Wiki perhaps, etc.

I'm afraid I'm one of those who neglected that thread. I thought I saw Even
answer is so I was happy to keep out of it.

I can't speak to how the files should be described in GRASS. Generally
I'd hope that the value of using GDAL is that GRASS wouldn't need to
talk too much about specific formats.

On the GDAL side we often have special info in trac wiki pages under the
BuildHints page, but in this case the issues are more usage rather than
building so there is no obvious place for user contributions. The format
pages for NITF do have quite a bit of info. It is (I think) the only driver
with an "advanced" page. However, there are many kinds of NITF
files and thus usage patterns that are an issue so it isn't clear how to
handle that in the user docs.

I've skimmed the rest of your post, but I don't have any particular
advice or anything to add. I do think it is reasonable to correct
such images into a nicer form with gdalwarp before running r.in.gdal.

Best regards,
Frank

Now, the long story is that,

I 've been through the list (grass-user, gdal-dev), manuals (grass
importing
related stuff & gdal's driver's documentation), well known and easy to find
pages in GRASS-Wiki, GDAL tutorials, etc. I haven't found any concrete
examples on working with NITF images and GFOSS. [**See various links in the
end]

I would appreciate any hints (while I am actively searching How To do this
best) to save time in the quest to find an optimal solution which can be
scripted to handle various images from the same source. Unfortunately,
there
is no access in GeoTiffs (so far/yet)!

Set aside the internal compression "obstacle" (JPEG2000) which one can
overcome using the OpenJPEG driver (or, necessarily some proprietary (?)
bound
solution [0, 1, 2, 3]), NITF files that contain multiple SUBDATASETS may
look
like:

# code ###################################################
gdalinfo -nomd -norat -noct -nofl 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf

Driver: NITF/National Imagery Transmission Format
Files: 10APR13WV020600013APR10075059-P1BS-500060446050_05_P003.ntf
Size is 35840, 27648
Coordinate System is `'
GCP Projection =
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"]]
GCP[ 0]: Id=UpperLeft, Info=
          (0.5,0.5) -> (43.8088888888889,-23.3544444444444,0)
GCP[ 1]: Id=UpperRight, Info=
          (35839.5,0.5) -> (43.6219444444444,-23.3344444444444,0)
GCP[ 2]: Id=LowerRight, Info=
          (35839.5,27647.5) -> (43.6225,-23.1958333333333,0)
GCP[ 3]: Id=LowerLeft, Info=
          (0.5,27647.5) -> (43.8083333333333,-23.2152777777778,0)
Subdatasets:
  SUBDATASET_1_NAME=NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
  SUBDATASET_1_DESC=Image 1 of 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
  SUBDATASET_2_NAME=NITF_IM:1:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
  SUBDATASET_2_DESC=Image 2 of 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,27648.0)
Upper Right (35840.0, 0.0)
Lower Right (35840.0,27648.0)
Center (17920.0,13824.0)
Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Gray
  Overviews: 17920x13824, 8960x6912, 4480x3456, 2240x1728, 1120x864
  Overviews: arbitrary
###################################################

The 1st SUBDATASET is reported as:

# code ###################################################
gdalinfo -sd 1 -nomd -norat -noct -nofl 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf

[ -- cut -- some Warnings -- ]

Driver: NITF/National Imagery Transmission Format
Files: 10APR13WV020600013APR10075059-P1BS-500060446050_05_P003.ntf
Size is 35840, 27648
Coordinate System is `'
GCP Projection =
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"]]
GCP[ 0]: Id=UpperLeft, Info=
          (0.5,0.5) -> (43.8088888888889,-23.3544444444444,0)
GCP[ 1]: Id=UpperRight, Info=
          (35839.5,0.5) -> (43.6219444444444,-23.3344444444444,0)
GCP[ 2]: Id=LowerRight, Info=
          (35839.5,27647.5) -> (43.6225,-23.1958333333333,0)
GCP[ 3]: Id=LowerLeft, Info=
          (0.5,27647.5) -> (43.8083333333333,-23.2152777777778,0)
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,27648.0)
Upper Right (35840.0, 0.0)
Lower Right (35840.0,27648.0)
Center (17920.0,13824.0)
Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Gray
  Overviews: 17920x13824, 8960x6912, 4480x3456, 2240x1728, 1120x864
  Overviews: arbitrary

[ -- cut -- some Warnings -- ]
###################################################

The 2nd SUBDATASET is reported as:

# code ###################################################
gdalinfo -sd 2 -nomd -norat -noct -nofl 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf

[ -- cut -- some Warnings -- ]

Driver: NITF/National Imagery Transmission Format
Files: 10APR13WV020600013APR10075059-P1BS-500060446050_05_P003.ntf
Size is 97, 78
Coordinate System is:
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"]]
GeoTransform =
  43.8096587433111, -0.001954571759259283, -1.803751803795437e-06
  -23.35525135093495, 0.0002068865740741814, 0.001823593073593144
Corner Coordinates:
Upper Left ( 43.8096587, -23.3552514) ( 43d48'34.77"E, 23d21'18.90"S)
Lower Left ( 43.8095181, -23.2130111) ( 43d48'34.26"E, 23d12'46.84"S)
Upper Right ( 43.6200653, -23.3351834) ( 43d37'12.24"E, 23d20' 6.66"S)
Lower Right ( 43.6199246, -23.1929431) ( 43d37'11.73"E, 23d11'34.60"S)
Center ( 43.7147917, -23.2740972) ( 43d42'53.25"E, 23d16'26.75"S)
Band 1 Block=97x1 Type=Byte, ColorInterp=Undefined
###################################################

Likewise, reporting on a Multi-Spectral NTF container, the 2nd SUBDATASET
is
the Spatial Reference System information.

Unfortunately, though expected, r.in.gdal complains about the missing
projection info for the 1st SUBDATASET.

# code ###################################################
r.in.gdal in=NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf out=m_sd_1
ERROR: Projection of dataset does not appear to match current location.

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

       Import dataset PROJ_INFO is:
       cellhd.proj = 0 (unreferenced/unknown)

       You can use the -o flag to r.in.gdal to override this check and use
       the location definition for the dataset.
       Consider generating a new location from the input dataset using the
       'location' parameter.

[ -- cut -- some Warnings -- ]
###################################################

Adding the "-o" flag, gets to the next error:

# code ###################################################
r.in.gdal in=NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf out=m_sd_1 -o

WARNING: Over-riding projection check
ERROR: Illegal latitude for North

[ -- cut -- some Warnings -- ]
###################################################

Adding the "-l" switch, takes too much time "waiting" in 0% for the import
process. Now, the GCP's are given in the 2nd SUBDATASET(s). Is there any
way
to avoid usnig "r.region" and make the long story short, just by using some
gdal-* magic?

Using gdalwarp to go NorthUp:

#################################################### code #
gdalwarp -of NITF -co "ICORDS=G" NITF_IM:0:10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003.ntf 10APR13WV020600013APR10075059-
P1BS-500060446050_05_P003_NorthUp.ntf
###################################################

allows for imoprting the image via r.in.gdal. However, this image is
obviously and heavily skewed/shifted. Using r.region doesn't make any
sense...

since the image boundaries coincide!

#################################################### code #
gdalinfo -nomd -noct -norat -nofl -nogcp 10APR13WV020600013APR10075059-
M1BS-500060446050_05_P003_NorthUp.tif

Corner Coordinates:
Upper Left ( 43.6218649, -23.1956837) ( 43d37'18.71"E, 23d11'44.46"S)
Lower Left ( 43.6218649, -23.3543246) ( 43d37'18.71"E, 23d21'15.57"S)
Upper Right ( 43.8088360, -23.1956837) ( 43d48'31.81"E, 23d11'44.46"S)
Lower Right ( 43.8088360, -23.3543246) ( 43d48'31.81"E, 23d21'15.57"S)
###################################################

and minor differences when gdal-warping in to a NITF file

#################################################### code #
Corner Coordinates:
Upper Left ( 43.6219339, -23.1955450) ( 43d37'18.96"E, 23d11'43.96"S)
Lower Left ( 43.6219339, -23.3544550) ( 43d37'18.96"E, 23d21'16.04"S)
Upper Right ( 43.8088994, -23.1955450) ( 43d48'32.04"E, 23d11'43.96"S)
Lower Right ( 43.8088994, -23.3544550) ( 43d48'32.04"E, 23d21'16.04"S)
###################################################

After importing, one band reports

#################################################### code #
r.info -g M1BS_from_TIF.1

north=-23.1956836847222
south=-23.3543246272222
east=43.8088360363889
west=43.6218648544444
###################################################

The corner coordinates are practically identical. Hence, it doesn't look
like
boundary (re-)definition/correction is required for these rasters.

So, how do you go about your NITF images?

The only "nice-working" attempt so far, is when I extacted a specific
fragment
out of the Multi-Spectral container, for example, using the "-te" switch
while
warping, i.e. (using values, of course, in place of West, South, North and
East):

gdalwarp -s_srs 'epsg:4326' -te West South North East
03APR13WV010500013APR03073141-P1BS-500060446050_03_P003.ntf test.tif

These small fragments coincide perfectly with some vectorised area of
intere
st.

And, finally, and fortunately, adding the "-rpc" magic, seems to do the
trick.
Importing an rpc-warped image in GRASS, looks, after some cross-checks,
nice
:slight_smile:

Any other solutions?

Thanks, Nikos

---
[0]
<
http://www.jpeg.org/faq.phtml?action=show_answer&question_id=q3f042a68b1081
>
[1] <http://trac.osgeo.org/gdal/wiki/ECW&gt;
[2] <http://trac.osgeo.org/gdal/wiki/JasPer&gt;
[3] <http://trac.osgeo.org/gdal/wiki/JP2KAK&gt;
[4] <http://trac.osgeo.org/gdal/wiki/MrSID&gt;

NITF and GDAL

[5] <http://www.gdal.org/frmt_nitf.html&gt;
[6] <http://www.gdal.org/frmt_nitf_advanced.html&gt;

NITF Related posts (even a bit)

[7] <http://lists.maptools.org/pipermail/fwtools/2008-March/001193.html&gt;,
2008, refers also to license stuff.
[8] <
http://osgeo-org.1560.x6.nabble.com/gdal-dev-Can-t-read-NITF-image-td5064590.html
>
[9] <http://lists.osgeo.org/pipermail/gdal-dev/2012-April/032431.html&gt;
[10] <http://web.archiveorange.com/archive/v/H2CPnHaFxXHV9S1zxqmE&gt;
[11] <http://lists.osgeo.org/pipermail/grass-user/2000-May/003502.html&gt;,
very
old!
[12] <http://lists.osgeo.org/pipermail/grass-user/2006-March/033353.html&gt;,
old
[13] <http://lists.osgeo.org/pipermail/grass-user/2006-March/033300.html&gt;,
old
[14] <
http://osgeo-org.1560.x6.nabble.com/gdal-dev-NITF-with-MultiDataSet-td5050649.html&gt;,
"Try building a vrt with the NITF datasets."
[15] <http://www.mail-archive.com/gdal-dev@lists.osgeo.org/msg02376.html&gt;,
refers to "-co".
[16] <http://marc.info/?l=gdal-dev&m=121323053222339&w=2&gt;
[17] <http://osdir.com/ml/gdal-development-gis-osgeo/2009-01/msg00044.html
>,
about compression/compressed output.

Advanced coding discussions

[18] <
http://osgeo-org.1560.x6.nabble.com/NITF-image-metadata-domain-field-in-GDAL-SUBDATASETS-td3762940.html
>
[19] <
http://osgeo-org.1560.x6.nabble.com/gdal-dev-How-to-use-RPC-Metadata-from-NITF-Files-in-a-MapServer-TileIndex-td4981951.html&gt;,
"The images also
have rpc metadata and using the -rpc option in gdalwarp reprojects them
quite
well."
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam,
warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Software Developer

3

Nikos Alexandris:

> I am trying to handle ntf files (NITF) as easy as possible -- working with
> some WorldView 01/02 and QuickBird imagery.

> 1) The question is: how do you correctly import in GRASS-GIS QuickBird &
> WorldView imagery from N(I)TF containers?

Frank Warmerdam wrote:

Your working solution with gdalwarp seems reasonable. Is there a problem
with it?

No (except that I left this option as the "final attempt" of many tests --
working with some >1GB images here, so it took me hours... :D). Your
confirmation, though, is highly useful!

Perhaps you were hoping to do the RPC based warp within GRASS?

Yes, I was hopping to say something like:

--%<---
r.in.gdal in=NITF.ntf out=NITF sds=1 band=1,2,3,4 transform=rpc[,tps]
[resampling=n,b,c,cs,l,a,m]
--->%--

:smiley: -- I know, I am asking for too much. I have read the following post
already: <http://fwarmerdam.blogspot.gr/2009/12/death-by-complexity.html&gt;\.

Also, there is already a (kind of related) grass-ticket:
<http://trac.osgeo.org/grass/ticket/798&gt;\.

> 2) A second, of less importance, question, is: how important are the
> warnings like: "Warning 1: Unable to save auxilary information in
> /vsisubfile/3884_471349721,10APR13WV020600013APR10075059-
> P1BS-500060446050_05_P003.ntf.aux.xml." ?

This is not important. I presume you are getting this with JPEG2000
compressed image streams in NITF?

Correct

If you file a bug I might be able to clear this up but there are layers of
complexity that make it challening.

I will (take the time later to) file a ticket. I was hopping, anyhow, in case
such warnings are not a real big deal, to have a way to supress them. That
would ease off, me thinks at least, scripting. Some option like the "-q" one
-- though, in this case it doesn't silence these Warnings.

> My working solution (seems to be), is to gdal-warp the SUBDATASET that
> contains the raster spectral bands _and_ by using the -rpc switch
> (mentioned in some past post in GDAL-dev's ML [19]). Otherwise, the
> resulting warped image is heavily skewed and not ready for analysis.
> E.g.:
> --%<---
> gdalwarp -rpc NITF_IM:0:10APR13WV020600013APR10075059-
> P1BS-500060446050_05_P003.ntf 10APR13WV020600013APR10075059-
> P1BS-500060446050_05_P003_NorthUp.tif
> --->%--

> This post is, also, sort of a BUMP related to my previous latest message
> in the GDAL-dev-ML's thread entitled "Can't read NITF images" :-! I
> think that N(I)TF files, no matter whatsoever the status of the related
> drivers are, as well as the "feelings" towards favouring it or not,
> deserves some more examples, a page in GRASS-Wiki perhaps, etc.

I'm afraid I'm one of those who neglected that thread. I thought I saw Even
answer is so I was happy to keep out of it.

I can't speak to how the files should be described in GRASS. Generally
I'd hope that the value of using GDAL is that GRASS wouldn't need to
talk too much about specific formats.

It certainly is so (in general).

On the GDAL side we often have special info in trac wiki pages under the
BuildHints page, but in this case the issues are more usage rather than
building so there is no obvious place for user contributions. The format
pages for NITF do have quite a bit of info. It is (I think) the only driver
with an "advanced" page. However, there are many kinds of NITF
files and thus usage patterns that are an issue so it isn't clear how to
handle that in the user docs.

I am (slowly) getting to see that there is a "bigger" picture.

I've skimmed the rest of your post, but I don't have any particular
advice or anything to add. I do think it is reasonable to correct
such images into a nicer form with gdalwarp before running r.in.gdal.

Your reply as a whole is very useful for me.

Thank you Frank, Nikos
--

[rest deleted]

4

On Mon, Jul 15, 2013 at 1:52 AM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

Nikos Alexandris:

Perhaps you were hoping to do the RPC based warp within GRASS?

Yes, I was hopping to say something like:

--%<---
r.in.gdal in=NITF.ntf out=NITF sds=1 band=1,2,3,4 transform=rpc[,tps]
[resampling=n,b,c,cs,l,a,m]
--->%--

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.

Markus M

5

Frank Warmerdam:

>> Perhaps you were hoping to do the RPC based warp within GRASS?

Nikos Alexandris:

> Yes, I was hopping to say something like:

( Yeah, right... hopping :smiley: )

> --%<---
> r.in.gdal in=NITF.ntf out=NITF sds=1 band=1,2,3,4 transform=rpc[,tps]
> [resampling=n,b,c,cs,l,a,m]
> --->%--

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.

Ok -- On/Off: would it make sense to be able to say

grass -c Some_NITF_Container /grassdb/Location_Based_On_NITF

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

?

Nikos

6

Nikos A:

> No (except that I left this option as the "final attempt" of many tests --
> working with some >1GB images here, so it took me hours... :D). Your
> confirmation, though, is highly useful!

Eli Adam wrote:

Sometimes to avoid materializing large intermediate datasets on disk, a VRT
can be used. I'm not sure if this would be useful in your case.

Right, in theory it is extremely useful. Dunno if it'll work. In my To Do
list... :slight_smile:

Frank W:

> On the GDAL side we often have special info in trac wiki pages under the
> BuildHints page, but in this case the issues are more usage rather than
> building so there is no obvious place for user contributions. The format
> pages for NITF do have quite a bit of info. It is (I think) the only
> driver with an "advanced" page. However, there are many kinds of NITF
> files and thus usage patterns that are an issue so it isn't clear how to
> handle that in the user docs.

I and others have put GDAL usage notes here,
http://trac.osgeo.org/gdal/wiki/CodeSnippets

Thank you, Nikos

7

On Tue, Jul 16, 2013 at 12:34 PM, Nikos Alexandris
<nik@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.

Ok -- On/Off: would it make sense to be able to say

grass -c Some_NITF_Container /grassdb/Location_Based_On_NITF

You probably mean

grass -c Some_GDAL/OGR_recognized_Container
/grassdb/Location_Based_On_GDAL/OGR_recognized_dataset

which is already 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,
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.

Markus M

8

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

9

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

10

Markus Metz:

[almost all previous discussion erased]

Following a test using an IKONOS-NTF set of files – they do not contain SUBDATASETS. Actually, the same files as in my last/previous post. Will re-test also for a NITF container that includes SUBDATASETs.

What says g.proj georef=your_geofile -p?

That’s a very nice parameter which I never thought of to combine with “-p” :slight_smile:

,—8<–

g.proj -p georef=15JUN11IK0101000po_697515_pan_0000000.ntf

Trying to open with OGR…

…succeeded.

WARNING: Read of file 15JUN11IK0101000po_697515_pan_0000000.ntf was

successful, but it did not contain projection information. 'XY

(unprojected)’ will be used

XY location (unprojected)

`–>%—

This is however 1) strange – it misses a message like “Trying to open with GDAL…” and 2) (successively) not true because gdalinfo reports:

,—8<–

gdalinfo -proj4 -nomd -noct -nogcp 15JUN11IK0101000po_697515_pan_0000000.ntf | grep proj

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

`–>%—

Exactly the same happens (as above) for the multi-spectral bands. A bug? Using this file in

“grass7 -c 15JUN11IK0101000po_697515_pan_0000000.tif /grassdb/test_Location”

creates the following Location:

,—%<–

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

`–>%—

What reports r.in.gdal when trying to import in a location that

definitively does not match the SRS?

Now, r.in.gdal seems to correctly recognise the SRS!

,—8<–

r.in.gdal in=15JUN11IK0101000po_697515_pan_0000000.ntf out=import_test_15JUN11IK0101000po_697515_pan_0000000.ntf

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

Location PROJ_INFO is:

name: Lat/Lon

proj: ll

datum: wgs84

ellps: wgs84

no_defs: defined

Dataset PROJ_INFO is:

name: Universal Transverse Mercator

proj: utm

datum: wgs84

ellps: wgs84

zone: 37

south: defined

towgs84: 0,0,0,0,0,0,0

no_defs: defined

ERROR:

You can use the -o flag to r.in.gdal to override this check and use

the location definition for the dataset.

Consider generating a new location from the input dataset using the

‘location’ parameter.

`–>%—

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

subdataset of the NITF file.

Not applicable in this case, no SUBDATASETs in these file(s).

Nikos

11

On Sat, Jul 20, 2013 at 9:20 PM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

1)

What says g.proj georef=your_geofile -p?

g.proj -p georef=15JUN11IK0101000po_697515_pan_0000000.ntf

Trying to open with OGR...
...succeeded.
WARNING: Read of file 15JUN11IK0101000po_697515_pan_0000000.ntf was
successful, but it did not contain projection information. 'XY
(unprojected)' will be used
XY location (unprojected)

This is however 1) strange -- it misses a message like "Trying to open with
GDAL..." and 2) (successively) not true because gdalinfo reports:

g.proj first tries to open the geofile with OGR, if that fails, it
tries to open it with GDAL. Apparently OGR found a driver with which
it could successfully open the geofile. This might be a bug in OGR
(Frank?).

Using this file in

"grass7 -c 15JUN11IK0101000po_697515_pan_0000000.tif /grassdb/test_Location"

creates the following Location:

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

grass7 -c uses g.proj to create a location, thus the resulting
projection info of grass7 -c and g.proj georef= should be identical.

2)

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

Now, r.in.gdal seems to correctly recognise the SRS!

Good.

So it is either a bug in g.proj (even though the call to OGROpen()
succeeds) or in OGR.

Markus M

12

Nikos Alexandris:

> > I am trying to handle ntf files (NITF) as easy as possible -- working
> > with some WorldView 01/02 and QuickBird imagery.

[..]

> > 2) A second, of less importance, question, is: how important are the
> > warnings like: "Warning 1: Unable to save auxilary information in
> > /vsisubfile/3884_471349721,10APR13WV020600013APR10075059-
> > P1BS-500060446050_05_P003.ntf.aux.xml." ?

Frank W:

> This is not important. I presume you are getting this with JPEG2000
> compressed image streams in NITF?

Nikos A:

Correct

> If you file a bug I might be able to clear this up but there are layers of
> complexity that make it challening.

I will (take the time later to) file a ticket.

Done as http://trac.osgeo.org/gdal/ticket/5173

Nikos

13

Should I file for this a ticket in GRASS' trac first?

Nikos

--%<---

1)

Markus Metz wrote:

>> What says g.proj georef=your_geofile -p?

Nikos Alexandris:

> g.proj -p georef=15JUN11IK0101000po_697515_pan_0000000.ntf

> Trying to open with OGR...
> ...succeeded.
> WARNING: Read of file 15JUN11IK0101000po_697515_pan_0000000.ntf was
> successful, but it did not contain projection information. 'XY
> (unprojected)' will be used XY location (unprojected)

> This is however 1) strange -- it misses a message like "Trying to open
> with GDAL..." and 2) (successively) not true because gdalinfo reports:

Markus M:

g.proj first tries to open the geofile with OGR, if that fails, it
tries to open it with GDAL. Apparently OGR found a driver with which
it could successfully open the geofile. This might be a bug in OGR
(Frank?).

Nikos A:

> Using this file in
> "grass7 -c 15JUN11IK0101000po_697515_pan_0000000.tif
> /grassdb/test_Location"

> creates the following Location:
> 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

Markus M:

grass7 -c uses g.proj to create a location, thus the resulting
projection info of grass7 -c and g.proj georef= should be identical.

2)

Markus M:

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

Nikos A:

> Now, r.in.gdal seems to correctly recognise the SRS!

Markus M:

Good. So it is either a bug in g.proj (even though the call to OGROpen()
succeeds) or in OGR.

14

Yet another NITF related "problem". An ntf file contains (no SUBDATASETS)

gdalinfo -nogcp -nomd -mm M1BS.ntf

,--%<---
Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Blue
    Computed Min/Max=183.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
Band 2 Block=1024x1024 Type=UInt16, ColorInterp=Green
    Computed Min/Max=252.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
Band 3 Block=1024x1024 Type=UInt16, ColorInterp=Red
    Computed Min/Max=86.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
Band 4 Block=1024x1024 Type=UInt16, ColorInterp=Alpha
    Computed Min/Max=50.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
`--->%--

Attempting to warp the data with the -rpc switch and convert to GeoTIFF,
begins with

"Using band 4 of source image as alpha."

This results in an "empty" NIR band in the GeoTIFF file.

,--%<---
Band 1 Block=6728x1 Type=UInt16, ColorInterp=Gray
    Computed Min/Max=0.000,2047.000
Band 2 Block=6728x1 Type=UInt16, ColorInterp=Undefined
    Computed Min/Max=0.000,2047.000
Band 3 Block=6728x1 Type=UInt16, ColorInterp=Undefined
    Computed Min/Max=0.000,2047.000
Band 4 Block=6728x1 Type=UInt16, ColorInterp=Undefined
    Computed Min/Max=0.000,0.000
`--->%--

No success using '-co "ALPHA=NO"'. Same "problem" when using gdal_translate
on the whole ntf file. However, instructing

gdal_translate -b 4 M1BS.ntf b4.tif

gives

gdalinfo -nogcp -nomd -mm b4.tif

,--%<---
Band 1 Block=6876x1 Type=UInt16, ColorInterp=Gray
    Computed Min/Max=50.000,2047.000
`--->%--

But, this is not "fine". How is this to be done correctly, RPC-based warping
and conversion in one go?

Thanks, Nikos

15

On Wednesday 31 of July 2013 11:55:27 Nikos Alexandris wrote:

Yet another NITF related "problem". An ntf file contains (no SUBDATASETS)

gdalinfo -nogcp -nomd -mm M1BS.ntf

,--%<---
Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Blue
    Computed Min/Max=183.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
Band 2 Block=1024x1024 Type=UInt16, ColorInterp=Green
    Computed Min/Max=252.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
Band 3 Block=1024x1024 Type=UInt16, ColorInterp=Red
    Computed Min/Max=86.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
Band 4 Block=1024x1024 Type=UInt16, ColorInterp=Alpha
    Computed Min/Max=50.000,2047.000
  Overviews: 3438x3448, 1719x1724, 859x862, 429x431, 214x215
  Overviews: arbitrary
`--->%--

Attempting to warp the data with the -rpc switch and convert to GeoTIFF,
begins with

"Using band 4 of source image as alpha."

This results in an "empty" NIR band in the GeoTIFF file.

,--%<---
Band 1 Block=6728x1 Type=UInt16, ColorInterp=Gray
    Computed Min/Max=0.000,2047.000
Band 2 Block=6728x1 Type=UInt16, ColorInterp=Undefined
    Computed Min/Max=0.000,2047.000
Band 3 Block=6728x1 Type=UInt16, ColorInterp=Undefined
    Computed Min/Max=0.000,2047.000
Band 4 Block=6728x1 Type=UInt16, ColorInterp=Undefined
    Computed Min/Max=0.000,0.000
`--->%--

No success using '-co "ALPHA=NO"'. Same "problem" when using gdal_translate
on the whole ntf file. However, instructing

gdal_translate -b 4 M1BS.ntf b4.tif

gives

gdalinfo -nogcp -nomd -mm b4.tif

,--%<---
Band 1 Block=6876x1 Type=UInt16, ColorInterp=Gray
    Computed Min/Max=50.000,2047.000
`--->%--

But, this is not "fine". How is this to be done correctly, RPC-based warping
and conversion in one go?

Related thread/solution:

http://lists.osgeo.org/pipermail/gdal-dev/2012-December/034984.html

Translating to a VRT, band by band -- is it the only solution?

Nikos