1

I downloaded some rasters from the USGS High Resolution National Hydrography Dataset Plus. In particular, I am interested in the flow accumulation images ("fac.tif" in the downloaded folders). I am uploading these images to Earth Engine, but curiously about 1/4 of the images give me an error warning:

Error: Projection for fac.tif could not be determined. Make sure both CRS and affine transform are present.

These problem TIFs look identical to the non-problem TIFs and I am not sure how to proceed. I tried importing them to ArcMap and then exporting them as TIFs, thinking maybe the projections would be corrected, but then the size of the exported image ballooned from 400MB to 5GB. However, the problem TIFs obviously have CRS info because ArcMap placed them correctly on a map.

Any ideas for how I might fix this or at least diagnose the problem? I was thinking maybe something with GDAL to set the CRS manually but I'm not very experienced with it.

P.S. Here is a link to a working TIF: https://www.sciencebase.gov/catalog/item/5d30c28be4b01d82ce84a9ec

And a non-working TIF: https://www.sciencebase.gov/catalog/item/5d30c28de4b01d82ce84a9fe

EDIT:

Here are the gdalinfo results for both rasters.

Problem Raster

Driver: GTiff/GeoTIFF
Files: fac.tif
       fac.tif.ovr
       fac.tif.aux.xml
Size is 36774, 42216
Coordinate System is:
PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",
    GEOGCS["GCS_North_American_1983",
        DATUM["D_North_American_1983",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Albers"],
    PARAMETER["false_easting",0.0],
    PARAMETER["false_northing",0.0],
    PARAMETER["central_meridian",-96.0],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_origin",23.0],
    UNIT["Meter",1.0],
    VERTCS["NAVD_1988",
        VDATUM["North_American_Vertical_Datum_1988"],
        PARAMETER["Vertical_Shift",0.0],
        PARAMETER["Direction",1.0],
        UNIT["Centimeter",0.01]]]
Origin = (-301195.000000000000000,1062005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -301195.000, 1062005.000)
Lower Left  ( -301195.000,  639845.000)
Upper Right (   66545.000, 1062005.000)
Lower Right (   66545.000,  639845.000)
Center      ( -117325.000,  850925.000)
Band 1 Block=128x128 Type=Int32, ColorInterp=Gray
  Min=0.000 Max=294046826.000
  Minimum=0.000, Maximum=294046826.000, Mean=41553.746, StdDev=2825681.435
  NoData Value=-2147483648
  Overviews: 18387x21108, 9194x10554, 4597x5277, 2299x2639, 1150x1320, 575x660, 288x330, 144x165
  Metadata:
    STATISTICS_COVARIANCES=7984475570911.16
    STATISTICS_MAXIMUM=294046826
    STATISTICS_MEAN=41553.745648755
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=2825681.434789

And for the working TIF

Driver: GTiff/GeoTIFF
Files: fac.tif
       fac.tif.ovr
       fac.tif.aux.xml
Size is 43748, 24574
Coordinate System is:
PROJCS["NAD_1983_Albers",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",23],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-53125.000000000000000,1487025.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -53125.000, 1487025.000) ( 96d35'52.93"W, 36d25'46.73"N)
Lower Left  (  -53125.000, 1241285.000) ( 96d34'52.04"W, 34d14' 2.11"N)
Upper Right (  384355.000, 1487025.000) ( 91d40'34.26"W, 36d21'11.04"N)
Lower Right (  384355.000, 1241285.000) ( 91d47'53.97"W, 34d 9'33.66"N)
Center      (  165615.000, 1364155.000) ( 94d 9'45.35"W, 35d19'10.20"N)
Band 1 Block=128x128 Type=Int32, ColorInterp=Gray
  Min=0.000 Max=401738279.000
  Minimum=0.000, Maximum=401738279.000, Mean=35430.142, StdDev=2923924.998
  NoData Value=-2147483647
  Overviews: 21874x12287, 10937x6144, 5469x3072, 2735x1536, 1368x768, 684x384, 342x192
  Metadata:
    RepresentationType=ATHEMATIC
    STATISTICS_COVARIANCES=8549337391801.696
    STATISTICS_MAXIMUM=401738279
    STATISTICS_MEAN=35430.141633088
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=2923924.9976362
2
  • 1
    Use gdalinfo for working and not-working tiffs and add the results. Ballooning size means that the original tiff was compressed but you created an uncompressed version.
    – user30184
    May 6, 2020 at 20:44
  • @user30184 Okay I added the results. Yeah it looks like they're both LZW compressed. I suppose I could re-export from ArcMap with LZW compression and see if that works, but that would be cumbersome with many images. If I could try something with GDAL, that would be ideal.
    – la_leche
    May 6, 2020 at 22:03

1 Answer 1

2

gdalwarp will do the job of converting the CRS as well as compressing the output.

Something like:

gdalwarp -t_srs EPSG:5072 -of GTiff -co "COMPRESS=DEFLATE" -co "PREDICTOR=2" infile.tif outfile.tif

This will use "deflate" compression with horizontal prediction, which I've found to be better than LZW. Depending on your gdal installation and where you plan to use the resulting file you could try ZSTD for the compression, but in my limited experience it's not as widely supported.

2
  • 1
    EPSG:9001 in the place where it is in the working WKT means the unit, metre epsg-registry.org/export.htm?gml=urn:ogc:def:uom:EPSG::9001 so that is not the EPSG code of that projection. Gdalwarp command is otherwise good. I would add also -co tiled=yes.
    – user30184
    May 7, 2020 at 7:04
  • 1
    @user30184 - you're right, of course. I was too quick to grab the first EPSG number I saw. I've edited the response to use EPSG:5072 which is an Albers projection for CONUS. If the OP wants to dig deeper, it's possible to provide a WKT in a separate file. This WKT could exactly match the projection in the file that works.
    – Llaves
    May 7, 2020 at 14:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.