[GRASSLIST:6860] Re: Goode in Grass

Mon May 23 17:47:19 EDT 2005

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 at 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 at 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 at epf.pl>
To: "Markus Neteler" <neteler at itc.it>
Cc: <grasslist at 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 www
>>
>> Looking 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.HTM
>>
>>>From 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.tif
>
> This 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.tif
>
> 1. 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.tif
>
> Yours badly puzzled,
> Maciek
-------------- next part --------------
A non-text attachment was scrubbed...
Name: goodes_sphere.png
Type: image/png
Size: 124294 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20050523/45a20051/goodes_sphere.png