Skip to content

Commit

Permalink
USGSDEM: remove write support
Browse files Browse the repository at this point in the history
  • Loading branch information
rouault committed Jan 28, 2025
1 parent 76a65b8 commit 35c05e9
Show file tree
Hide file tree
Showing 7 changed files with 5 additions and 1,834 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ Supported Formats: (ro:read-only, rw:read-write, +:update, v:virtual-I/O s:subda
NOAA_B -raster- (rov): NOAA GEOCON/NADCON5 .b format (*.b)
NSIDCbin -raster- (rov): NSIDC Sea Ice Concentrations binary (.bin) (*.bin)
RIK -raster- (rov): Swedish Grid RIK (.rik) (*.rik)
USGSDEM -raster- (rwv): USGS Optional ASCII DEM (and CDED) (*.dem)
USGSDEM -raster- (rov): USGS Optional ASCII DEM (and CDED) (*.dem)
GXF -raster- (rov): GeoSoft Grid Exchange Format (*.gxf)
BAG -raster,multidimensional raster,vector- (rw+v): Bathymetry Attributed Grid (*.bag)
S102 -raster,multidimensional raster- (rovs): S-102 Bathymetric Surface Product (*.h5)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ Supported Formats: (ro:read-only, rw:read-write, +:update, v:virtual-I/O s:subda
BYN -raster- (rov): Natural Resources Canada's Geoid (*.byn, *.err)
NOAA_B -raster- (rov): NOAA GEOCON/NADCON5 .b format (*.b)
RIK -raster- (rov): Swedish Grid RIK (.rik) (*.rik)
USGSDEM -raster- (rwv): USGS Optional ASCII DEM (and CDED) (*.dem)
USGSDEM -raster- (rov): USGS Optional ASCII DEM (and CDED) (*.dem)
GXF -raster- (rov): GeoSoft Grid Exchange Format (*.gxf)
KEA -raster- (rw+v): KEA Image Format (.kea) (*.kea)
BAG -raster,multidimensional raster,vector- (rw+v): Bathymetry Attributed Grid (*.bag)
Expand Down
144 changes: 1 addition & 143 deletions autotest/gdrivers/usgsdem.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
import gdaltest
import pytest

from osgeo import gdal, osr
from osgeo import osr

pytestmark = pytest.mark.require_driver("USGSDEM")

Expand Down Expand Up @@ -75,148 +75,6 @@ def test_usgsdem_3():
)


###############################################################################
# Test CreateCopy()


def test_usgsdem_4():

tst = gdaltest.GDALTest(
"USGSDEM",
"usgsdem/39079G6_truncated.dem",
1,
61424,
options=["RESAMPLE=Nearest"],
)
tst.testCreateCopy(check_gt=1, check_srs=1, vsimem=1)


###############################################################################
# Test CreateCopy() without any creation options


@pytest.mark.require_driver("DTED")
def test_usgsdem_5():

ds = gdal.Open("data/n43.dt0")
ds2 = gdal.GetDriverByName("USGSDEM").CreateCopy(
"tmp/n43.dem", ds, options=["RESAMPLE=Nearest"]
)

if ds.GetRasterBand(1).Checksum() != ds2.GetRasterBand(1).Checksum():
print(ds2.GetRasterBand(1).Checksum())
print(ds.GetRasterBand(1).Checksum())
ds2 = None
print(open("tmp/n43.dem", "rb").read())
pytest.fail("Bad checksum.")

gt1 = ds.GetGeoTransform()
gt2 = ds2.GetGeoTransform()
for i in range(6):
if gt1[i] != pytest.approx(gt2[i], abs=1e-5):
print("")
print("old = ", gt1)
print("new = ", gt2)
pytest.fail("Geotransform differs.")

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
assert ds2.GetProjectionRef() == srs.ExportToWkt(), "Bad SRS."

ds2 = None


###############################################################################
# Test CreateCopy() without a few creation options. Then create a new copy with TEMPLATE
# creation option and check that both files are binary identical.


@pytest.mark.require_driver("DTED")
def test_usgsdem_6():

ds = gdal.Open("data/n43.dt0")
ds2 = gdal.GetDriverByName("USGSDEM").CreateCopy(
"tmp/file_1.dem",
ds,
options=[
"PRODUCER=GDAL",
"OriginCode=GDAL",
"ProcessCode=A",
"RESAMPLE=Nearest",
],
)

ds3 = gdal.GetDriverByName("USGSDEM").CreateCopy(
"tmp/file_2.dem", ds2, options=["TEMPLATE=tmp/file_1.dem", "RESAMPLE=Nearest"]
)

del ds2
del ds3

f1 = open("tmp/file_1.dem", "rb")
f2 = open("tmp/file_2.dem", "rb")

# Skip the 40 first bytes because the dataset name will differ
f1.seek(40, 0)
f2.seek(40, 0)

data1 = f1.read()
data2 = f2.read()

assert data1 == data2

f1.close()
f2.close()


###############################################################################
# Test CreateCopy() with CDED50K profile


@pytest.mark.require_driver("DTED")
def test_usgsdem_7():

ds = gdal.Open("data/n43.dt0")

# To avoid warning about 'Unable to find NTS mapsheet lookup file: NTS-50kindex.csv'
with gdal.quiet_errors():
ds2 = gdal.GetDriverByName("USGSDEM").CreateCopy(
"tmp/000a00DEMz",
ds,
options=[
"PRODUCT=CDED50K",
"TOPLEFT=80w,44n",
"RESAMPLE=Nearest",
"ZRESOLUTION=1.1",
"INTERNALNAME=GDAL",
],
)

assert ds2.RasterXSize == 1201 and ds2.RasterYSize == 1201, "Bad image dimensions."

expected_gt = (
-80.000104166666674,
0.000208333333333,
0,
44.000104166666667,
0,
-0.000208333333333,
)
got_gt = ds2.GetGeoTransform()
for i in range(6):
if expected_gt[i] != pytest.approx(got_gt[i], abs=1e-5):
print("")
print("expected = ", expected_gt)
print("got = ", got_gt)
pytest.fail("Geotransform differs.")

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("NAD83")
assert ds2.GetProjectionRef() == srs.ExportToWkt(), "Bad SRS."

ds2 = None


###############################################################################
# Test truncated version of http://download.osgeo.org/gdal/data/usgsdem/various.zip/39109h1.dem
# Undocumented format
Expand Down
117 changes: 0 additions & 117 deletions doc/source/drivers/raster/usgsdem.rst
Original file line number Diff line number Diff line change
Expand Up @@ -33,127 +33,10 @@ VTP.
Driver capabilities
-------------------

.. supports_createcopy::

.. supports_georeferencing::

.. supports_virtualio::

Creation Issues
---------------

GDAL supports export of geographic (and UTM) USGS DEM and CDED data
files, including the ability to generate CDED 2.0 50K products to
Canadian federal government specifications.

Input data must already be sampled in a geographic or UTM coordinate
system. By default the entire area of the input file will be output, but
for CDED50K products the output file will be sampled at the production
specified resolution and on product tile boundaries.

If the input file has appropriate coordinate system information set,
export to specific product formats can take input in different
coordinate systems (i.e. from Albers projection to NAD83 geographic for
CDED50K production).

|about-creation-options|
The following creation options are supported:

- .. co:: PRODUCT
:choices: DEFAULT, CDED50K

When CDED50K is specified, the output
file will be forced to adhere to CDED 50K product specifications. The
output will always be 1201x1201 and generally a 15 minute by 15
minute tile (though wider in longitude in far north areas).

- .. co:: TOPLEFT
:choices: <long\,lat>

For CDED50K products, this is used to specify
the top left corner of the tile to be generated. It should be on a 15
minute boundary and can be given in decimal degrees or degrees and
minutes (eg. TOPLEFT=117d15w,52d30n).

- .. co:: RESAMPLE
:choices: Nearest, Bilinear, Cubic, CubicSpline
:default: Bilinear

Set the resampling
kernel used for resampling the data to the target grid. Only has an
effect when particular products like CDED50K are being produced.

- .. co:: DEMLevelCode
:choices: 1, 2, 3
:default: 1

DEM Level (1, 2 or 3 if set).

- .. co:: DataSpecVersion
:choices: <integer>

Data and Specification version/revision
(eg. 1020)

- .. co:: PRODUCER
:choices: <text>

Up to 60 characters to be put into the producer
field of the generated file.

- .. co:: OriginCode
:choices: <text>

Up to 4 characters to be put into the origin
code field of the generated file (YT for Yukon).

- .. co:: ProcessCode
:choices: <character>

One character to be put into the process code
field of the generated file (8=ANUDEM, 9=FME, A=TopoGrid).

- .. co:: TEMPLATE
:choices: <filename>

For any output file, a template file can be
specified. A number of fields (including the Data Producer) will be
copied from the template file if provided, and are otherwise left
blank.

- .. co:: ZRESOLUTION
:default: 1.0

DEM's store elevation information as positive
integers, and these integers are scaled using the "z resolution." By
default, this resolution is written as 1.0. However, you may specify
a different resolution here, if you would like your integers to be
scaled into floating point numbers.

- .. co:: NTS
:choices: <name>

NTS Mapsheet name, used to derive TOPLEFT. Only has an
effect when particular products like CDED50K are being produced.

- .. co:: INTERNALNAME
:choices: <name>

Dataset name written into file header. Only
has an effect when particular products like CDED50K are being
produced.

Example: The following would generate a single CDED50K tile, extracting
from the larger DEM coverage yk_3arcsec for a tile with the top left
corner -117w,60n. The file yk_template.dem is used to set some product
fields including the Producer of Data, Process Code and Origin Code
fields.

::

gdal_translate -of USGSDEM -co PRODUCT=CDED50K -co TEMPLATE=yk_template.dem \
-co TOPLEFT=-117w,60n yk_3arcsec 031a01_e.dem

--------------

NOTE: Implemented as :source_file:`frmts/usgsdem/usgsdemdataset.cpp`.
Expand Down
2 changes: 1 addition & 1 deletion frmts/usgsdem/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
add_gdal_driver(TARGET gdal_USGSDEM SOURCES usgsdem_create.cpp usgsdemdataset.cpp PLUGIN_CAPABLE NO_DEPS)
add_gdal_driver(TARGET gdal_USGSDEM SOURCES usgsdemdataset.cpp PLUGIN_CAPABLE NO_DEPS)
gdal_standard_includes(gdal_USGSDEM)
target_include_directories(gdal_USGSDEM PRIVATE $<TARGET_PROPERTY:alg,SOURCE_DIR> ${GDAL_RASTER_FORMAT_SOURCE_DIR}/mem)
Loading

0 comments on commit 35c05e9

Please sign in to comment.