Resample MODIS latitude/longitude at 1km resolution

[nevermind]

Hey everybody,

I am trying to resample some MODIS level-2 data (cloud effective radius from MOD06_L2.A2021299.1620.061.2021300211624) and export it as a georeferenced geoTIFF. Generally, I got it to work, but the resolution drops of my effective radius geoTIFF drops down from 1km to 5km per pixel. I supposed that  this is because the latitude/longitude datasets in the MODIS hdf-file only have a 5km resolution, so I would like to resample them to a resolution of 1km, before I use them to georeference the effective radius data. This is my code:

========================

from satpy import Scene
from pyresample.geometry import AreaDefinition

file = ["E:/Jasper/Studium/BA_Thesis/MODIS_data/MOD06_L2.A2021299.1620.061.2021300211624.hdf"]
modis = {"modis_l2": file}
channel = Scene(filenames=modis)
channel.load(["cloud_effective_radius", "latitude", "longitude"])
lats = channel["latitude"].compute().data
lons = channel["longitude"].compute().data

area_id = "Ecuador"
proj_id = "WGS84"
description = "Galapagos Archipelago WGS84"
projection = "+init=EPSG:4326"
width = lats.shape[1]
height = lats.shape[0]
area_extent = (lats.min(), lats.max(), lons.min(), lons.max())

my_area = AreaDefinition(area_id, description, proj_id, projection, width, height, area_extent)
print(my_area)

channel_res = channel.resample(my_area, method="nearest")
channel_res.save_datasets(writer="geotiff",
fill_value=0,
datasets=["cloud_effective_radius"],
filename="{name}_{start_time:%Y%m%d_%H%M%S}_georef.tif",
base_dir="E:/Jasper/Studium/BA_Thesis/MODIS_data/output")

========================

I have attached both the georeferenced and ungeoreferenced TIFFs.

Any help is greatly appreciated!

  - Jasper

[lobsiger...@gmail.com]

Jasper,

whithout understanding your code but looking at my own satpy scripts I see this:

When I plot MODIS ( l1b ) data from Aqua/Terra I have a special load statement:

...

scn = Scene(filenames = bestfiles, reader = 'modis_l1b')

scn.load([composite], resolution=1000)

...

IIRC this is the only satpy script where I indicate the resolution=xxxx.

I'm not sure this has anything to do with your problem though :-(.

Regards,

Ernst

[nevermind]

Hey Ernst,

Thanks for the idea, but I have already tried this and sadly it did not work.

- Jasper

[lobsiger...@gmail.com]

Hi Jasper,

not sure the nearest resampler has also a resolution keyword? But have you also tried reduce_data=False ?

At least this keyword sounds like your problem :-\.

Ernst

[David Hoese]

Jasper,

If the MODIS data you are providing is only available at 5km then that's
the only resolution the reader/Scene will be able to load. The
`resolution` passed to Scene.load as Ernst pointed out is only for
choosing between the resolutions available from the data files (ex. L1b
files are typically available in 250m, 500m, 1km for bands 1 and 2).

"reduce_data" is passed to Scene.resample and will not have an effect
here. Data reduction is used for gridded data (ex. geostationary) and
stops pre-resampling processing which slices away the input chunks that
Satpy knows won't appear in the output. It has no effect for swath-based
data like MODIS (unless you have a gridded L2 or L3 product which I
don't think our readers currently support).

Resampling is the only way to convert the resolution of the data to
something higher resolution *but* you are still only starting with 5km
data so it won't add any information that you don't already have. The
only other reason to do this would be to compare to other data at that
resolution, but in that case you want to resample to the Swath or Area
of that other data. In your code you are creating a new AreaDefinition
based on your existing swath extents and size which might not be the
best way to preserve resolution. You may want to try creating your
AreaDefinition by specifying the resolution in `create_area_def` instead
of the shape as described here:

https://pyresample.readthedocs.io/en/latest/geometry_utils.html

Lastly, you should make sure you update the python-geotiepoints. I just
released a new version last week that should speed up MODIS geolocation
interpolation by 30-60 seconds.

Dave

> --
> You received this message because you are subscribed to the Google
> Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pytroll+u...@googlegroups.com
> <mailto:pytroll+u...@googlegroups.com>.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/pytroll/847a4067-493a-450e-a5ab-e38e2f1a8f98n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/847a4067-493a-450e-a5ab-e38e2f1a8f98n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[nevermind]

Hi Dave,

Thanks for the advice! I have created a new area definition using this code:

=======================

from satpy import Scene
from pyresample.geometry import create_area_def

file = ["E:/Jasper/Studium/BA_Thesis/MODIS_data/MOD06_L2.A2021299.1620.061.2021300211624.hdf"]
modis = {"modis_l2": file}
channel = Scene(filenames=modis)
channel.load(["cloud_effective_radius", "latitude", "longitude"])
lats = channel["latitude"].compute().data
lons = channel["longitude"].compute().data

area_id = "Ecuador"
proj_id = "WGS84"
description = "Galapagos Archipelago WGS84"

proj_string = "+init=EPSG:4326"
area_extent = (lons.min(), lats.min(), lons.max(), lats.max())

my_area = create_area_def("Ecuador",
{'proj': 'longlat', 'datum': 'WGS84'},
area_extent=area_extent,
resolution=0.00904372,
units="degrees",
description=description)

print(my_area)

channel_res = channel.resample(my_area, method="nearest")
channel_res.save_datasets(writer="geotiff",
fill_value=0,
datasets=["cloud_effective_radius"],
filename="{name}_{start_time:%Y%m%d_%H%M%S}_georef.tif",
base_dir="E:/Jasper/Studium/BA_Thesis/MODIS_data/output")

=======================

It seems to be working and I am getting a georeferenced TIFF, but there are still some problems that I am unsure of how to solve:

1) The resolution value is just my conversion of 1km pixel size into degree. Sadly, every MODIS scene has a slightly different resolution, so I would like to find a way to import the pixel size from the metadata, just like the area_extent part, but I don't know how to do this/where the data is stored in the HDF file. I could probably make it work with just this standard value, but it would probably be better to always use the actual pixel size.

2) Using the "nearest neighbor" resampling method leads to "blurry" pixels, which makes sense due to how the algorithm is operating, but I was wondering if there was a way to fix this. Otherwise, I would be overestimating the cloud cover. I have attached both the georeferenced and ungeoreferenced TIFFs to illustrate this.

- Jasper