New to Satpy

297 views
Skip to first unread message

Bathobile Maseko

unread,
Oct 6, 2020, 8:40:49 AM10/6/20
to pytroll
Dear All

I have used the CF writer in Satpy to convert the SEVIRI hrit data file to NetCDF format, however it writes out the data in map units (dims: (y,x)). How can I get the data in (lon, lat) instead?

My script:


Thanking you in advance for your assistance 

Stephan Finkensieper

unread,
Oct 6, 2020, 10:45:02 AM10/6/20
to pyt...@googlegroups.com
Dear Bathobile,

SEVIRI HRIT data come in the geostationary satellite projection. This is an irregular grid, that means you cannot represent it with one-dimensional longitude and latitude vectors. Instead, the (lon, lat) coordinates of each grid cell are unique, so that you end up with two-dimensional latitude and longitude coordinates. You should find these (variables "longitude" and "latitude") in the netCDF file produced by the CF writer.

Alternatively, you can resample the scene to a regular grid in "latlong" projection before saving it. Then the projection coordinates y/x are equivalent to latitude/longitude and you can rename them accordingly. See the example below.

Hope this helps,
Stephan

import glob
import satpy
import pyresample.geometry


# Read SEVIRI data.
filenames = glob.glob('/tmp/seviri/*')
scene = satpy.Scene(filenames=filenames,
                    reader='seviri_l1b_hrit')
scene.load(['IR_108'])

# Resample to regular lat/lon grid.
area = pyresample.geometry.create_area_def(
    area_id='mygrid',
    projection={'proj': 'latlong', 'lon_0': 0},
    width=180,  # increase as needed
    height=180,
    area_extent=(-90, -90, 90, 90)
)
resampled = scene.resample(area)

# Projection coordinates y/x are now equivalent to lat/lon,
# so rename them accordingly.
resampled['IR_108'] = resampled['IR_108'].rename(
    {'y': 'lat', 'x': 'lon'})
resampled['IR_108']['lon'].attrs = {'units': 'degrees_east',
                                    'standard_name': 'longitude'}
resampled['IR_108']['lat'].attrs = {'units': 'degrees_north',
                                    'standard_name': 'latitude'}

# Save to netCDF.
resampled.save_datasets(
    datasets=['IR_108'],
    writer='cf',
    filename='test.nc',
    flatten_attrs=True,
    exclude_attrs=['raw_metadata'],
    include_lonlats=False  # note this one here!
)



--
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.
To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/25a0fb92-2d68-4533-aaab-f67b8fdaaec6n%40googlegroups.com.

David Hoese

unread,
Oct 6, 2020, 11:19:06 AM10/6/20
to pyt...@googlegroups.com
Bathobile,

I think the main takeaway here is that you should already have longitude
and latitude variables in your NetCDF file. Is this not the case?

Dave
>     filename='test.nc <http://test.nc>',
>     flatten_attrs=True,
>     exclude_attrs=['raw_metadata'],
>     include_lonlats=False  # note this one here!
> )
> |
>
>
>
> On Tue, Oct 6, 2020 at 2:40 PM Bathobile Maseko <bmase...@gmail.com
> <mailto:bmase...@gmail.com>> wrote:
>
> Dear All
>
> I have used the CF writer in Satpy to convert the SEVIRI hrit data
> file to NetCDF format, however it writes out the data in map units
> (dims: (y,x)). How can I get the data in (lon, lat) instead?
>
> My script:
>
>
> Thanking you in advance for your assistance
>
> --
> 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>.
> <https://groups.google.com/d/msgid/pytroll/25a0fb92-2d68-4533-aaab-f67b8fdaaec6n%40googlegroups.com?utm_medium=email&utm_source=footer>.
>
> --
> 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/CACn-jjZ-3aRfCO0qcbqz1AZjRsLwEhamjSiJXe0f7cpVt99Z8g%40mail.gmail.com
> <https://groups.google.com/d/msgid/pytroll/CACn-jjZ-3aRfCO0qcbqz1AZjRsLwEhamjSiJXe0f7cpVt99Z8g%40mail.gmail.com?utm_medium=email&utm_source=footer>.

Bathobile Maseko

unread,
Oct 12, 2020, 3:36:08 AM10/12/20
to pytroll
Hi Stephan,

The netCDF file produced by the CF writer did have longitudes and latitudes. However, they were not in a manner that made sense to me, as I was getting a lot of INFINITY "values" which I didn't know what it was or what it meant. So I thought that the problem might be due to the dimensions being in map units (y,x).

Below is a screenshot of what I mean:


Resampling the scene to a regular grid in "latlong" projection as you suggested really helped me.  Now my dimensions are in the latitude and longitude coordinates, and my lat and lon make sense according to my area extent. 
 

Thank you for your assistance.

Stephan Finkensieper

unread,
Oct 14, 2020, 4:31:16 AM10/14/20
to pyt...@googlegroups.com
Hi Bathobile,

glad you found a solution! The infinity values belong to pixels where the instrument "looks" into space and therefore don't have a valid geo-location. However, these fill values don't match the _FillValue attribute of the netCDF variable (NaN). I'll fix that. Thanks for the report!

Stephan

Reply all
Reply to author
Forward
0 new messages