pySW4.utils.geo module¶

Python module for handling DEM/GDEM/DTM and other raster data readable with gdal.

author:	Shahar Shani-Kadmiel s.shanikadmiel@tudelft.nl
copyright:	Shahar Shani-Kadmiel
license:	This code is distributed under the terms of the GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)

class GeoTIFF[source]¶

Bases: object

Class for handling GeoTIFF files.

In the GeoTIFF file format, it is assumed that the origin of the grid is at the top-left (north-west) corner. Hence, dx is positive in the x or lon direction and dy is negative in the y or lat direction.

Note

This class should be populated by read_GeoTIFF() or get_dem().

_read(filename, rasterBand, substring='*_dem.tif', verbose=False)[source]¶: Private method for reading a GeoTIFF file.

copy()[source]¶: Return a deepcopy of the GeoTIFF object.

dtype¶

e¶

elev¶

elevation_profile(xs, ys)[source]¶: Get the elevation values ar points in xs and ys.

extent¶

get_intensity(azimuth=315.0, altitude=45.0, scale=None, smooth=None, normalize=False)[source]¶

This is a method to create an intensity array that can be used as illumination to create a shaded relief map or draping data over. It is assumed that the grid origin is at the upper-left corner. If that is not the case, add 90 to azimuth.

Parameters:

relief : a 2d ndarray

Topography or other data to calculate intensity from.

azimuth : float

Direction of light source, degrees from north.

altitude : float

Height of light source, degrees above the horizon.

scale : float

Scaling value of the data.

smooth : float

Number of cells to average before intensity calculation.

normalize : bool

By default the intensity is clipped to the [0,1] range. If set to True, intensity is normalized to [0,1].

Returns:

intensity : ndarray

a 2d array with illumination data clipped in the [0,1] range. Optionally, can be normalized (instead of clipped) to [0,1]. Same size as relief.

keep(w, e, s, n)[source]¶: Deprecated! Will dissapear in future vertions. Please use set_new_extent() method.

nx¶

ny¶

reproject(epsg=None, proj4=None, match=None, spacing=None, interpolation='nearest', error_threshold=0.125, target_filename=None, verbose=False)[source]¶

Reproject the data from the current projection to the specified target projection epsg or proj4 or to match an existing GeoTIFF file. Optionally, the grid spacing can be changed as well.

To keep the same projection and only resample the data leave epsg, proj4, and match to be None and specify spacing.

Warning

Watch Out! This operation is performed in place on the actual data. The raw data will no longer be accessible afterwards. To keep the original data, use the copy() method to create a copy of the current object.

Parameters:

epsg : int

The target projection EPSG code. See the Geodetic Parameter Dataset Registry for more information.

proj4 : str

The target Proj4 string. If the EPSG code is unknown or a custom projection is required, a Proj4 string can be passed. See the Proj4 documentation for a list of general Proj4 parameters.

match : str or GeoTIFF instance

Path (relative or absolute) to an existing GeoTIFF file or GeoTIFF instance (already in memory) to match size and projection of. Current data is resampled to match the shape and number of pixels of the existing GeoTIFF file or object. It is assumed that both GeoTIFF objects cover the same extent.

spacing : int or float

Target spacing in the units of ther target projection.

interpolation : {‘nearest’, ‘bilinear’, ‘cubic’, ‘lanczos’}

Resampling algorithm:

‘nearest’ : Nearest neighbour

‘bilinear’ : Bilinear (2x2 kernel)

‘cubic’ : Cubic Convolution Approximation (4x4 kernel)

‘lanczos’ : Lanczos windowed sinc interpolation (6x6 kernel)

error_threshold : float

Error threshold for transformation approximation in pixel units. (default is 0.125 for compatibility with gdalwarp)

target_filename : str

If a target filename is given then the reprojected data is saved as target_filename and read into memory replacing the current data. This is faster than reprojecting and then saving. Otherwise the write_GeoTIFF() method can be used to save the data at a later point if further manipulations are needed.

Notes

Examples of some EPSG codes and their equivalent Proj4 strings are:

4326   -> '+proj=longlat +datum=WGS84 +no_defs'

32636  -> '+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs'

102009 -> '+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40
           +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m
           +no_defs'

and so on … See the Geodetic Parameter Dataset Registry for more information.

resample(by=None, to=None, order=0)[source]¶

Method to resample the data either by a factor or to the specified spacing.

This method uses zoom().

Warning

Watch Out! This operation is performed in place on the actual data. The raw data will no longer be accessible afterwards. To keep the original data, use the copy() method to create a copy of the current object.

Parameters:

by : float

If by=1, nothing happens.

If by<1, the grid is sub-sampled by the factor.

if by>1, the grid is super-sampled byt the factor.

to : float

The specified spacing to which the grid is resampled to.

order : int

The order of the spline interpolation, default is 0. The order has to be in the range 0-5.

0 - nearest (fastest)

1 - bi-linear

3 - cubic (slower)

…

Note that this may result in slightly different spacing than desired and more importantly, may cause a discrepency between x and y spacing.

s¶

set_new_extent(w, e, s, n, fill_value=None, mask=False)[source]¶

Change the extent of a GeoTIFF file, trimming the boundaries or padding them as needed.

Warning

Watch Out! This operation is performed in place on the actual data. The raw data will no longer be accessible afterwards. To keep the original data, use the copy() method to create a copy of the current object.

Parameters:

w, e, s, n : float

