[gdal-dev] MODIS reprojection

[laura0]

Hi,
I'm working with GDAL tool in order to use MODIS (MOD14.A2 and MYD14.A2)
data but I have problems with reprojection in UTM 32. The lat/lon pixel
coordinates I get are not correct (similar but a bit shifted).
This is what I have done:
In order to geolocate swaths, there is a separate product that contains two
arrays, Latitude and Longitude. Each row, column in these arrays contains
the lat and lon for each corresponding data pixel. I have used gdal to get
from this product, the following information:
Coordinate System is `'
GCP[ 0]: Id=, Info=
(0.5,0.5) -> (36.5436744689941,25.0500869750977,0)
GCP[ 1]: Id=, Info=
(123.5,0.5) -> (32.4521675109863,24.7417182922363,0)
GCP[ 2]: Id=, Info=
....
....
GCP[143]

and then I added ALL these rows as option to the command gdal_translate -of
Gtiff -co "tfw=yes" -gcp 0.5 0.5 36.5436744689941 25.0500869750977 -gcp
123.5 0.5 32.4521675109863 24.7417182922363 ...
MOD14.A2013078.2035.005.NRT.hdf_001_0.gaiag.tif temp.tif

with gdalwarp -t_srs "+proj=utm +zone=32 +datum=WGS84" temp.tif output.tif I
get an image that looks ok but the lat/lon coordinates are still not precise
(particularly for longitude).

Can someone help me?
Thanks

Laura



--
View this message in context: http://osgeo-org.1560.x6.nabble.com/MODIS-reprojection-tp5047700.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.

[Rutger]

Hey Laura,

Since the MxD14.A2 products are already projected in the sinusoidal
projection there is no need to warp using GCP's. Fitting a polynomial
through some GCP's will always be less accurate then using the projection
the data is already in.

Somewhere in the metadata in the HDF there should be a boundingbox (corner
coordinates, in sinusoidal projection) for the specific tile you are using,
you can then assign these with the projection using gdal_translate. Once
assigned, you can use gdalwarp in the same way as you already described.

Youre gdal_translate command should look somewhat like this:
gdal_translate -of "VRT" -a_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=6371007.181 +units=m +no_defs" -a_ullr 2223901.039333
-3335851.559000 3335851.559000 -4447802.078667 <input HDF subdataset>
<output vrt file>

You can get the correct subdataset reference by running gdalinfo on you HDF.
The coordinates i used are for tile h20v12, so make sure you use the ones
from your tile. I used a vrt as output which you can use as an input for
gdal_translate, but its also fine to use Geotiff or anything else you
prefer.

Hope this helps.

Regards,
Rutger










--
View this message in context: http://osgeo-org.1560.x6.nabble.com/MODIS-reprojection-tp5047700p5047887.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.

[Etienne Tourigny]

I use a similar command, without the -a_ullr argument

e.g.

gdal_translate -a_srs '+proj=sinu +R=6371007.181 +nadgrids=@null +wktext'
<input HDF subdataset> <output tiff file>

Post by Rutger
Hey Laura,
Since the MxD14.A2 products are already projected in the sinusoidal
projection there is no need to warp using GCP's. Fitting a polynomial
through some GCP's will always be less accurate then using the projection
the data is already in.
Somewhere in the metadata in the HDF there should be a boundingbox (corner
coordinates, in sinusoidal projection) for the specific tile you are using,
you can then assign these with the projection using gdal_translate. Once
assigned, you can use gdalwarp in the same way as you already described.
gdal_translate -of "VRT" -a_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=6371007.181 +units=m +no_defs" -a_ullr 2223901.039333
-3335851.559000 3335851.559000 -4447802.078667 <input HDF subdataset>
<output vrt file>
You can get the correct subdataset reference by running gdalinfo on you HDF.
The coordinates i used are for tile h20v12, so make sure you use the ones
from your tile. I used a vrt as output which you can use as an input for
gdal_translate, but its also fine to use Geotiff or anything else you
prefer.
Hope this helps.
Regards,
Rutger
--
http://osgeo-org.1560.x6.nabble.com/MODIS-reprojection-tp5047700p5047887.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20130418/3593f26d/attachment-0001.html>

[Rutger]

Hey,

Post by Etienne Tourigny
I use a similar command, without the -a_ullr argument

And does it then capture the correct corner coordinates? If i run your
command it does assign the sinusoidal projection, but takes the pixel
coordinates (eg 0-2048) as the corner coordinates instead of the sinusoidal
ones.

@Laura,
Indeed we are talking about different products, sorry for the confusion. You
mentioned in your first email the products MYD14.A2 and MOD14.A2. I thought
you meant the "MYD14 A2" product, but you attached a "MYD14" product, which
is different. Its obvious where the confusion comes from. :)

These are all the fire products:
MYD14 Aqua Thermal Anomalies & Fire Swath 1000m 5 min
MOD14A2 Terra Thermal Anomalies & Fire Tile 1000m 8 day
MYD14A2 Aqua Thermal Anomalies & Fire Tile 1000m 8 day
MOD14A1 Terra Thermal Anomalies & Fire Tile 1000m Daily
MYD14A1 Aqua Thermal Anomalies & Fire Tile 1000m Daily
MOD14 Terra Thermal Anomalies & Fire Swath 1000m 5 min

The "MOD14" & "MYD14" products are indeed swath, and should be used together
with the MOD03 with the same timestamp.

I think the best way to process them is by using "MRT Swath" which you can
download at:
https://lpdaac.usgs.gov/tools

MRT Swath is way beyond the scope of this mailing list, but feel free to
send me a direct email about it if you have questions.

I have used Pyresample as an alternative which works very well, but requires
the use of Python:
http://code.google.com/p/pyresample/

To end with an actual GDAL related solution, there is an option to warp an
image using gdalwarp and so-called "geolocation arrays". I can get it to
work with the MODIS swaths, but i dont like the results at all. It is
previously discussed at:
http://osgeo-org.1560.x6.nabble.com/gdal-dev-reproject-python-numpy-binary-swath-lat-lon-td4978609.html

The MODIS swath products arent equally gridded and still contain what is
called the 'bow-tie' effect. (overlapping scans at high zenith angles). I
dont know if the GDAL geoloc option is supposed to handle this properly. The
image a posted earlier in the thread mentioned above clearly shows artifacts
which i cannot explain or resolve.


Regards,
Rutger








--
View this message in context: http://osgeo-org.1560.x6.nabble.com/MODIS-reprojection-tp5047700p5048019.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.

[Etienne Tourigny]

Post by Rutger
Hey,

Post by Etienne Tourigny
I use a similar command, without the -a_ullr argument

And does it then capture the correct corner coordinates? If i run your
command it does assign the sinusoidal projection, but takes the pixel
coordinates (eg 0-2048) as the corner coordinates instead of the sinusoidal
ones.

Actually, I used this with the level-3 tiled products such as MCD12Q1 or
MCD45A1. Never tried with level-2 or swath products.
Yes the corners are fine, I can convert a few of them to WGS84 and then
stitch them together with gdalwarp.

Etienne

Post by Rutger
@Laura,
Indeed we are talking about different products, sorry for the confusion. You
mentioned in your first email the products MYD14.A2 and MOD14.A2. I thought
you meant the "MYD14 A2" product, but you attached a "MYD14" product, which
is different. Its obvious where the confusion comes from. :)
MYD14 Aqua Thermal Anomalies & Fire Swath 1000m 5 min
MOD14A2 Terra Thermal Anomalies & Fire Tile 1000m 8
day
MYD14A2 Aqua Thermal Anomalies & Fire Tile 1000m 8
day
MOD14A1 Terra Thermal Anomalies & Fire Tile 1000m
Daily
MYD14A1 Aqua Thermal Anomalies & Fire Tile 1000m
Daily
MOD14 Terra Thermal Anomalies & Fire Swath 1000m 5 min
The "MOD14" & "MYD14" products are indeed swath, and should be used together
with the MOD03 with the same timestamp.
I think the best way to process them is by using "MRT Swath" which you can
https://lpdaac.usgs.gov/tools
MRT Swath is way beyond the scope of this mailing list, but feel free to
send me a direct email about it if you have questions.
I have used Pyresample as an alternative which works very well, but requires
http://code.google.com/p/pyresample/
To end with an actual GDAL related solution, there is an option to warp an
image using gdalwarp and so-called "geolocation arrays". I can get it to
work with the MODIS swaths, but i dont like the results at all. It is
http://osgeo-org.1560.x6.nabble.com/gdal-dev-reproject-python-numpy-binary-swath-lat-lon-td4978609.html
The MODIS swath products arent equally gridded and still contain what is
called the 'bow-tie' effect. (overlapping scans at high zenith angles). I
dont know if the GDAL geoloc option is supposed to handle this properly. The
image a posted earlier in the thread mentioned above clearly shows artifacts
which i cannot explain or resolve.
Regards,
Rutger
--
http://osgeo-org.1560.x6.nabble.com/MODIS-reprojection-tp5047700p5048019.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20130418/0115def9/attachment.html>

[Jonas Eberle]

Hi Rutger,

Post by Rutger
I have used Pyresample as an alternative which works very well, but requires
http://code.google.com/p/pyresample/

do you have an example how to reproject MODIS Swath data with this pyresample software?

Thanks in adavance.

Regards,
Jonas

[lauramail]

Hi Rutger,
thank you for your reply.
The metadata in the HDF show me these bounding box:
EASTBOUNDINGCOORDINATE=36.5436731054265
NORTHBOUNDINGCOORDINATE=43.0619355408597
SOUTHBOUNDINGCOORDINATE=21.8438333503548
WESTBOUNDINGCOORDINATE=6.83616802377613

and the first test has been:
gdal_translate -of Gtiff -co "tfw=yes" -a_ullr WESTBOUNDINGCOORDINATE
NORTHBOUNDINGCOORDINATE EASTBOUNDINGCOORDINATE SOUTHBOUNDINGCOORDINATE input.
tif temp.tif
gdalwarp -t_srs "+proj=utm +zone=32 +datum=WGS84" temp.tif output.tif

but with bad results. Then I read that in order to geolocate swaths, there is
a separate product called MOD03 with lat/lon matrix and I tried as my previous
post. It is not clear for me now which values should I put here (I copy from
your reply)
" +a=6371007.181 +b=6371007.181 +units=m +no_defs" -a_ullr 2223901.039333
-3335851.559000 3335851.559000 -4447802.078667"
your coordinates in option -a_ullr are totally different from mine, I guess
you have them in sinusoidal projection but I do not see anything similar.

Regards,
Laura

----Messaggio originale----
Da: kassies at gmail.com
Data: 18/04/2013 11.10
A: <gdal-dev at lists.osgeo.org>
Ogg: Re: [gdal-dev] MODIS reprojection
Hey Laura,
Since the MxD14.A2 products are already projected in the sinusoidal
projection there is no need to warp using GCP's. Fitting a polynomial
through some GCP's will always be less accurate then using the projection
the data is already in.
Somewhere in the metadata in the HDF there should be a boundingbox (corner
coordinates, in sinusoidal projection) for the specific tile you are using,
you can then assign these with the projection using gdal_translate. Once
assigned, you can use gdalwarp in the same way as you already described.
gdal_translate -of "VRT" -a_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=6371007.181 +units=m +no_defs" -a_ullr 2223901.039333
-3335851.559000 3335851.559000 -4447802.078667 <input HDF subdataset>
<output vrt file>
You can get the correct subdataset reference by running gdalinfo on you HDF.
The coordinates i used are for tile h20v12, so make sure you use the ones
from your tile. I used a vrt as output which you can use as an input for
gdal_translate, but its also fine to use Geotiff or anything else you
prefer.
Hope this helps.
Regards,
Rutger
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/MODIS-

reprojection-tp5047700p5047887.html

Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[Rutger]

Hey Laura,

The coordinates you found seem to be geographic (latitude/longitude), not
sinusoidal.

The MOD03 file is still unprojected, sometimes referred to as a swath or
granule. It contains the original sensor scans which dont yet have a defined
projection. It should be used only with other, lower level, swath products
like the MOD02HKM, MOD02QKM etc. I cant be used together with your MxD14A2
product.

The table on this website shows the 'raster type' for most MODIS Land
products:
https://lpdaac.usgs.gov/products/modis_products_table

Your higher level product is already projected (Sinusoidal tiles). I have
downloaded a MYD14A2 hdf and somewhere (with gdalinfo towards the end) there
should be two tags like:

UpperLeftPointMtrs=(1111950.519667,5559752.598333)
LowerRightMtrs=(2223901.039333,4447802.078667)

These are for tile h19v04. If you mention which tile you are using, or point
me to a file (on ftp://e4ftl01.cr.usgs.gov/MOLA/MYD14A2.005/), i wouldnt
mind looking it up for you. Although for your own confidence it would good
if you can replicate it yourself.

Let me know how it works out.

Regards,
Rutger





--
View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-R-Re-MODIS-reprojection-tp5047949p5047956.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.

[laura0]

If I try your command I get
ERROR 1: Unable to compute a transformation between pixel/line
and georeferenced coordinates for temp.tif.There is no affine transformation and no GCPs.


----Messaggio originale----

Da: ml-node+s1560n5047999h80 at n6.nabble.com

Data: 18/04/2013 16.18

A: "laura0"<lauramail at iol.it>

Ogg: Re: MODIS reprojection





I use a similar command, without the -a_ullr argument
e.g.
gdal_translate -a_srs '+proj=sinu +R=6371007.181 +nadgrids=@null +wktext' <input HDF subdataset> <output tiff file>

On Thu, Apr 18, 2013 at 6:10 AM, Rutger <[hidden email]> wrote:

Hey Laura,



Since the MxD14.A2 products are already projected in the sinusoidal

projection there is no need to warp using GCP's. Fitting a polynomial

through some GCP's will always be less accurate then using the projection

the data is already in.



Somewhere in the metadata in the HDF there should be a boundingbox (corner

coordinates, in sinusoidal projection) for the specific tile you are using,

you can then assign these with the projection using gdal_translate. Once

assigned, you can use gdalwarp in the same way as you already described.



Youre gdal_translate command should look somewhat like this:

gdal_translate -of "VRT" -a_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0

+a=<a href="tel:[phone]" value="[phone]">6371007.181 +b=<a href="tel:[phone]" value="[phone]">6371007.181 +units=m +no_defs" -a_ullr 2223901.039333

-3335851.559000 3335851.559000 -4447802.078667 <input HDF subdataset>

<output vrt file>



You can get the correct subdataset reference by running gdalinfo on you HDF.

The coordinates i used are for tile h20v12, so make sure you use the ones

from your tile. I used a vrt as output which you can use as an input for

gdal_translate, but its also fine to use Geotiff or anything else you

prefer.



Hope this helps.



Regards,

Rutger





















--

View this message in context: http://osgeo-org.1560.x6.nabble.com/MODIS-reprojection-tp5047700p5047887.html

Sent from the GDAL - Dev mailing list archive at Nabble.com.

_______________________________________________

gdal-dev mailing list

[hidden email]

http://lists.osgeo.org/mailman/listinfo/gdal-dev




_______________________________________________

gdal-dev mailing list

[hidden email]

http://lists.osgeo.org/mailman/listinfo/gdal-dev










If you reply to this email, your message will be added to the discussion below:
http://osgeo-org.1560.x6.nabble.com/MODIS-reprojection-tp5047700p5047999.html



To unsubscribe from MODIS reprojection, click here.

NAML








--
View this message in context: http://osgeo-org.1560.x6.nabble.com/R-Re-MODIS-reprojection-tp5048000.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20130418/0e5e8b4b/attachment.html>