Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

reprojection function creates incorrect reprojection #28

Open
cmosig opened this issue Oct 26, 2024 · 4 comments
Open

reprojection function creates incorrect reprojection #28

cmosig opened this issue Oct 26, 2024 · 4 comments

Comments

@cmosig
Copy link

cmosig commented Oct 26, 2024 •
edited
Loading

steps to reproduce: tcd-reproject 20211001_FVA_Walddrohnen_Totholz_4_ortho.tif

original tif
image

reprojected tif overlayed:
image

src CRS: was 4326 here. can send file upon request
@jveitchmichaelis
Copy link
Collaborator

Do you know what the output reprojection is? (Ie what does GDAL think) could be we're defaulting to something and not respecting the projection of the source image.

@cmosig
Copy link
Author

cmosig commented Jan 2, 2025

gdalinfo of the input image:

Driver: GTiff/GeoTIFF
Files: 20211001_FVA_Walddrohnen_Totholz_4_ortho.tif
Size is 29960, 20427
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    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]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (8.049522352437883,47.968616561879600)
Pixel Size = (0.000000334476611,-0.000000334476611)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   8.0495224,  47.9686166) (  8d 2'58.28"E, 47d58' 7.02"N)
Lower Left  (   8.0495224,  47.9617842) (  8d 2'58.28"E, 47d57'42.42"N)
Upper Right (   8.0595433,  47.9686166) (  8d 3'34.36"E, 47d58' 7.02"N)
Lower Right (   8.0595433,  47.9617842) (  8d 3'34.36"E, 47d57'42.42"N)
Center      (   8.0545328,  47.9652004) (  8d 3'16.32"E, 47d57'54.72"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
  Unit Type: metre
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
  Unit Type: metre
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
  Unit Type: metre
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Unit Type: metre

gdalinfo output of the output image:

Driver: GTiff/GeoTIFF
Files: 20211001_FVA_Walddrohnen_Totholz_4_ortho_proj_10.tif
Size is 11155, 11325
Coordinate System is:
PROJCRS["WGS 84 / World Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["World Mercator",
        METHOD["Mercator (variant A)",
            ID["EPSG",9804]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Very small scale conformal mapping."],
        AREA["World between 80°S and 84°N."],
        BBOX[-80,-180,84,180]],
    ID["EPSG",3395]]
Data axis to CRS axis mapping: 1,2
Origin = (896068.729402458644472,6069881.303413412533700)
Pixel Size = (0.100000000000000,-0.100000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
  JPEGTABLESMODE=1
  JPEG_QUALITY=75
Corner Coordinates:
Upper Left  (  896068.729, 6069881.303) (  8d 2'58.28"E, 47d58' 7.02"N)
Lower Left  (  896068.729, 6068748.803) (  8d 2'58.28"E, 47d57'42.42"N)
Upper Right (  897184.229, 6069881.303) (  8d 3'34.36"E, 47d58' 7.02"N)
Lower Right (  897184.229, 6068748.803) (  8d 3'34.36"E, 47d57'42.42"N)
Center      (  896626.479, 6069315.053) (  8d 3'16.32"E, 47d57'54.72"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  NoData Value=0

@jveitchmichaelis
Copy link
Collaborator

jveitchmichaelis commented Jan 2, 2025 •
edited
Loading

OK I guess 3395 is a default setting (in metric units). The idea would be that we transform to a fixed 10cm resolution, predict, and then transform the predictions back into the original projection. In theory you'd never really need the resampled image, it's a convenience for the model. You could also see if QGIS will reproject the 10cm image into your original coordinate space?

Are the predictions also incorrect? I think it's more important that those are aligned with the input - though typically with a different GSD.

@cmosig
Copy link
Author

cmosig commented Jan 2, 2025 •
edited
Loading

Until now, for inference I first did a reprojection with gdalwarp into EPSG:8857 with 0.1m resolution and then applied tcd-predict as I didn't trust the tcd-reproject :)

When I run the prediction on the raw image, I get a different error (see output below). Should tcd-predict automatically do a reprojection internally?

(tcd) > $ tcd-predict "restor/tcd-segformer-mit-b5" 20211001_FVA_Walddrohnen_Totholz_4_ortho.tif prediction                   
INFO:tcd_pipeline.models.model:Device: cuda
INFO:tcd_pipeline.pipeline:Saving results to prediction
{'data': {'root': '???', 'train': 'train.json', 'validation': 'val.json', 'test': 'test.json', 'images': 'images', 'output': 'prediction', 'name': 'tcd', 'gsd': 0.1, 'num_workers': 16, 'scale_factor': 1, 'tile_size': 1024, 'tile_overlap': 512, 'classes': ['canopy', 'tree'], 'grayscale': False}, 'model': {'device': 'cuda', 'task': 'semantic_segmentation', 'num_workers': 16, 'tta': True, 'batch_size': 1, 'checkpoint': None, 'name': 'segformer', 'backbone': 'nvidia/mit-b5', 'revision': 'main', 'weights': 'restor/tcd-segformer-mit-b5', 'in_channels': 3, 'num_classes': 2, 'learning_rate': 0.0001, 'loss': 'focal', 'datamodule': None, 'augment': True, 'pretrained': 'imagenet', 'learning_rate_schedule_patience': 5, 'resume': False, 'trainer': {'gpus': 1, 'min_epochs': 1, 'max_epochs': 75, 'auto_lr_find': False, 'auto_scale_batch_size': False, 'precision': 16, 'early_stopping_patience': 10, 'max_time': '00:48:00:00', 'debug_run': False}, 'wandb': {'project_name': '???'}}, 'preprocess': {'bounds_check': True, 'remove_border': {'radius': 10}, 'min_size': {'tree': 30, 'canopy': 30}, 'max_size': 100, 'merge_touching': {'threshold': 2}, 'merge_enclosed': True}, 'postprocess': {'iou_threshold': 0.7, 'iou_tiles': 0, 'min-size': {'tree': 30, 'canopy': 30}, 'max-size': 100, 'cache_folder': './temp', 'stateful': True, 'cleanup': True, 'dissolve': True, 'debug_images': False, 'cache_format': 'geotiff', 'use_nms': True, 'confidence_threshold': 0.2, 'segmentation_merge_pad': 32}, 'evaluate': {}, 'input': '???', 'output': '???'}
INFO:tcd_pipeline.data.imagedataset:Geographic information present, loading as a geo dataset
INFO:tcd_pipeline.data.tiling:Source resolution is 0.0000
INFO:tcd_pipeline.data.imagedataset:Dataset has 1 tiles of size 1024x1024 px.
INFO:tcd_pipeline.cache.semantic:Caching to 9 tiles, approximately 0.90GB needed for temporary storage during inference
  0%|                                                                                                    | 0/1 [00:00<?, ?it/s]INFO:tcd_pipeline.cache.semantic:Caching to 9 tiles, approximately 0.90GB needed for temporary storage during inference
  0%|                                                                                                    | 0/1 [00:00<?, ?it/s]
Traceback (most recent call last):
  File "/net/home/cmosig/miniconda3/envs/tcd/bin/tcd-predict", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/net/home/cmosig/gitRepos/tcd/src/tcd_pipeline/scripts/predict.py", line 93, in main
    res = pipeline.predict(input, warm_start=args.resume)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/net/home/cmosig/gitRepos/tcd/src/tcd_pipeline/pipeline.py", line 204, in predict
    return self.model.predict_tiled(image, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/net/home/cmosig/gitRepos/tcd/src/tcd_pipeline/models/model.py", line 157, in predict_tiled
    for _, batch in progress_bar:
  File "/net/home/cmosig/miniconda3/envs/tcd/lib/python3.12/site-packages/tqdm/std.py", line 1181, in __iter__
    for obj in iterable:
  File "/net/home/cmosig/miniconda3/envs/tcd/lib/python3.12/site-packages/torch/utils/data/dataloader.py", line 631, in __next__
    data = self._next_data()
           ^^^^^^^^^^^^^^^^^
  File "/net/home/cmosig/miniconda3/envs/tcd/lib/python3.12/site-packages/torch/utils/data/dataloader.py", line 675, in _next_data
    data = self._dataset_fetcher.fetch(index)  # may raise StopIteration
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/net/home/cmosig/miniconda3/envs/tcd/lib/python3.12/site-packages/torch/utils/data/_utils/fetch.py", line 51, in fetch
    data = [self.dataset[idx] for idx in possibly_batched_index]
            ~~~~~~~~~~~~^^^^^
  File "/net/home/cmosig/gitRepos/tcd/src/tcd_pipeline/data/imagedataset.py", line 44, in __getitem__
    return self.tiles[idx]
           ~~~~~~~~~~^^^^^
  File "/net/home/cmosig/gitRepos/tcd/src/tcd_pipeline/data/tiling.py", line 416, in __getitem__
    data, window = self._get_data(idx)
                   ^^^^^^^^^^^^^^^^^^^
  File "/net/home/cmosig/gitRepos/tcd/src/tcd_pipeline/data/tiling.py", line 393, in _get_data
    self.dataset.read(
  File "rasterio/_io.pyx", line 627, in rasterio._io.DatasetReaderBase.read
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 333. PiB for an array with shape (4, 306149957, 306149957) and data type uint8

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@jveitchmichaelis @cmosig