1. Topics
  2. Any idea why I'd be getting a missing region when reprojecting?

Any idea why I'd be getting a missing region when reprojecting?

hi@...
 

I'll try my best to describe the issue although it's really confusing me so I'm not sure how well I can explain it.

I'm trying to reproject some satellite imagery in the form of a raster that covers the entire world using the CRS of the imagery to EPSG:4326 using the reproject function. The raster before being reprojected uses a sinusoidal projection with meters as units and looks like this:


When I try to reproject it, I get the following:


As you can see, anything to the right of a certain horizontal value in the original raster is just missing entirely.

I really can't figure out why this is happening and I've tried changing dozens of things in my code to fix it with very little success. I've provided some of my code below.

Sorry for the code formatting, couldn't see any option for better formatting.

for label in layers.keys():
    with rasterio.open(os.path.join(OUTPUT_DIR, f'{label}_intermediate.tif')) as intermediate_file:
        # Compute transform needed for output
        final_transform, final_width, final_height = calculate_default_transform(
            intermediate_file.crs,
            FINAL_CRS,
            intermediate_file.width,
            intermediate_file.height,
            *intermediate_file.bounds,
            resolution = 0.05
        )
 
        # Extract parameters from intermediate file
        kwargs = intermediate_file.meta.copy()
 
        # Update using the final parameters
        kwargs.update({
            'crs': FINAL_CRS,
            'transform': final_transform,
            'width': final_width,
            'height': final_height
        })
 
        # Compute reprojection and write to file
        with rasterio.open(os.path.join(OUTPUT_DIR, f'{label}_final.tif'), 'w', **kwargs) as final_file:
            dst, tsf = reproject(
                source = rasterio.band(intermediate_file, 1),
                destination = rasterio.band(final_file, 1),
                src_transform = intermediate_file.transform,
                src_crs = intermediate_file.crs,
                dst_transform = final_transform,
                dst_crs = FINAL_CRS,
                dst_resolution = 0.05,
                resampling = Resampling.average,
                warp_mem_limit = 1024 * REPROJECTION_MEMORY_GB,
                kwargs = {'COMPRESSION': 'LZW', 'TILED': 'YES'}
            )

So, like I've said, I feel like I've tried everything and with very little success. Here's a summary.
 
  • I tried changing the ordering of width and height since rasterio's functions aren't consistent in the ordering but this didn't help.
  • I tried changing the resampling method... nothing.
  • I thought it might be a memory problem so increased/decreased the memory of the reproject function... nothing. Also downscaled the original raster to reduce memory but again nothing.
  • I tried processing a vertical slice of the data covering only the affected region but the output raster was still hidden. However, when I removed the resolution parameters so it would maintain the original resolution the region did manage to appear properly. But, when I then try to repeat this on the global raster I just get a blank output which seems to happen when I process files that are too large although I'm not sure why.
  • Continuing with the vertical slice... if I set the resolution parameters back to how they are in the code but then modify the output's transform I'm able to have the regions appear again. I modify the transform by rounding the upper left corner coordinates so that the bound is from -180 to 180 rather than what is originally set which has some floating point error and is more like -179.999995 to 180.000005. However, again, when I run this on the global raster I just get an empty raster. I don't think it's a memory issue since it's using the same resolution used to generate the global raster with the missing region.
  • I was able to slightly increase the amount of data that was visible but I can't recall what exactly I changed to do this. I know it involved increasing/decreasing one of the width variables or something along those lines which slightly extended the "visible" region. Not very much, but enough that there was maybe another degree that was visible. I tried further increasing/decreasing this variable but that didn't help.
I'm definitely feeling fairly lost. There's a good chance I'm simply doing something that I don't understand well as I'm new to this stuff. Let me know if there's any additional info I can provide. Any help would be greatly appreciated :)
hi@...
 

Just an update: I set the env variable CHECK_WITH_INVERT_PROJ to YES and was able to generate the output using the vertical slice but again, when I try to do this with the global raster I get a blank output. I feel like I may have solved the original problem I posted about but am now running into an issue with generating the global raster.
hi@...
 

Another update, for posterity. I tried doing the same reprojection using the gdal.Warp and it worked perfectly with far less effort around setting up transforms and whatnot. I filed an issue with rasterio: github.com/mapbox/rasterio/issues/2074
1 - 3 of 3