output images are in unconvential south-north orientation

[johntruckenbrodt]

... causing incompatibilities when working with GDAL.

current:

Pixel Size = (10.000000000000000,10.000000000000000)
...
Corner Coordinates:
Upper Left  (  600000.000, 4890240.000) ( 10d15' 2.22"E, 44d 9'30.67"N)
Lower Left  (  600000.000, 5000040.000) ( 10d16'19.35"E, 45d 8'48.31"N)
Upper Right (  709800.000, 4890240.000) ( 11d37'21.87"E, 44d 8' 6.90"N)
Lower Right (  709800.000, 5000040.000) ( 11d40' 3.44"E, 45d 7'21.62"N)
Center      (  654900.000, 4945140.000) ( 10d57'11.72"E, 44d38'34.50"N)

expected:

Pixel Size = (10.000000000000000,-10.000000000000000)
...
Corner Coordinates:
Upper Left  (  600000.000, 5000040.000) ( 10d16'19.35"E, 45d 8'48.31"N)
Lower Left  (  600000.000, 4890240.000) ( 10d15' 2.22"E, 44d 9'30.67"N)
Upper Right (  709800.000, 5000040.000) ( 11d40' 3.44"E, 45d 7'21.62"N)
Lower Right (  709800.000, 4890240.000) ( 11d37'21.87"E, 44d 8' 6.90"N)
Center      (  654900.000, 4945140.000) ( 10d57'11.72"E, 44d38'34.50"N)

Temporary workaround: swap the pixels and coordinates using gdalwarp:

gdalwarp gt.tif gt_swp.tif

[alexamici]

@johntruckenbrodt The projection and orientation of the output image should be identical to the ones of the input DEM, so my guess is that your input DEM has the south-north orientation.

If you confirm my guess I intend to keep the behaviour as it is and mark the issue wontfix. In the case the output has a different orientation than the input, then I need more info on the two data, like the output of gdalinfo on the input DEM and the output image.

[priit111]

Here is gdalinfo of a processed s1 GRD GTC

Driver: GTiff/GeoTIFF
Files: /home/priit/PycharmProjects/run_sarsen/saaremaa_grd_s1.tif
Size is 6000, 6000
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (20.000000000000000,55.000000000000000)
Pixel Size = (0.000833333333333,0.000833333333333)
Metadata:
  AREA_OR_POINT=Area
  ascending_node_time=2025-04-07T15:57:18.435758
  azimuth_pixel_spacing=10.0
  azimuth_time_interval=0.001499360625527246
  family_name=SENTINEL-1
  geospatial_bounds=POLYGON((24.57285511361589 58.15520719687188,24.13074328927226 59.64648688336036,19.5124495010717 59.20790100848891,20.1447481025113 57.72332109659177,24.57285511361589 58.15520719687188))
  geospatial_lat_max=59.64648688336036
  geospatial_lat_min=57.72332109659177
  geospatial_lon_max=24.57285511361589
  geospatial_lon_min=19.5124495010717
  image_slant_range_time=0.005334924041781893
  incidence_angle_mid_swath=38.55607134423448
  long_name=beta nought
  mission_data_take_id=475925
  mode=IW
  number=A
  orbit_number=58653
  pass=ASCENDING
  product_first_line_utc_time=2025-04-07T16:13:06.028060
  product_last_line_utc_time=2025-04-07T16:13:31.026900
  product_type=GRD
  radar_frequency=5.40500045433435
  range_pixel_spacing=10.0
  range_sampling_rate=64345238.12571428
  relative_orbit_number=131
  start_time=2025-04-07T16:13:06.028060
  stop_time=2025-04-07T16:13:31.026900
  swaths=['IW']
  TIFFTAG_DATETIME=2025:04:07 17:08:48
  TIFFTAG_IMAGEDESCRIPTION=Sentinel-1A IW GRD HR L1
  TIFFTAG_SOFTWARE=Sentinel-1 IPF 003.91
  transmitter_receiver_polarisations=['VV', 'VH']
  units=m2 m-2
Image Structure Metadata:
  COMPRESSION=ZSTD
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  20.0000000,  55.0000000) ( 20d 0' 0.00"E, 55d 0' 0.00"N)
Lower Left  (  20.0000000,  60.0000000) ( 20d 0' 0.00"E, 60d 0' 0.00"N)
Upper Right (  25.0000000,  55.0000000) ( 25d 0' 0.00"E, 55d 0' 0.00"N)
Lower Right (  25.0000000,  60.0000000) ( 25d 0' 0.00"E, 60d 0' 0.00"N)
Center      (  22.5000000,  57.5000000) ( 22d30' 0.00"E, 57d30' 0.00"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  Description = beta nought

And here is the gdalinfo of the srtm used

➜  ~ gdalinfo /home/priit/SAR/srtm/srtm_41_01.tif                       
Driver: GTiff/GeoTIFF
Files: /home/priit/SAR/srtm/srtm_41_01.tif
Size is 6000, 6000
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (20.000000000000000,60.000000000000000)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  20.0000000,  60.0000000) ( 20d 0' 0.00"E, 60d 0' 0.00"N)
Lower Left  (  20.0000000,  55.0000000) ( 20d 0' 0.00"E, 55d 0' 0.00"N)
Upper Right (  25.0000000,  60.0000000) ( 25d 0' 0.00"E, 60d 0' 0.00"N)
Lower Right (  25.0000000,  55.0000000) ( 25d 0' 0.00"E, 55d 0' 0.00"N)
Center      (  22.5000000,  57.5000000) ( 22d30' 0.00"E, 57d30' 0.00"N)
Band 1 Block=6000x1 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

As you can see, it seems to be flipped.

[priit111]

Additionally, if I print out dem_ecef.coords in apps.py, I get something like this:

Coordinates:
  * axis         (axis) int64 24B 0 1 2
  * y            (y) float64 48kB 55.0 55.0 55.0 55.0 ... 60.0 60.0 60.0 60.0
  * x            (x) float64 48kB 20.0 20.0 20.0 20.0 ... 25.0 25.0 25.0 25.0
    spatial_ref  int64 8B dask.array<chunksize=(), meta=np.ndarray>

Process finished with exit code 0```

[alexamici]

@priit111 oh, I see the point now. Thanks!

I've become used to open the DEM and the image with Xarray where the best practice is to have all increasing coordinates.

And I also recall now adding these two lines just for conveninece:

https://github.com/bopen/sarsen/blob/main/sarsen/scene.py#L21

I wonder if I can remove them or I need to add an option 😢