Alignment issues when source cog is in non-native projection

[staffordsmith83]

We are having problems where COGs are not positioned correctly on the cesium globe when they are reprojected by TIFFImageryProvider:

This COG is in EPSG:32756, available here: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/56/J/NP/2023/4/S2A_56JNP_20230410_0_L2A/TCI.tif

[hongfaqiu]

I don't know what's going on either, but there may be some issues with the simple bbox's coordinate system conversion. Do you know how to improve it?

[staffordsmith83]

Hi @hongfaqiu no I dont know how to improve it - but if you point me to some lines in the codebase, I could take a look?

[hongfaqiu]

Thank you for your help! All the code related to reprojection is here: https://github.com/hongfaqiu/TIFFImageryProvider/blob/main/packages/TIFFImageryProvider/src/TIFFImageryProvider.ts#L200-L227
I have just simply obtained the bounding box of the TIFF, converted it to the bounding box in the 4326 coordinate system, and passed it into the rectangle. May be I should rewrite the entire TilingScheme?

[staffordsmith83]

Im not sure, Ill have a look at this tomorrow, perhaps your method is not sufficient. Thanks for your work so far! We are trying to use your code in our open source tool TerriaJS

[hongfaqiu]

Thank you for your recognition! I really love TerriaJS, and have learned a lot from it.

[staffordsmith83]

I have confirmed that reprojecting a cog to EPSG:4326 fixes the alignment problem:

Datasets
Original cog that appears misaligned: https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/56/J/NP/2023/4/S2A_56JNP_20230410_0_L2A/TCI.tif
Reprojected with gdal and uploaded to S3: https://stafs-cogs.s3.ap-southeast-2.amazonaws.com/Whian_Whian_TCI_4326_cog.tif

[hongfaqiu]

Cool! I have found a repo geotiff-geokeys-to-proj4 that can perform coordinate system transformation on TIFF files. It seems to perform the transformation on every pixel rather than just the bounding box. I'm not sure if this will be effective.

[staffordsmith83]

Yes I dont think reprojecting the bounding box is sufficient. This is how OpenLayers handles it: https://openlayers.org/doc/tutorials/raster-reprojection.html

[hongfaqiu]

Great news! I solved this problem by creating a custom TilingScheme, you can pass projFunc to handle it.

import proj4 from 'proj4-fully-loaded'; 

TIFFImageryProvider.fromUrl(YOUR_TIFF_URL, {
  projFunc: (code) => {
    if (![4326].includes(code)) {
      {
        try {
          let prj = proj4(`EPSG:${code}`, "EPSG:4326")
          let unprj = proj4("EPSG:4326", `EPSG:${code}`)
          if (prj && unprj) return {
            project: prj.forward,
            unproject: unprj.forward
          }
        } catch (e) {
          console.error(e);
        }
      }
    }
  },
});

But the demo performance of nextjs was abnormal. I created a vite example and it worked very well. I will solve this problem in the afternoon.

[staffordsmith83]

@hongfaqiu this is great news, thanks for your efforts! I will check it out now

[staffordsmith83]

Dear @hongfaqiu , I have incorporated your changes into my terriajs test branch, and it resolved the alignment issue, thanks so much for your work! However, it does present a different problem: the COGs are sometimes displayed in the wrong location when zoomed out:

[staffordsmith83]

Dear @hongfaqiu , our cesium expert suggests that you can fix this by changing how this._rectangle is computed, to account for the northwest and southwest corners having different longitudes - I think this code should do it:

 const swMeters = new Cartesian3();
    options.rectangleSouthwestInMeters.clone(swMeters);
    const neMeters = new Cartesian3();
    options.rectangleNortheastInMeters.clone(neMeters);
    const seMeters = new Cartesian3(neMeters.x, swMeters.y);
    const nwMeters = new Cartesian3(swMeters.x, neMeters.y);

    const southwest = this.projection.unproject(swMeters);
    const southeast = this.projection.unproject(seMeters);
    const northwest = this.projection.unproject(nwMeters);
    const northeast = this.projection.unproject(neMeters);

 this._rectangle = new Rectangle(
      Math.min(
        southwest.longitude,
        southeast.longitude,
        northwest.longitude,
        northeast.longitude
      ),
      Math.min(
        southwest.latitude,
        southeast.latitude,
        northwest.latitude,
        northeast.latitude
      ),
      Math.max(
        southwest.longitude,
        southeast.longitude,
        northwest.longitude,
        northeast.longitude
      ),
      Math.max(
        southwest.latitude,
        southeast.latitude,
        northwest.latitude,
        northeast.latitude
      )
    );

[hongfaqiu]

Dear @staffordsmith83, thank you and your team for the work you have done. It's working great now!
There is still a lot of room for improvement in this library, and I look forward to working with you to build it together.
--edited--
Correction: There are still some issues, and trying to fix them.

[hongfaqiu]

When I project the non-3857 coordinate system into the 4326 coordinate system, I noticed that the four vertices constructed using the bbox attribute of the geotiff do not form a regular rectangle. This also causes the failure of constructing a Cesium.Rectangle.
The TilingScheme only performs projection transformation and mapping on regular rectangular tile boundaries, and it cannot handle tile seams.
I am considering implementing the coordinate transformation of projections in a different way in the frontend. It seems that the TilingScheme approach is incorrect.

[hongfaqiu]

I have developed a front-end raster reprojection method, which produces the following results as shown in the gif. It addresses some issues but there are still some remaining problems to be resolved.

[staffordsmith83]

Thanks for your work on this @hongfaqiu , we are following with great interest

[staffordsmith83]

Dear @hongfaqiu how did you go resolving these issues? Any luck?

[hongfaqiu]

Dear @staffordsmith83, I was busy with other things recently, but today I had some time to update the code. However, the results are not satisfactory, and there are areas for improvement.

Originally, my idea was to make the implementation of the tileXYToRectangle function consistent with the GeographicTilingScheme, where I would request additional tile data to ensure complete coverage after reprojection. However, this approach encountered some issues, and I had to switch to a solution similar to the WebMercatorTilingScheme.

In any case, the reprojection function is currently executing sequentially on the CPU, which makes it inefficient and blocks the rendering process. That's why I have marked the reprojection method as experimental.

[staffordsmith83]

Dear @hongfaqiu thanks somuch for your efforts! Yes its a difficult issue, the library that we are using for 2d mode, georaster-layer-for leaflet also has issues reprojecting. Also busy with other things, but can hopefully come back to this soon. Thanks again!

[hongfaqiu]

@staffordsmith83 Thank you very much for your ongoing support. I will also continue to follow your progress in terriaJS and look forward to finding a more effective general solution. It would be even better if this could help advance Cesium's capabilities in reprojection.

[staffordsmith83]

Dear @hongfaqiu did you make any progress with the reprojection? I notice in your changelog you have implemented some web based reprojection that is slow but works? Perhaps GeoWarp is an option? https://github.com/DanielJDufour/geowarp

[hongfaqiu]

Dear @hongfaqiu did you make any progress with the reprojection? I notice in your changelog you have implemented some web based reprojection that is slow but works? Perhaps GeoWarp is an option? https://github.com/DanielJDufour/geowarp

Excellent repo! I will try it next week. Thanks your sharing.

[hongfaqiu]

I have looked at DanielJDufour's code, and its reprojection is also based on CPU computation. Furthermore, it performs poorly on non-rectangular boundaries as it clips a portion of the content.
If there is no solution that offers performance improvements of 1-2 orders of magnitude, the feasibility of browser-based reprojection remains low.
I still have other work tasks this week, so I hope to have some free time next week to investigate this issue.