gdalinfo segfault after gdalwarp of HDF4/MODIS

MODIS gdalpwarp'ed GeoTIFF crashing gdalinfo

Markus,

First, I would suggest you rebuild without using the internal PROJ.4 in
OGDI.  You can configure it with the --with-proj option to use an external
PROJ.4, instead of the very very old copy of the code that lives inside the
OGDI code base.  I think the crash is related to the fact that at least parts
of the OGDI PROJ.4 are getting used. 

When I run gdalinfo against the file it does not crash, but it does report
an effective eccentricity of 1 and a PROJ.4 definition that looks like:

  +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=0 +units=m +no_defs 

(well technically, I get this error with OpenEV, gdalinfo reports no errors
but also doesn't report projected coordinates). 

The gdalinfo does report a very odd coordinate system:

PROJCS["unnamed",
    GEOGCS["Unknown datum based upon the custom spheroid",
        DATUM["unknown",
            SPHEROID["unnamed",6371007.181,1]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

In particular, the SPHEROID has an inverse flattening of 1 which
makes no sense. 

I don't quite understand how the sinusoidal projection got selected
in your gdalwarp.  With my Modis input scenes I am only producing 4
GCPs and they are in WGS84.  So a directly gdalwarp like yours would
result in a lat/long scene, not Sinusoidal.

I would be interested in access to your MODIS scene if that is not a problem.

Frank,

I am using the external PROJ in OGDI:

./configure \
        --with-proj=yes \
        --with-projlib=/usr/local/lib/libproj.so \
        --with-projinc=/usr/local/include \
        --with-zlib

...
checking whether byte ordering is bigendian... no
checking for pj_init in -lproj... yes
checking for projects.h... yes
Using PROJ_LIB=/usr/local/lib/libproj.so for PROJ.4 library.
Using PROJ_INCLUDE=-I/usr/local/include for PROJ.4 includes.
checking for deflateInit_ in -lz... yes
checking for zlib.h... yes
Using internal Expat implementation.
updating cache ./config.cache
...


Oh, I just realize that this is not the case:
[neteler@dandre2 ogdi]$ ldd /usr/local/lib/libogdi.so
        linux-gate.so.1 =>  (0xffffe000)
        libdl.so.2 => /lib/libdl.so.2 (0x40050000)
        libzlib_ogdi31.so => /usr/local/lib/libzlib_ogdi31.so (0x40053000)
        libexpat_ogdi31.so => /usr/local/lib/libexpat_ogdi31.so (0x40060000)
        libm.so.6 => /lib/tls/libm.so.6 (0x4007b000)
        libc.so.6 => /lib/tls/libc.so.6 (0x4009e000)
        /lib/ld-linux.so.2 => /lib/ld-linux.so.2 (0x80000000)
[neteler@dandre2 ogdi]$ ldd /usr/local/lib/libogdi
libogdi31.so  libogdi.so
[neteler@dandre2 ogdi]$ ldd /usr/local/lib/libogdi31.so
        linux-gate.so.1 =>  (0xffffe000)
        libdl.so.2 => /lib/libdl.so.2 (0x40050000)
        libzlib_ogdi31.so => /usr/local/lib/libzlib_ogdi31.so (0x40053000)
        libexpat_ogdi31.so => /usr/local/lib/libexpat_ogdi31.so (0x40060000)
        libm.so.6 => /lib/tls/libm.so.6 (0x4007b000)
        libc.so.6 => /lib/tls/libc.so.6 (0x4009e000)
        /lib/ld-linux.so.2 => /lib/ld-linux.so.2 (0x80000000)

Mysterious. Shouldn't it list PROJ4 as well?

You can grab the scene here:
http://mpa.itc.it/markus/tmp/frank/
(also as tar.gz)

Thanks for your quick response!

 Markus

Markus, 

The external libproj should be shown as a dependency of the libogdi.so.  That
it is not, suggests to me that it is using an internal copy. You can also 
see this clearly in the earlier gdb traceback which shows pj_init() from
libogdi being called. 

I don't know why your configure options aren't taking.  Perhaps you need a 
really clean build?  Perhaps there are other environment variables (like TARGET)
interfering with the configuration selection? 

On further analysis I see that your file is an HDF_GRID, and isn't using
the most recent work with the GRing at all as that is only for swaths. 

However, your grids are picking up georeferencing from normal HDF EOS Grid
projection info.  This was working fine, but it turns out the GeoTIFF writer
was writing spherical earth models as semimajor + 0 inverse flattening.  
The geo_normalize.c code in libgeotiff on read back did not properly recognise
this case and so the earth model was being corrupted. 

I have made two corrections:
 o libgeotiff - geo_normalize.c - an inverse flattening of 0 is interpreted
   as a sphere.  The semi-major axis is set equal to the semi minor. 
 o GDAL - gt_wkt_srs.cpp: Emit user defined spherical earth models as 
   semi-major + semi-minor instead of semi-major + 0 inverse flattening. 

For your purposes, you can rebuild GDAL from CVS and rewarp to GeoTIFF. You
don't need the libgeotiff change. 

I'm closing this bug, but feel free to re-open if the problem persists.

Frank,

ok, now OGDI uses the external PROJ (I had TARGET still enables).
Two aside suggestions for OGDI:

- configure.in
dnl Disable configure caching ... it causes lots of hassles.
dnl define([AC_CACHE_LOAD], )
dnl define([AC_CACHE_SAVE], )

- maybe removal of the interal proj to avoid such problems?

#------------------------

I relaunched gdalwarp, then tried to warp to LatLong:

gdalwarp -t_srs 'EPSG:4326' test2.tif test_ll.tif
ERROR 6: Failed to initialize PROJ.4 with `+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=0 +units=m +no_defs '.
effective eccentricity = 1.
Creating output file that is 2400P x 2400L.
ERROR 6: Failed to initialize PROJ.4 with `+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=0 +units=m +no_defs '.
effective eccentricity = 1.
:0...10...20...30...40...50...60...70...80...90...100 - done.

[neteler@dandre2 egu2005_snowmelt]$ gdalinfo test_ll.tif
Driver: GTiff/GeoTIFF
Size is 2400, 2400
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 = (0.000000,5559752.598333)
Pixel Size = (463.31271653,-463.31271653)
Corner Coordinates:
Upper Left  (       0.000, 5559752.598) (  0d 0'0.00"E,5559752d35'21474836534.00"N)
Lower Left  (       0.000, 4447802.079) (  0d 0'0.00"E,4447802d 4'17179869227.20"N)
Upper Right ( 1111950.520, 5559752.598)
(1111950d31'4294967306.80"E,5559752d35'21474836534.00"N)
Lower Right ( 1111950.520, 4447802.079) (1111950d31'4294967306.80"E,4447802d
4'17179869227.20"N)
Center      (  555975.260, 5003777.338)
(555975d15'35.40"E,5003777d20'17179869202.60"N)
Band 1 Block=2400x3 Type=Byte, ColorInterp=Gray

This should be identical to your former results.

#######################
[Recompilation]

gdalinfo test_ll.tif
Driver: GTiff/GeoTIFF
Size is 3211, 2064
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 = (0.000000,50.189228)
Pixel Size = (0.00484483,-0.00484483)
Corner Coordinates:
Upper Left  (   0.0000000,  50.1892278) (  0d 0'0.00"E, 50d11'21.22"N)
Lower Left  (   0.0000000,  40.1894911) (  0d 0'0.00"E, 40d11'22.17"N)
Upper Right (  15.5567608,  50.1892278) ( 15d33'24.34"E, 50d11'21.22"N)
Lower Right (  15.5567608,  40.1894911) ( 15d33'24.34"E, 40d11'22.17"N)
Center      (   7.7783804,  45.1893595) (  7d46'42.17"E, 45d11'21.69"N)
Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray

Voila'! Works!

Comparing the result to VMAP0 coastline I still observe a shift.
I'll investigate if the "MODIS Reprojection Tool" MRT does some
sort of geometric correction. Probably it's also due to a datum
transformation error in my VMAP0 data, but it looks a bit too much
for that.

Attached a screen shot MODIS vs VMAP0 coastline.

Anyway, it's a great step forward to have geocoding for MODIS
data now!

Best regards

 Markus

MODIS/VMAP0 spatial shift after Grid -> Sinusoidal -> LatLong?

Markus, 

I don't know why you are seeing this significant error.  I didn't really
have the impression that MODIS geolocation information was all that accurate
though. 

Frank,

in general MODIS geocoding is very good. In times of ISIN I observed an
error which was later fixed in the MRT 3.1. The error with a current
MRT version is in the subpixel range (of course it's hard to test
with 1km pixel resolution).

I am not yet well familiar with the snow product. But since it
comes from the same satellite, just processed in a different way,
there should not be such an error.

I'll reopen this bug in case that MRT performs better (and will
post why it performs better if I find out). I'll also compare
to other high-res data as VMAP0 might not be the best test case.

Best regards

 Markus

Frank,

two observations:

* in frmts/hdf4/hdf4imagedataset.cpp I found

  oSRS.SetWellKnownGeogCS( "WGS84" );

Perhaps this could be conditionalized to "SPHERE" for Sinusoidal
projection (if I understand the code)?

* My MODIS data set MYD10A1.A2004244.h18v04.004.2004256140521.hdf
has
GRingPointLatitude:  39.8198,50.0070,49.9990,39.8144
GRingPointLongitude: 0.0001,-0.0087,15.5724,13.0379

while I get after warping to LatLong:

gdalinfo snow_cover_ll2.tif
Driver: GTiff/GeoTIFF
Size is 3211, 2064
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 = (0.000000,50.189228)
Pixel Size = (0.00484483,-0.00484483)
Corner Coordinates:
Upper Left  (   0.0000000,  50.1892278) (  0d 0'0.00"E, 50d11'21.22"N)
Lower Left  (   0.0000000,  40.1894911) (  0d 0'0.00"E, 40d11'22.17"N)
Upper Right (  15.5567608,  50.1892278) ( 15d33'24.34"E, 50d11'21.22"N)
Lower Right (  15.5567608,  40.1894911) ( 15d33'24.34"E, 40d11'22.17"N)
Center      (   7.7783804,  45.1893595) (  7d46'42.17"E, 45d11'21.69"N)
Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray

Maybe this explains the shift? I don't know actually where it comes
from. The reported corners are different from those of the GRing.

Don't beat me :-) I assume that my MODIS Grid data are not what you
intended today. But it seems that it's almost resolved.

Best regards

 Markus

Markus, 

There are several different data paths in the hdf4 code.  Some, such 
as the one deriving the corners from the GRing are explicitely based
on WGS84.  Others, such as this grid get their spheroid definition from
USGS parameters encoded in the file. 

On the other issue, warping to lat/long cannot be expected to give the same
corner coordinates as the GRing corners on the original grid in sinusoidal
coordinates.  

Frank,

for the sake of precision I comment again:

Frank wrote:
> On the other issue, warping to lat/long cannot be expected to give the same
> corner coordinates as the GRing corners on the original grid in sinusoidal
> coordinates.  

This is not clear to me: sinusoidal coordinates are metric while GRing are
LatLong. 

Here the MRT results:

MODIS Reprojection Tool (v3.2a July 2004)
Start Time:  Wed Mar 16 10:40:10 2005
Input image and reprojection info
---------------------------------
input_filename:         
/nfsmnt/levi0/ssi/satellite_data/modis_snow/MYD10A1.A2004244.h18v04.004.2004256140521.hdf
output_filename:         /nfsmnt/levi0/ssi/satellite_data/modis_snow/snow.tif
input_filetype:          HDF-EOS
output_filetype:         GEOTIFF
input_projection_type:   SIN
input_datum:             WGS84
output_projection_type:  GEO
output_datum:            WGS84
resampling_type:         NN
input projection parameters:  6371007.18 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 
output projection parameters: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 

input image corners (lat/lon):
    UL:  50.00 0.00 
    UR:  50.00 15.56 
    LL:  40.00 0.00 
    LR:  40.00 13.05 

input image spatial subset corners (lat/lon):
    UL:  50.00 0.00 
    UR:  50.00 15.56 
    LL:  40.00 0.00 
    LR:  40.00 15.56 

     band          select   type lines smpls pixsiz      min    max   fill
  1) Day_Tile_Snow_Cover   1  UINT8  2400  2400 463.3127      0    255    255
  2) Snow_Spatial_QA    1  UINT8  2400  2400 463.3127      0    255    255
  3) Day_Tile_Snow_Albedo   0  UINT8  2400  2400 463.3127      0    100    255

SINUSOIDAL PROJECTION PARAMETERS:
   Radius of Sphere:     6371007.181000 meters
   Longitude of Center:     0.000000 degrees
   False Easting:      0.000000 meters 
   False Northing:     0.000000 meters 

NNResample : processing band Day_Tile_Snow_Cover
% complete (2065 rows): 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%

NNResample : processing band Snow_Spatial_QA
% complete (2065 rows): 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%

Output image info
-----------------
output image corners (lat/lon):
    UL:  49.999999995507 0.000000000000 
    UR:  49.999999995507 15.561593960000 
    LL:  39.995426045507 0.000000000000 
    LR:  39.995426045507 15.561593960000 

output image corners (X-Y projection units):
    UL:  0.000000000000 49.999999995507 
    UR:  15.561593960000 49.999999995507 
    LL:  0.000000000000 39.995426045507 
    LR:  15.561593960000 39.995426045507 

     band               type lines smpls pixsiz      min    max   fill
  1) Day_Tile_Snow_Cover  UINT8  2065  3212 0.0048      0    255    255
  2) Snow_Spatial_QA   UINT8  2065  3212 0.0048      0    255    255

End Time:  Wed Mar 16 10:40:25 2005
Finished processing!


gdalinfo snow.Day_Tile_Snow_Cover.tif
Driver: GTiff/GeoTIFF
Size is 3212, 2065
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 = (0.002422,49.997578)
Pixel Size = (0.00484483,-0.00484483)
Metadata:
  TIFFTAG_SOFTWARE=MODIS Reprojection Tool  v3.2a July 2004
Corner Coordinates:
Upper Left  (   0.0024224,  49.9975776) (  0d 0'8.72"E, 49d59'51.28"N)
Lower Left  (   0.0024224,  39.9930036) (  0d 0'8.72"E, 39d59'34.81"N)
Upper Right (  15.5640164,  49.9975776) ( 15d33'50.46"E, 49d59'51.28"N)
Lower Right (  15.5640164,  39.9930036) ( 15d33'50.46"E, 39d59'34.81"N)
Center      (   7.7832194,  44.9952906) (  7d46'59.59"E, 44d59'43.05"N)
Band 1 Block=3212x1 Type=Byte, ColorInterp=Gray

I have overlayed this new geotiff to my vector coastline map,
it fit's well.

For comparison:

> 49.999999995507-49.9975776
[1] 0.002422396
> 0.002422396 * 2
[1] 0.004844792

Are we off by half a pixel (corner versus cell center)?

Best regards

 Markus

Frank wrote:
> On the other issue, warping to lat/long cannot be expected to give the same
> corner coordinates as the GRing corners on the original grid in sinusoidal
> coordinates.  

Markus wrote:
>This is not clear to me: sinusoidal coordinates are metric while GRing are
>LatLong. 

My point is that if you take a sinusoidal image, and reproject it to lat/long
you are going to get lat/long corners that are "further out" than the corners
of the input image because the output file is effectively somewhat rotated
and skewed compared to the input (due to the different projection).  

I'm afraid I am not clear on what I am supposed to pick up from the rest
of your (extensive) text.  Is the point that the Modis tool can reproject the
images to Lat/long in a geotiff file, and you get good agreement while 
this is not true if the reprojection is done with gdalwarp?

There could be off-by-half pixel errors in the interpretetion of MODIS grid
georeferencing.  If you want me (or Andrey) to follow up on that possibility
I would suggest you file that as a new bug report.  This one is getting 
unmanagably large to read through. 

If you think there is an error in the gdalwarp handling of the sinusoidal 
projection (or perhaps the HDF driver conversion of USGS Sinusoidal 
parameters to WKT/PROJ.4 format) then I would also suggest that be 
represented as a new bug report.