Hi
Solved.
The problem was that the folks who had geocoded MODIS data available at GLCF
UMIACS didn't pay enough attention to the difference between WGS84 and
SPHERE, though I don't really understand what was their error (anybody?).
Anyway, the header of these MODIS imagery is rubbish as to earth model:
[pingwin@localhost decomp]$ gdalinfo Goodes.EUAS.2000321.band3.tif
Driver: GTiff/GeoTIFF
Size is 27600, 17140
Coordinate System is:
PROJCS["unnamed",
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"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-215500.000000,8673500.000000)
Pixel Size = (500.00000000,-500.00000000)
Metadata:
AREA_OR_POINT=Point
TIFFTAG_SOFTWARE=hdf2geotiff v0.1.0
Corner Coordinates:
Upper Left ( -215500.000, 8673500.000)
(12347240d29'42949672966.25"W,496954943d36'1791001362456.19"N)
Lower Left ( -215500.000, 103500.000)
(12347240d29'42949672966.25"W,5930113d10'21474836526.57"N)
Upper Right (13584500.000, 8673500.000)
(778334516d47'2800318677035.68"E,496954943d36'1791001362456.19"N)
Lower Right (13584500.000, 103500.000)
(778334516d47'2800318677035.68"E,5930113d10'21474836526.57"N)
Center ( 6684500.000, 4388500.000) (382993638d
9'1378684502034.72"E,251442528d23'906238099491.38"N)
Band 1 Block=27600x1 Type=UInt16, ColorInterp=Gray
So when warping this one into e.g. latlong/wgs84, the sphere has to be
forced. But, which is strange and the reason of all the confussion, is has
to be forced for *both* the -s_srs and -t_srs, otherwise it is no good. See:
gdalwarp -s_srs '+proj=goode +ellps=sphere +lon_0=30E
+x_0=3335846.22854' -t_srs '+proj=latlong +ellps=sphere' -te 7 33 34 59
Goodes.EUAS.2000321.band3.tif modis_ll_wgs84.tif
The result we get:
[pingwin@localhost decomp]$ gdalinfo modis_ll_wgs84.tif
Driver: GTiff/GeoTIFF
Size is 8833, 8506
Coordinate System is:
GEOGCS["Normal Sphere (r=6370997)",
DATUM["unknown",
SPHEROID["unnamed",6370997,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (7.000000,59.000000)
Pixel Size = (0.00305663,-0.00305663)
Corner Coordinates:
Upper Left ( 7.0000000, 59.0000000) ( 7d 0'0.00"E, 59d 0'0.00"N)
Lower Left ( 7.0000000, 33.0002745) ( 7d 0'0.00"E, 33d 0'0.99"N)
Upper Right ( 33.9992447, 59.0000000) ( 33d59'57.28"E, 59d 0'0.00"N)
Lower Right ( 33.9992447, 33.0002745) ( 33d59'57.28"E, 33d 0'0.99"N)
Center ( 20.4996223, 46.0001373) ( 20d29'58.64"E, 46d 0'0.49"N)
Band 1 Block=8833x1 Type=UInt16, ColorInterp=Gray
Although gdalinfo says our output is on sphere, *actually* it is on wgs84 -
as you can see on the screendump, where the VMAP0 (projected ll/wgs84 as we
know) is overlayed and both fit perfect. The last thing remaining is to fix
the output's header with gdal_translate.
Maciek
----- Original Message ----- From: "Maciek Sieczka" <werchowyna@epf.pl>
To: "Markus Neteler" <neteler@itc.it>
Cc: <grasslist@baylor.edu>
Sent: Wednesday, April 13, 2005 10:33 PM
Subject: [GRASSLIST:6451] Re: Goode in Grass
Markus
Thank you for the links. I did some trial and error. Sorry it took so
long.
Please comment.>>I'd like to work on MODIS imagery of Poland, Eastern Europe. I
>>couldn't
>>figure out the appropriate projection definition by myself.>Most MODIS maps are in Sinusoidal (on MODIS sphere).
>Where did you get the map?say
ftp://ftp.glcf.umiacs.umd.edu/modis/500m/Eurasia/Goodes.EUAS.2000321/
via the GLCF UMIACS wwwLooking at
http://xserve.flids.com/pipermail/gdal-dev/2004-March/005095.html
"... cutting the image into pieces which I seperately reprojected using
either a
sinusoidal or a mollweide projection ..."which contains the link to
http://gis.esri.com/library/userconf/proc98/PROCEED/TO850/PAP844/P844.HTMFrom this I tried:
/*
/* Northern Eurasia
/* REGION = "02 "
4, 30.0000000, 65.0000000
-40.0000000, 40.7366111
-40.0000000, 60.0000000
-10.0000000, 60.0000000
-10.0000000, 90.0000000
180.0000000, 90.0000000
180.0000000, 40.7366111
end#####################################
dd.r02:input
projection GEOGRAPHIC
units DD
parameters
output
projection MOLLWEIDE
units METERS
spheroid SPHERE
xshift 3335846.22854
yshift -336410.83237
parameters
30 00 00
end
######################################note that we cannot cut out on the fly:
gdalwarp -s_srs '+proj=moll +ellps=sphere +lon_0=30E \
+y_0=-336410.83237 +x_0=3335846.22854' \
Goodes.EUAS.2000321.band4.tif d02_europe_EUAS.2000321.band4.moll.tifThis step seems redundant to me. Is there a reason you don't reproject to
wgs84 at once? Whether I apply the first step or not, the final wgs84
layer is the same.#Since it's only defined for a subregion, extract the valid part
#(not completely true for part of greenland, I think):
gdalwarp -s_srs '+proj=moll +ellps=sphere +lon_0=30E \
+y_0=-336410.83237 +x_0=3335846.22854' \
-t_srs '+init=epsg:4326' -te -40 40 180 90 \
d02_europe_EUAS.2000321.band4.moll.tif
d02_europe_EUAS.2000321.band4.LL.tif1. These settings you found give more-less good result. But, according to
that ESRI document, do you understand how the x,y shifts were derived?2. As we move east or west from the central meridian in Mollweide's
projection, there is an growing E distortion. How to calculate an optimal
x_0 for Mollweide's, so the E distortion is minimalized in the area of my
interest (a part of Poland)?3. I noticed that +y_o=-316410.83237 (20000 more than in this ESRI
document)
gives a better match in the y axis (when the VMAP0 "COASTLINE" is
overlayed,
go figure). Is it me, the ESRI document or this MODIS raster in error ? Or
maybe the VMAP0 (20 km?!). I don't have any other global data to compare
to.
But also when some wgs84 Landsat TM of Poland is overlayed on the MODIS
reprojected to wgs84, they fit better when "my" +y_o=-316410.83237 is
used.4. There is Goode's projection in proj 4.4.9 (+proj=goode) which takes
care
of all the "zones" in some automated way. In case of the Eurasian MODIS we
discuss here it works, as expected, similarly to +proj=moll, but:a. It handles the N and S utmost parts of my MODIS image properly while
Mollweide's introduces a vast distorion there.
b. y_0 is intended not to be provided in the proj definition for
+proj=goode
(I guess), but, following my observation of the 20km vertical shift, if
y_0=20000 is given, the resulting wgs84 MODIS fits the VMAP0 data *much*
better than when y_0 is omitted. Example:gdalwarp -s_srs '+proj=goode +ellps=sphere +lon_0=30E
+towgs84=0,0,0,0,0,0,0
+x_0=3335846.22854 +y_0=20000' -t_srs '+proj=latlong +ellps=wgs84
+towgs84=0,0,0,0,0,0,0' -te 7 33 34 59 Goodes.EUAS.2000321.band3.tif
modis_goode_esri_y0+20000.tifYours badly puzzled,
Maciek
(attachments)
