User-defined area is not positioned correctly (x_0and y_0 stere projection parameters ignored in SIFT scene graph, but still passed to Satpy for resampling)

[nedelceo]

When I try to plot satellite data using custom area definition, it is shifted to the east. But when I use Satpy without SIFT, it works as expected.

I was processing VIIRS data with SIFT_2.0.0b0 on a Windows machine.

I tried to run Satpy in the same environment as SIFT. To do this, I run python from the edited SIFT.bat file.

Area definition in areas.yaml:

EU_1160x800:
  description: EU
  projection:
    proj: stere
    lat_0: 90
    lat_ts: 90
    lon_0: 2.5
    k: 1
    x_0: 4235947.983429
    y_0: 1861952.316054
    datum: WGS84
    units: m
  shape:
    width: 1160
    height: 800
  area_extent:
    lower_left_xy: [0.000000, -5818127.123200]
    upper_right_xy: [8467532.711400, 0.000000]

SIFT scearnshot:

Image generated by Satpy:

Edited SIFT.bat:

@echo off
REM Initialize SIFT installation if necessary and run SIFT

set base_dir=%~p0

REM Activate the conda environment
call %base_dir%Scripts\activate

REM Create a signal file that we have run conda-unpack
set installed=%base_dir%.installed
if not exist "%installed%" goto install

set /p install_dir=< %installed%
if not %base_dir% == %install_dir:~0,-1% goto install

REM goto run_sift

:install
  echo Running one-time initialization of SIFT installation...
  conda-unpack
  echo %base_dir% > %installed%

REM :run_sift

REM echo Running SIFT...

REM python -m uwsift %*


pip install pycoast
python

cmd /k

[djhoese]

My guess is at least one of these applications is ignoring the x_0 and y_0 in your projection. Could you try maybe changing the x_0 and see if that shift is reflected in SIFT? If not, then could you try removing the x_0 and y_0 and adjusting your extents to match the shift they were doing?

[nedelceo]

You're right, changing or removing x_0 and y_0 has no effect. Do you have any idea how to easily but to precisely change extents to compensate for the missing x_0 and y_0.

[djhoese]

I forget if x_0/y_0 are added or subtracted from the coordinates, but you should be able to find out relatively easily. So you have extents:

    lower_left_xy: [0.000000, -5818127.123200]
    upper_right_xy: [8467532.711400, 0.000000]

So you'd take the 0.0 (an x coordinate) and subtract x_0 (4235947.983429) and the same for the other x coordinate 8467532.711400. Then do the same (subtract) for the y coordinates but with y_0. My guess is SIFT is not using x_0/y_0 on the map (because I'm guessing I was lazy and chose to ignore them) and Satpy resampling is.

[nedelceo]

I subtracted this as you advised. At first glance it seemed to work. The area with no data seemed fine. But when I loaded the data, I got image bellow. Data was placed to the same wrong place.

[djhoese]

The resampled data may be cached. I'm not sure. @ameraner would know better, but he may not be working today. If you have another data case I would guess that that should work (hopefully). Or, there may be a flag you can provided to clear out the cache. --clear-workspace maybe?

[nedelceo]

Unfortunately even with another data, the result is the same.

[ameraner]

--clear-workspace will indeed clear the cache - but cached data arrays will be generated only when the inventory/caching is activated (using the config keyword

storage:
    use_inventory_db: True

(see https://sift.readthedocs.io/en/latest/configuration/storage.html )

This is not the current default, so this is probably not the reason for the observed issue.... as just confirmed by using different data while I was writing this reply

[ameraner]

I could reproduce the issue; however, after applying the subtraction, and setting x_0 and y_0 to 0, the issue is solved. @nedelceo could it be that you forgot to set x_0 and y_0 to 0 (or to remove them) after the subtraction? I think that x_0 and y_0 are ignored for the map visualisation in SIFT, but they are still passed to satpy/pyresample for the resampling (it's easy to get tricked/confused here).

This is what I get with some AVHRR data:

area:

EU_1160x800:
  description: EU
  projection:
    proj: stere
    lat_0: 90
    lat_ts: 90
    lon_0: 2.5
    k: 1
    x_0: 0
    y_0: 0
    datum: WGS84
    units: m
  shape:
    width: 1160
    height: 800
  area_extent:
    lower_left_xy: [-4235947.983429, -7680079.439254001]
    upper_right_xy: [4231584.727971001, -1861952.316054]

[nedelceo]

Great, that is it. I forgot to removed x_0 and y_0. Thank you both.

[nedelceo]

Is it possible to do some similar subtraction trick for area defined below, if there are no x_0 and y_0?

CE_1160x800:
  description:   Central Europe 1160x800
  projection:
    proj: merc
    lat_0: 50
    lon_0: 14.5
    ellps: WGS84  
  shape:
    width: 1160
    height: 800  
  area_extent:
    lower_left_xy: [-1805152.04,  5172213.35]
    upper_right_xy: [1811322.74,  7660679.05] 

[ameraner]

Hmm.. trying this area I only see a black screen (I suppose this is also what you see?) Resampling also doesn't seem to work correctly.. I'll open a separate issue for this, I suspect this is a different problem.

[ameraner]

Issue for the mercator thing: #374
I'm also renaming this issue to point to the found problem.