Subdividing large files due to gdalwarp - too big image size

Question

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

Other approaches:

• 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.