Creating area definition by arrays of latitudes and longitudes

[anikfal]

Description: We have generated a PNG composite image and want to overlay the coast lines on this PNG image by using pycoast.

Question: The latitudes and longitudes of the pixels of this PNG image could be taken as numpy arrays from its source NetCDF file. How could we use the lat/long arrays in create_area_def to create area definition?

sample Code:

from pyresample import create_area_def
area_def = create_area_def('my_area',
...                            {'proj': 'latlong', 'lon_0': 0},
...                            area_extent=[-180, -90, 180, 90],
...                            resolution=1,
...                            units='degrees',
...                            description='Global 1x1 degree lat-lon grid')

I expect that by using lat/long arrays, some other parameters such as area_extent and resolution would be of no use.

[djhoese]

Pycoast requires that things be gridded (uniformly spaced). There are some cases like polar-orbiting satellite data where this is not the case and we describe this data with longitude/latitude 2D arrays in a pyresample SwathDefinition. SwathDefinitions do not work with pycoast.

An array of longitude/latitude that are uniformly spaced in lon/lat space (evenly spaced degrees) is a special case of the AreaDefinition object. You can use the corner points of the lon/lat arrays to define your area_extent in degrees and set the proj_dict as you have above in your example code ({'proj': 'latlong', 'datum': 'WGS84'}). I mentioned in another issue that area extent is usually in meters but for the latlong projection it should be in degrees.

So you could either create your AreaDefinition directly or use create_area_def. You know the size of the area (the shape of the image) and you have the center coordinates of each pixel. You can use those center coordinates to determine what the area extent (our edge of the pixels instead of the center points) should be.

IF your image is a projected image then you will need to know what that projection is and use it in create_area_def. create_area_def will then do the work for you to transform the lon/lat extents to coordinates in the projection space.

[anikfal]

SwathDefinitions do not work with pycoast. IF your image is a projected image then you will need to know what that projection is and use it in create_area_def.

As far as I understood, images with curvilinear grid (2D lat/lon) which are not projected, could not be overlaid by pycoast. So I think images from MSG SEVIRI data are considered as an unprojected swath, and hence could not be overlaid by pycoast.

[djhoese]

Let's see what @mraspaud and @pnuu say, they have more experience with SEVIRI data (I'm in the US).

[djhoese]

Talked with @adybbroe a little bit and he verified that SEVIRI data is projected to a geostationary projection/grid so as long as you have the definition for that projection it should work fine.

[pnuu]

MSG/SEVIRI data are in a well defined grid in geos projection, so the overlaying should work.

If at all possible in your workflow, you could add the grid directly when saving from satpy by defining the configuration as a dictionary. As an example, this adds 10x10 degree lon/lat grid over the saved image: scn.save_datasets(..., overlay={'grid': {'major_lonlat': (10, 10), 'write_text': False}}).

Unfortunately I haven't used this functionality. You can find more options to pass from pycoast documentation: https://pycoast.readthedocs.io/en/latest/

[anikfal]

Thank you for the explanations. The geos projection is supported by PROJ Software, and I hope I can overlay SEVIRI data by pycoast without much difficulty.

[sjoro]

No need to create AreaDefinition for SEVIRI. Prime service (seviri_0deg) and IODC (seviri_iodc) definitions can be found in the areas.yaml-file. So you could use get_area_def-method from satpy.resamle-module to easily get the definition you need.

[anikfal]

@sjoro Would you please add a few lines of codes?

[pnuu]

I hope this works. I'm on my phone and can't test.

from satpy.resample import get_area_def

adef = get_area_def('seviri_0deg')

EDIT: fix the import

[anikfal]

I forgot to mention that my source MSG/SEVIRI files were subset domains of the main data. Anyway, the code below can overlay the coastlines over the RGB image (PNG file) of the SEVIRI subset data:

from PIL import Image
from pycoast import ContourWriterAGG
img = Image.open('colorized_ir_clouds_20190401_120010.png')  #RGB image from MSG/SEVIRI data
proj4_string = '+proj=geos +a=6378169.0 +b=6356583.800000001 +h=35785831.0 +lon_0=41.5 +units=m +type=crs'
area_extent = (-663100, 1650000, 2583000, 4675000)
area_def = (proj4_string, area_extent)
cw = ContourWriterAGG('/home/ah/MSG8/gshhg-shp-2.3.7')
cw.add_coastlines(img, area_def, resolution='l', level=4, outline='red', width=3)
img.show()

I think seviri_0deg and seviri_iodc are usable for the complete SEVIRI data, not the subsets. The values of the area_extent of the PNG file could be found from the metadata of its source NetCDF file. The scales and meaning of area_extent are now clear to me, but I'm looking for a document with thorough explanations about the calculation of area_extent and/or geostationary projection.

[djhoese]

@anikfal Does this help at all: https://proj.org/operations/projections/geos.html

The X/Y coordinates in the 'geos' projection are in meters. The extents are the outer edges of the Area in this X/Y plane.

[anikfal]

@djhoese Thank you. Seems to be a good document.

[sjoro]

@anikfal from satpy.resample import get_area_def adef = get_area_def('seviri_0deg')