I can't figure out how to warp a large file (2GB+) with no specified to CRS, to a specified CRS - I keep getting a "Too big image size" error.
For context, I want to take large NITF files, warp them to a set coordinate system (EPSG 3857, 4326, etc.), create tiles via gdal_retile, and stick that in a GeoServer image pyramid to use for a wms in Leaflet. Everything in this workflow is great as long as the file is under a certain size (~770 MB or so). Unfortunately, my files can be over 2 GB. My example file is 791 MB, and is from U Arkansas CAST's Corona data set.
When I run gdalwarp on my large file, I get the following error:
$ gdalwarp -t_srs EPSG:4326 -of NITF -co "ICORDS=D" -rpc ds1042-2184df057.ntf ds1042-2184df057.w.ntf Creating output file that is 1635251P x 522168L. ERROR 1: Unable to create file ds1042-2184df057.w.ntf, Too big image size : 853875744168
Since the file is too big to gdalwarp in its native state, I need to cut it up into pieces, warp the subparts individually, and then use gdal_merge to put them back together (I think).
I've investigated gdal2tiles.py and gdal_translate, but gdal_retile seems like a good bet to have a minimal amount of files to preprocess and merge. My (possibly not the greatest) plan is to cut the too-big original images into 6 small parts, so that even if the size is slightly over 2GB, the image sub-parts will not be so big that gdalwarp can't handle them.
I get the image size from gdalinfo:
$ gdalinfo -nomd -noct -norat ds1042-2184df057.ntf Driver: NITF/National Imagery Transmission Format Files: ds1042-2184df057.ntf Size is 107957, 7590 Coordinate System is `' GCP Projection = GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9108"]], AUTHORITY["EPSG","4326"]] GCP[ 0]: Id=UpperLeft, Info= (0.5,0.5) -> (32.5269444444444,34.9538888888889,0) GCP[ 1]: Id=UpperRight, Info= (107956.5,0.5) -> (35.4597222222222,35.3908333333333,0) GCP[ 2]: Id=LowerRight, Info= (107956.5,7589.5) -> (35.5452777777778,35.2005555555556,0) GCP[ 3]: Id=LowerLeft, Info= (0.5,7589.5) -> (32.5286111111111,34.7463888888889,0) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 7590.0) Upper Right (107957.0, 0.0) Lower Right (107957.0, 7590.0) Center (53978.5, 3795.0) Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
Then, I divide the width by 3 and the height by 2 (rounding up if needed), so that I get 6 image sub-parts. 107957/3 rounds up to 35986, and 7590/2 = 3795. Then, I try to subdivide the image, which is met with total failure:
$ gdal_retile.py -r bilinear -levels 1 -ps 35986 3795 -targetDir subpart_test ds1042-2184df057.ntf ERROR 1: Unable to compute a transformation between pixel/line and georeferenced coordinates for TEMP. There is no affine transformation and no GCPs. Reprojection failed for subpart_test\1\ds1042-2184df057_1_1.tif, error 3
- I tried using gdal_translate on the original file to both turn it into a GeoTiff and manually add the gcps shown in gdalinfo, and then ran gdal_retile, but that had no effect - still got a gdal_retile error about affine transformation/ no GCPs .
- gdal2tiles.py doesn't look like I can specify the image size, and I'll just end up with a ton of small files to loop through.
- I read about gdal_translate in this issue (Splitting raster into smaller chunks using GDAL?), but had hoped gdal_retile would perform the same task on the command line, without writing an extra script to loop through multiple pixel windows.
- I tried using gdaltransform instead of gdalwarp to set the CRS, but afterward, gdal_retile still complains about a lack of affine transformation and GCPs on files that are the result of gdaltransform.
So, what is the best way to take a large file (2GB+) that lacks a specified Coordinate System, and warp it to a specified CRS? I'm new to GDAL, so am hopefully missing something obvious.