[GRASS-user] Working with NTF files
Nikos Alexandris
nik at nikosalexandris.net
Sun Jul 14 14:53:18 PDT 2013
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
:-)
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>
[2] <http://trac.osgeo.org/gdal/wiki/JasPer>
[3] <http://trac.osgeo.org/gdal/wiki/JP2KAK>
[4] <http://trac.osgeo.org/gdal/wiki/MrSID>
NITF and GDAL
[5] <http://www.gdal.org/frmt_nitf.html>
[6] <http://www.gdal.org/frmt_nitf_advanced.html>
NITF Related posts (even a bit)
[7] <http://lists.maptools.org/pipermail/fwtools/2008-March/001193.html>,
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>
[10] <http://web.archiveorange.com/archive/v/H2CPnHaFxXHV9S1zxqmE>
[11] <http://lists.osgeo.org/pipermail/grass-user/2000-May/003502.html>, very
old!
[12] <http://lists.osgeo.org/pipermail/grass-user/2006-March/033353.html>, old
[13] <http://lists.osgeo.org/pipermail/grass-user/2006-March/033300.html>, old
[14] <http://osgeo-org.1560.x6.nabble.com/gdal-dev-NITF-with-MultiDataSet-td5050649.html>, "Try building a vrt with the NITF datasets."
[15] <http://www.mail-archive.com/gdal-dev@lists.osgeo.org/msg02376.html>,
refers to "-co".
[16] <http://marc.info/?l=gdal-dev&m=121323053222339&w=2>
[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>, "The images also
have rpc metadata and using the -rpc option in gdalwarp reprojects them quite
well."
More information about the grass-user
mailing list