The west-, east-, south-, and north-most coordinate to keep (may also be the xmin, xmax, ymin, ymax coordinate).

fill_value : float or None

Value to pad the data with incase extent is expanded beyond the data. By default this is set to the nodata value.

mask : bool

If set to True, data is masked using the fill_value.

set_nodata(value)[source]¶: Replace the current nodata value with the new value.

write_GeoTIFF(filename, nodata=None)[source]¶

Write a GeoTIFF file.

Parameters:

filename : str

Path (relative or absolute) to the saves file.

nodata : int or float or NaN

Value to register as the no-data value in the GeoTIFF file header.

write_topo_file(filename)[source]¶: Write an ASCII Lon, Lat, z file for SW4 topography.

x¶: Returns x coordinates vector of the data. This usually corresponds to Longitude.

xy2d¶: Returns 2 2d ndarray s with x and y coordinates of the data. This is the result of meshgrid() of the x and y vectors of the data.

y¶: Returns y coordinates vector of the data. This usually corresponds to Latitude.

_get_tiles(path, lonmin, lonmax, latmin, latmax, tilename='ASTGTM2_{}{:02d}{}{:03d}', verbose=False)[source]¶

This is a helper function for get_dem() function.

It generates a list of tile filenames that are needed for the stitching process.

Note

This function should not be used directly by the user. See get_dem().

gdalopen(filename, substring='*_dem.tif')[source]¶

Wrapper function around gdal.Open().

gdal.Open() returns None if no file is read. This wrapper function will try to read the filename directly but if None is returned, several other options, including compressed file formats (‘zip’, ‘gzip’, ‘tar’) and url retrieval are attempted.

get_dem(path, lonmin, lonmax, latmin, latmax, tilename='ASTGTM2_{}{:02d}{}{:03d}', rasterBand=1)[source]¶

This function reads ASTER-GDEM GeoTIFF tiles into memory, stitches them together if more than one file is read, and cuts to the desired extent given by lonmin, lonmax, latmin, latmax in decimal degrees.

Parameters:

path : str

Path (relative or absolute) to where the ASTER-GDEM zipped tiles are stored.

lonmin, lonmax : float

The west- and east-most coordinate of the desired extent (may also be the xmin and xmax coordinate).

latmin, latmax : float

The south- and north-most coordinate of the desired extent (may also be the ymin and ymax coordinate).

tilename : str

String template of the tile names.

rasterBand : int

The band number to read from each tile (defaults to 1).

Returns:

GeoTIFF

A populated GeoTIFF instance with the stitched elevation data.

read_GeoTIFF(filename, rasterBand=1, substring='*_dem.tif', verbose=False)[source]¶

Read a single GeoTIFF file.

Parameters:

filename : str

Path (relative or absolute) to a GeoTIFF file.

Note

It is possible to read a GeoTIFF file directly from within a compressed archive without extracting it by using the /vsiPREFIX/.

For instance, to read ‘my.tif’ from a subdirectory within ‘my.zip’ use: /vsizip/my.zip/subdirectory/my.tif if path to ‘my.zip’ is relative or /vsizip//full_path_to/my.zip/subdirectory/my.tif if using the absolute path. Note the extra / after /vsizip/.

See http://www.gdal.org/gdal_virtual_file_systems.html for more details.

rasterBand : int

The band number to read.

substring : str

Compressed file formats may archive more than one file. The filename needs to be explicit as described in the note above. However, if the filename is not explicit, this function will try to find the correct file based on substring or regular expression matching.

Note

Being explicit is faster (but not by much)

Returns:

GeoTIFF

A populated (or empty) GeoTIFF instance with the elevation data.

save_GeoTIFF(filename, data, tlx, tly, dx, dy, epsg=None, proj4=None, dtype=<class 'numpy.int16'>, nodata=nan, rasterBand=1)[source]¶

Save data at a known projection as a GeoTIFF file.

Parameters:

filename : str

Path (relative or absolute) of the output GeoTIFF file.

data : ndarray

A 2d array of data to write to file.

tlx : float

Top-left-x (usually West-) coordinate of data.

tly : flaot

Top-left-y (usually North-) coordinate of data.

dx : float

Pixel size of data in the x direction, positive in the East direction.

dy : float

Pixel size of data in the y direction, positive in the North direction.

epsg : int

The target projection EPSG code. See the Geodetic Parameter Dataset Registry for more information.

proj4 : str

The target Proj4 string. If the EPSG code is unknown or a custom projection is required, a Proj4 string can be passed. See the Proj4 documentation for a list of general Proj4 parameters.

dtype : dtype

One of following dtypes should be used:

numpy.int16 (default)

numpy.int32

numpy.float32

numpy.float64

nodata : int or float or NaN

Value to register as the no-data value in the GeoTIFF file header. By default this is set to numpy.nan. Other common options are -9999, -12345 or any other value that suits your purpose.

rasterband : int

The band number to write. (default is 1)

Notes

Examples of some EPSG codes and their equivalent Proj4 strings are:

4326   -> '+proj=longlat +datum=WGS84 +no_defs'

32636  -> '+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs'

102009 -> '+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40
           +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m
           +no_defs'

and so on … See the Geodetic Parameter Dataset Registry for more information.

xy2gridpoints(x, y, grid_shape, grid_extent)[source]¶: Get the grid coordinates of the data coordinates in x, y

Previous topic

Next topic

This Page

pySW4.utils.geo module¶