search

openclimatefix / satellite_image_processing

Apache License 2.0
12 stars 3 forks source link

Reprojecting other channels #8

Open JackKelly opened 4 years ago

JackKelly commented 4 years ago

Consider using pangeo-pyinterp. Might be faster than PyResample?

Some discussion of reprojection:

dfulu commented 4 years ago

Hey @JackKelly,

I'm trying to extend the projection here to the non-HRV channels.

In cell 6 of convert_native_to_reprojected_netcdf.ipynb there are some hard coded numbers that are used to calculate the ground control points for the HRV channel. I think these may be different for the non-HRV channels. You comment that they came from scene['HRV'].area where

>>> scene['HRV'].area
Area ID: geos_seviri_hrv
Description: SEVIRI high resolution channel area
Projection ID: seviri_hrv
Projection: {'a': '6378169', 'b': '6356583.8', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 5568
Number of rows: 4176
Area extent: (-5569248.1217, 1394687.3495, 5566247.7186, 5569248.1217)

Can you explain a little more how these relate to each other as I don't see a perfect mapping?

ps. for convenience here are the contents of that cell. 

TYPICAL_EAST_COLUMN_PLANNED = 2733

def get_ground_control_points_for_seviri_hrv(
    src_height: int, 
    src_width: int, 
    east_column_planned: int=TYPICAL_EAST_COLUMN_PLANNED) -> List:
    """
    Args:
        east_column_planned: Taken from the Native file header:
          header['15_DATA_HEADER']['ImageDescription']['PlannedCoverageHRV']['LowerEastColumnPlanned']
          The HRV channel on Meteosat Second Generation satellites doesn't scan the full number of columns.
          The east boundary of the HRV channel changes (e.g. to maximise the amount of the image which
          is illuminated by sunlight.
    """
    x_offset = 1000 + ((TYPICAL_EAST_COLUMN_PLANNED - east_column_planned) * 1000)
    y_offset = 500

    left =  -2736372.0236 + x_offset
    right =  2832376.1893 + x_offset
    bottom = 1394687.3495 + y_offset # from scene['HRV'].area
    top =    5570248.4773 + y_offset

    top_left = GroundControlPoint(row=0, col=0, x=left, y=top, id='top_left')
    bottom_left = GroundControlPoint(row=src_height, col=0, x=left, y=bottom, id='bottom_left')
    bottom_right = GroundControlPoint(row=src_height, col=src_width, x=right, y=bottom, id='bottom_right')
    top_right = GroundControlPoint(row=0, col=src_width, x=right, y=top, id='top_right')

    ground_control_points = [top_left, bottom_left, bottom_right, top_right]
    return ground_control_points
JackKelly commented 4 years ago

Yeah, sorry that's not better commented! IIRC, it took me a while to find the right numbers! I think I used the area extent from the HRIT version of HRV (see compare_EUMETSAT_file_formats.ipynb cell 20), and manually tweaked the numbers a little until the image aligned with the coastlines plot! I can send you some HRIT files if that's useful (I think I downloaded HRIT files with all channels).

dfulu commented 4 years ago

That could be useful thanks