Opened 19 years ago
Last modified 19 years ago
#805 closed defect (fixed)
hdf driver/gdalwarp: MODIS Sinusoidal grid seems to be shifted by 0.5 pixel
Reported by: | Markus Neteler | Owned by: | warmerdam |
---|---|---|---|
Priority: | high | Milestone: | |
Component: | GDAL_Raster | Version: | unspecified |
Severity: | normal | Keywords: | |
Cc: |
Description
Frank, comparing the 'gdalwarp' results to those of the MODIS Reprojection tool (MRT), it appears that 'gdalwarp' results seem to be shifted by 0.5 pixel. Procedures: (1) MODIS Sinusoidal Grid -> MRT -> LatLong/WGS84: matches other coastline data (2) MODIS Sinusoidal Grid -> gdalwarp -> Sinusoidal -> gdalwarp -> LatLong/WGS84: shifted to other reference data In details (test file: MYD10A1.A2004244.h18v04.004.2004256140521.hdf) GRingPointLatitude: 39.8198,50.0070,49.9990,39.8144 GRingPointLongitude: 0.0001,-0.0087,15.5724,13.0379 #################################################### Case (1) MODIS Reprojection tool (MRT): 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 The result seen in GDAL: gdalinfo snow.Day_Tile_Snow_Cover.tif Size is 3212, 2065 ... 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 -> This looks ok. So the MRT produces correct results. #################################################### Case (2) rewarping with 'gdalwarp': - using gdalwarp on the HDF file, then again gdalwarp to produce LatLong: gdalinfo snow_cover_ll2.tif Driver: GTiff/GeoTIFF Size is 3211, 2064 ... 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 -> Size: we have lost 1 row, 1 column -> result of gdalwarp is shifted #################################################### Verification of the shift, looking at North value of upper left corner: 49.999999995507-49.9975776 = 0.002422396 0.002422396 * 2 = 0.004844792 -> suggests that gdalwarp operates on MODIS Sinusoida Grid with 0.5 pixel shift. Probably there is a mistake in confusing cell corner coordinates with cell center coordinates. Summary: The results of case (2),gdalwarp, should be identical to case (1), MRT. The MRT results are correct. Best regards Markus
Attachments (7)
Change History (30)
comment:2 by , 19 years ago
The MODIS scene here which I used for the experiments is here: http://mpa.itc.it/markus/tmp/frank/ (also as tar.gz) Thanks Markus
comment:3 by , 19 years ago
Markus, It seems you are right. HDF-EOS documentation is not clear on this topic. I have made required changes. Andrey
comment:4 by , 19 years ago
Andrey, thanks for this! I have tried again, but still a major north-south shift appears. I am not sure if it is on my side: gdalwarp -t_srs EPSG:4326 'HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Snow_Spatial_QA' MOD_Grid_Snow_500m_Day_Tile_Snow_Cover GDAL: Auto register /usr/local/lib/gdalplugins/gdal_GRASS.so using GDALRegister_GRASS. HDF4Image: List of dimensions in grid Snow_Spatial_QA: YDim,XDim HDF4Image: Grid projection: projection code: 16, zone code -1, sphere code -1 HDF4Image: Grid geolocation: top left X 0.000000, top left Y 5559752.598333, low right X 1111950.519667, low right Y 4447802.078667, cols 2400, rows 2400 GDAL: GDALOpen(HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Snow_Spatial_QA) succeeds as HDF4Image. GDAL: GDALOpen(MOD_Grid_Snow_500m_Day_Tile_Snow_Cover) succeeds as GTiff. GDAL: GDALWarpKernel()::GWKNearestNoMasksByte() Src=0,0,2400x2400 Dst=0,0,3211x2064 :0...10...20...30...40...50...60...70...80...90...100 - done. GDAL: GDALClose(MOD_Grid_Snow_500m_Day_Tile_Snow_Cover) GDAL: GDALClose(HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Snow_Spatial_QA) GDAL: GDALDeregister_GTiff() called. Observations: sphere code -1, isn't Sinusoidal on a sphere? http://nsidc.org/data/modis/faq.html#7 Maybe this introduces the shift when warping. I have attached modis_aqua_snow_qgis_vmap2.jpg to illustrate the shift (see the islands such as Elba or corsica). Best regards Markus
by , 19 years ago
Attachment: | modis_aqua_snow_qgis_vmap2.jpg added |
---|
Still shift to be observed in MODIS Grid
comment:5 by , 19 years ago
I made a mistake in my previous fix. Please, try the new version. Spheroid code -1 is normal, that means that spheroid parameters should be taken from the HDF file, not from the predefined table. Regards, Andrey
comment:6 by , 19 years ago
Andrey, it's still shifted... I'll run make distclean and try again. Best regards, Markus
comment:7 by , 19 years ago
Markus, Now I'm agree that there is a problem with vertical shift. The problem is not in missed half pixel deltas, the shift is much bigger. I don't know where we have a bug, I will dig into this issue more deeply. Andrey
comment:8 by , 19 years ago
Markus, Are you sure that this scene has the correct georeferencing? Do you have the same problem on other MODIS scenes? I just don't know how to solve this. I just fetch the parameters from HDF and pass them as is. Andrey
comment:9 by , 19 years ago
Andrey, well, I was able to correctly geocode it with MRT. Maybe they do tricks for SIN (they did once for the ISIN)? Basically their code uses gctp (when you download MRT after registration, you get all the source code as well). I have put it here: http://mpa.itc.it/markus/tmp/MRTLinux.zip 69MB Maybe their reading code contains something specific? Sorry for all these problems! Markus
comment:10 by , 19 years ago
Hmm, I have spent even more time digging into this issue. Now I'm not sure that we have a bug in GDAL. It is likely mismatch in ellipsoids parameters. What I did step by step: $ gdalwarp -t_srs "+proj=latlong +ellps=sphere" HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover gdal_out.tif $ resample -p mrt.prm where the mrt.prm is INPUT_FILENAME = MYD10A1.A2004244.h18v04.004.2004256140521.hdf OUTPUT_FILENAME = mrt.tif RESAMPLING_TYPE = BILINEAR OUTPUT_PROJECTION_TYPE = GEO OUTPUT_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 ) and I got the coincided results. Image sizes are different, but when I opened them in OpenEV all coast lines become the same. Could you try this out? Note, that you should choose the same ellipsoid types both in gdalwarp and in MRT projection parameters. Best regards, Andrey PS. Please, use the latest CVS snapshot.
comment:11 by , 19 years ago
Andrey, I have just updated from CVS, it currently fails to compile: ../.libs/libgdal.so: undefined reference to `HDF4Dataset::GetDataTypeSize(long) const' ../.libs/libgdal.so: undefined reference to `HDF4Dataset::GetDataType(long) const' collect2: ld returned 1 exit status make[1]: *** [ogrinfo] Fehler 1 make[1]: Leaving directory `/home/neteler/soft/gdal/ogr' make: *** [ogr-apps] Fehler 2 even after 'make distclean'. Do I err on my side? Thanks, Markus
comment:13 by , 19 years ago
Andrey, compiles, works, but still something odd. When using gdalwarp for reprojecting from Sinusoidal to LatLong/WGS84 I have to use +ellps=sphere (!). Then it almost (but not perfectly) matches my WGS84 vector data. I assume that the hdf driver is using the wrong sphere definition for MODIS Sinusoidal which not the PROJ4 sphere but a "private" MODIS sphere: GRS80 Authalic Sphere is r=6370997m, but the MODIS sphere has instead a radius of 6371007.181. (see http://edcdaac.usgs.gov/landdaac/tools/mrtswath/info/ReleaseNotes.pdf). I don't know how to verify this but would like to ask you to do it... Once MODIS SIN is assigned to the correct sphere I am sure that it finally will work. Best regards Markus
comment:14 by , 19 years ago
Andrey, compiles, works, but still something odd. When using gdalwarp for reprojecting from Sinusoidal to LatLong/WGS84 I have to use +ellps=sphere (!). Then it almost (but not perfectly) matches my WGS84 vector data. I assume that the hdf driver is using the wrong sphere definition for MODIS Sinusoidal which not the PROJ4 sphere but a "private" MODIS sphere: GRS80 Authalic Sphere is r=6370997m, but the MODIS sphere has instead a radius of 6371007.181. (see http://edcdaac.usgs.gov/landdaac/tools/mrtswath/info/ReleaseNotes.pdf). I don't know how to verify this but would like to ask you to do it... Once MODIS SIN is assigned to the correct sphere I am sure that it finally will work. Best regards Markus
comment:15 by , 19 years ago
Markus, HDF driver uses the _right_ sphere. $ gdalinfo HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover ... PROJCS["unnamed", GEOGCS["Unknown datum based upon the custom spheroid", DATUM["Not specified (based on custom spheroid)", SPHEROID["Custom spheroid",6371007.181,0]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433]], PROJECTION["Sinusoidal"], PARAMETER["longitude_of_center",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["Meter",1]] ... Note that spheroid is the same as you are mentioned (6371007.181). This is fetched from the HDF file in agreement with the GCTP reference. I just can't find nothing criminal in the georeferencing code. Have you tried the reprojection scheme I described previously? I think that the source of problem in your case is wrong spheroid assumptions for different layers. Regards, Andrey
comment:16 by , 19 years ago
> From Andrey Kiselev 2005-04-19 10:23 [reply] ------- > Markus, > HDF driver uses the _right_ sphere. Sorry for the confusion. Yes, the sphere is right. However, rewarped maps result shifted. No idea, why... Below the procedure (sorry, long, but hopefully useful): # from MODIS HDF to normal sphere: gdalwarp -t_srs "+proj=latlong +ellps=sphere" \ HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover \ gdal_normal_sphere.tif gdalinfo gdal_normal_sphere.tif | head -25 Driver: GTiff/GeoTIFF Size is 3211, 2064 Coordinate System is: GEOGCS["Normal Sphere (r=6370997)", DATUM["unknown", SPHEROID["unnamed",6370997,0]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433]] Origin = (0.000000,50.000000) Pixel Size = (0.00484489,-0.00484489) Corner Coordinates: Upper Left ( 0.0000000, 50.0000000) ( 0d 0'0.01"E, 50d 0'0.00"N) Lower Left ( 0.0000000, 40.0001409) ( 0d 0'0.01"E, 40d 0'0.51"N) Upper Right ( 15.5569513, 50.0000000) ( 15d33'25.02"E, 50d 0'0.00"N) Lower Right ( 15.5569513, 40.0001409) ( 15d33'25.02"E, 40d 0'0.51"N) Center ( 7.7784756, 45.0000705) ( 7d46'42.51"E, 45d 0'0.25"N) Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray # From normal sphere to WGS84 gdalwarp -t_srs EPSG:4326 gdal_normal_sphere.tif gdal_wgs84_LL.tif gdalinfo gdal_wgs84_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.00484485,-0.00484485) Corner Coordinates: Upper Left ( 0.0000000, 50.1892281) ( 0d 0'0.01"E, 50d11'21.22"N) Lower Left ( 0.0000000, 40.1894659) ( 0d 0'0.01"E, 40d11'22.08"N) Upper Right ( 15.5568006, 50.1892281) ( 15d33'24.48"E, 50d11'21.22"N) Lower Right ( 15.5568006, 40.1894659) ( 15d33'24.48"E, 40d11'22.08"N) Center ( 7.7784003, 45.1893470) ( 7d46'42.24"E, 45d11'21.65"N) Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray #Comparing to Vector map: ogrinfo -summary polbnda_europe_wgs84.shp polbnda_europe_wgs84 | head -18 INFO: Open of `polbnda_europe_wgs84.shp' using driver `ESRI Shapefile' successful. Layer name: polbnda_europe_wgs84 Geometry: Polygon Feature Count: 4378 Extent: (-31.265751, 17.224695) - (46.725971, 80.834053) Layer SRS WKT: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.01745329251994328, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]] ... should match, right? Unfortunately it doesn't, see attachment modis_snow_vmap0_wgs84_gdal.png No idea what's wrong. ######## This is what MRT does: 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 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 Output image info ----------------- output image corners (lat/lon): UL: 49.999999995507 0.000000000000 UR: 49.999999995507 15.562266608930 LL: 39.994423054211 0.000000000000 LR: 39.994423054211 15.562266608930 output image corners (X-Y projection units): UL: 0.000000000000 49.999999995507 UR: 15.562266608930 49.999999995507 LL: 0.000000000000 39.994423054211 LR: 15.562266608930 39.994423054211 band type lines smpls pixsiz min max fill 1) Day_Tile_Snow_Cover UINT8 1698 2641 0.0059 0 255 255 The MRT result is shown in the second attachment, which looks good: modis_snow_vmap0_wgs84_MRT.png File details: gdalinfo modis_snow_vmap0_wgs84_MRT.Day_Tile_Snow_Cover.tif Driver: GTiff/GeoTIFF Size is 2641, 1698 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.002946,49.997054) Pixel Size = (0.00589257,-0.00589257) Metadata: TIFFTAG_SOFTWARE=MODIS Reprojection Tool v3.2a July 2004 Corner Coordinates: Upper Left ( 0.0029463, 49.9970537) ( 0d 0'10.61"E, 49d59'49.39"N) Lower Left ( 0.0029463, 39.9914768) ( 0d 0'10.61"E, 39d59'29.32"N) Upper Right ( 15.5652129, 49.9970537) ( 15d33'54.77"E, 49d59'49.39"N) Lower Right ( 15.5652129, 39.9914768) ( 15d33'54.77"E, 39d59'29.32"N) Center ( 7.7840796, 44.9942652) ( 7d47'2.69"E, 44d59'39.35"N) Band 1 Block=2641x1 Type=Byte, ColorInterp=Gray Find also attached the gdalinfo output of GDALwarped file and MRTwarped file. Overall, I don't want to stay on your nerves but just help to find the problem. Best regards Markus
by , 19 years ago
Attachment: | modis_snow_vmap0_wgs84_gdal.png added |
---|
MODIS gdalwarp result vs VMAP0 (both WGS84)
by , 19 years ago
Attachment: | modis_snow_vmap0_wgs84_MRT.png added |
---|
MODIS MRT program result vs VMAP0 (both WGS84)
comment:17 by , 19 years ago
Markus, What about file with parameters for transforation with MRT? > Overall, I don't want to stay on your nerves but just help to find the problem. No problm, I want to solve this poblm too. I see the shift when I use different ellipsoids and long/lat "projection". As soon as I have transformed both images to the same ellipsoids the shift goes away. Andrey
comment:18 by , 19 years ago
> I see the shift when I use different ellipsoids and long/lat "projection". As > soon as I have transformed both images to the same ellipsoids the shift goes away. Yes, but according to the gdalinfo output the ellipsoids are identical Do you see a mistake in my procedure as posted earlier today? Thanks Markus
comment:19 by , 19 years ago
I do not see any mistake, also I heve read through the RT sources and found that they are calculate bounding coordinates identically to me (anyway, the half pixel error cannot be a reason of so huge shift). Andrey
comment:20 by , 19 years ago
Well, one more experiment. gdalinfo reports about source MODIS dataset the following: Size is 2400, 2400 PROJCS["unnamed", GEOGCS["Unknown datum based upon the custom spheroid", DATUM["Not specified (based on custom spheroid)", SPHEROID["Custom spheroid",6371007.181,0]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433]], PROJECTION["Sinusoidal"], PARAMETER["longitude_of_center",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["Meter",1]] Origin = (0.000000,5559752.598333) Pixel Size = (463.31271653,-463.31271653) .... Corner Coordinates: Upper Left ( 0.000, 5559752.598) ( 0d 0'0.01"E, 50d 0'0.00"N) Lower Left ( 0.000, 4447802.079) ( 0d 0'0.01"E, 40d 0'0.00"N) Upper Right ( 1111950.520, 5559752.598) ( 15d33'26.06"E, 50d 0'0.00"N) Lower Right ( 1111950.520, 4447802.079) ( 13d 3'14.66"E, 40d 0'0.00"N) Center ( 555975.260, 5003777.338) ( 7d 4'15.84"E, 45d 0'0.00"N) Band 1 Block=2400x1 Type=Byte, ColorInterp=Gray I do $ gdalwarp -t_srs EPSG:4326 HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover gdal_wgs84.tif And now I have 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) Metadata: AREA_OR_POINT=Area Corner Coordinates: Upper Left ( 0.0000000, 50.1892278) ( 0d 0'0.01"E, 50d11'21.22"N) Lower Left ( 0.0000000, 40.1894911) ( 0d 0'0.01"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 Note, that the warped image is 11' shifted to north! Thatis exactly where we have a problem. I don't know why this happens, probably this is an error in the datum shifting code which arises in cases of transformation to unprojected scene. I will do more researches in this direction Andrey
comment:21 by , 19 years ago
Loooks like this bug should be reassigned against PROJ.4 (or I just misunderstand something). Let's try the datum shifting to geographic coordinates: $ cat coords.lst 771180 4594849 $ cs2cs +proj=sinu +ellps=sphere +to +proj=latlong +ellps=sphere coords.lst 9d14'5.33"E 41d19'20.961"N 0.000 $ cs2cs +proj=sinu +ellps=sphere +to +proj=latlong +ellps=WGS84 coords.lst 9d14'5.33"E 41d30'48.042"N 2210.028 Well, there is 11' shift we are walking around. I'm not too knowledgeable in PROJ.4 code, but I shall try to find the error. Andrey
comment:22 by , 19 years ago
Andrey mentioned in IRC that he expected the lat/long location to be the same on both ellipsoids (sphere and WGS84). Internally PROJ.4 is converting the location on the sphere into a geocentric location, and then convert it back to lat/long on WGS84. Because their shapes are quite different the angles are also noticably different (as well as the big shift in Z).
comment:23 by , 19 years ago
Markus, I have discussed this issue with Frank and what we are decided for now: 1. PROJ.4 is correct in its transformation, it properly shifts coordinates between different ellipsoids. 2. MRT is not completely correct. It should shift the image too. 3. If you want use both PROJ.4 and MRT results simultaneously you should not warp the image between datums based on different ellipsoids, but just rewrite the projection definition, i.e. $ gdalwarp -t_srs "+proj=latlong +ellps=sphere" HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover gdal_sphere.tif $ gdal_translate -a_srs EPSG:4326 gdal_sphere.tif gdal_wgs84.tif gdal_wgs84.tif now should be compatible with the MRT output. Regards, Andrey
Note:
See TracTickets
for help on using tickets.