lib.geolib module

Geospatial functions for rasters, vectors.

pygeotools.lib.geolib.CE90(x_offset, y_offset)

pygeotools.lib.geolib.LE90(z_offset)

class pygeotools.lib.geolib.ProjBox(bbox, epsg)

Object containing bbox geom and srs

Used for automatic projection determination

class pygeotools.lib.geolib.Site(name, extent, srs=None, refdem_fn=None)

pygeotools.lib.geolib.applyGeoTransform(inX, inY, geoTransform)

pygeotools.lib.geolib.bbox2geom(bbox, a_srs=<Mock name='mock.SpatialReference()' id='140299428549904'>, t_srs=None)

pygeotools.lib.geolib.bilinear(px, py, band_array, gt)

Bilinear interpolated point at(px, py) on band_array

pygeotools.lib.geolib.block_stats(x, y, z, ds, stat='median', bins=None)

Compute points on a regular grid (matching input GDAL Dataset) from scattered point data using specified statistic

Wrapper for scipy.stats.binned_statistic_2d

Note: this is very fast for mean, std, count, but bignificantly slower for median

pygeotools.lib.geolib.block_stats_grid(x, y, z, ds, stat='median')

Fill regular grid (matching input GDAL Dataset) from scattered point data using specified statistic

pygeotools.lib.geolib.block_stats_grid_gen(x, y, z, res=None, srs=None, stat='median')

pygeotools.lib.geolib.cT_helper(x, y, z, in_srs, out_srs)

pygeotools.lib.geolib.clip_raster_by_shp(dem_fn, shp_fn)

pygeotools.lib.geolib.clip_shp(shp_fn, extent)

pygeotools.lib.geolib.copyproj(src_fn, dst_fn, gt=True)

Copy projection and geotransform from one raster file to another

pygeotools.lib.geolib.corner_extent(ul, ll, ur, lr)

Get min/max extent based on corner coord

pygeotools.lib.geolib.dd2dm(dd)

Convert decimal to degrees, decimal minutes

pygeotools.lib.geolib.dd2dms(dd)

Convert decimal degrees to degrees, minutes, seconds

pygeotools.lib.geolib.dem_geoid(dem_fn)

pygeotools.lib.geolib.dem_geoid_offsetgrid(dem_fn)

pygeotools.lib.geolib.dem_geoid_offsetgrid_ds(ds, out_fn=None)

pygeotools.lib.geolib.dm2dd(d, m)

Convert degrees, decimal minutes to decimal degrees

pygeotools.lib.geolib.dms2dd(d, m, s)

Convert degrees, minutes, seconds to decimal degrees

pygeotools.lib.geolib.dms2dd_str(dms_str, delim=' ', fmt=None)

pygeotools.lib.geolib.ds2mask(v_ds, r_ds)

pygeotools.lib.geolib.ds_IsEmpty(ds)

Check to see if dataset is empty after warp

pygeotools.lib.geolib.ds_cT(ds, x, y, xy_srs=<Mock name='mock.SpatialReference()' id='140299428549904'>)

Convert input point coordinates to map coordinates that match input dataset

pygeotools.lib.geolib.ds_extent(ds, t_srs=None)

Return min/max extent of dataset based on corner coordinates

xmin, ymin, xmax, ymax

If t_srs is specified, output will be converted to specified srs

pygeotools.lib.geolib.ds_fitplane(ds)

Fit a plane to values in GDAL Dataset

pygeotools.lib.geolib.ds_geom(ds, t_srs=None)

Return dataset bbox envelope as geom

pygeotools.lib.geolib.ds_geom_extent(ds, t_srs=None)

pygeotools.lib.geolib.ds_geom_intersection(ds_list, **kwargs)

pygeotools.lib.geolib.ds_geom_intersection_extent(ds_list, **kwargs)

pygeotools.lib.geolib.ds_geom_union(ds_list, **kwargs)

pygeotools.lib.geolib.ds_geom_union_extent(ds_list, **kwargs)

pygeotools.lib.geolib.ecef2ll(x, y, z)

pygeotools.lib.geolib.ell2geoid(lon, lat, z=0.0, geoid=<Mock name='mock.SpatialReference()' id='140299428549904'>)

pygeotools.lib.geolib.extent_compare(e1, e2, precision=0.001)

pygeotools.lib.geolib.extent_round(extent, precision=0.001)

pygeotools.lib.geolib.extract_profile(ds, geom, dl=None, km=False)

pygeotools.lib.geolib.fitPlaneLSQ(XYZ)

Fit a plane to input point data using LSQ

pygeotools.lib.geolib.fitPlaneSVD(XYZ)

Fit a plane to input point data using SVD

pygeotools.lib.geolib.fn2mask(fn, r_ds)

pygeotools.lib.geolib.gdaldem_aspect(fn)

pygeotools.lib.geolib.gdaldem_mem_ds(ds, processing='hillshade', returnma=False, computeEdges=False)

Wrapper for gdaldem functions

Uses gdaldem API, requires GDAL v2.1+

pygeotools.lib.geolib.gdaldem_mem_ma(ma, ds=None, res=None, extent=None, srs=None, processing='hillshade', returnma=False, computeEdges=False)

Wrapper to allow gdaldem calculations for arbitrary NumPy masked array input Untested, work in progress placeholder Should only need to specify res, can caluclate local gt, cartesian srs

pygeotools.lib.geolib.gdaldem_slope(fn)

pygeotools.lib.geolib.gdaldem_wrapper(fn, product='hs', returnma=True, verbose=True)

Wrapper for gdaldem functions

Note: gdaldem is directly avaialable through API as of GDAL v2.1

https://trac.osgeo.org/gdal/wiki/rfc59.1_utilities_as_a_library

This function is no longer necessry, and will eventually be removed.

pygeotools.lib.geolib.gen_proj_list()

Create list of projections with cascading preference

pygeotools.lib.geolib.geoid2ell(lon, lat, z=0.0, geoid=<Mock name='mock.SpatialReference()' id='140299428549904'>)

pygeotools.lib.geolib.geom2localortho(geom)

Convert existing geom to local orthographic projection

Useful for local cartesian distance/area calculations

pygeotools.lib.geolib.geom2mask(geom, r_ds)

pygeotools.lib.geolib.geom2mask_PIL(geom, ds)

pygeotools.lib.geolib.geom2shp(geom, out_fn, fields=False)

Write out a new shapefile for input geometry

pygeotools.lib.geolib.geom_dup(geom)

Create duplicate geometry

Needed to avoid segfault when passing geom around. See: http://trac.osgeo.org/gdal/wiki/PythonGotchas

pygeotools.lib.geolib.geom_extent(geom)

pygeotools.lib.geolib.geom_intersection(geom_list, **kwargs)

pygeotools.lib.geolib.geom_transform(geom, t_srs)

Transform a geometry in place

pygeotools.lib.geolib.geom_union(geom_list, **kwargs)

pygeotools.lib.geolib.geom_wh(geom)

Compute width and height of geometry in projected units

pygeotools.lib.geolib.getUTMsrs(geom)

pygeotools.lib.geolib.getUTMzone(geom)

Determine UTM Zone for input geometry

pygeotools.lib.geolib.get_GeoidOffset(lon, lat)

pygeotools.lib.geolib.get_HAE(lon, lat)

pygeotools.lib.geolib.get_MSL(lon, lat)

pygeotools.lib.geolib.get_NED(lon, lat)

pygeotools.lib.geolib.get_OpenElevation(lon, lat)

pygeotools.lib.geolib.get_center(ds, t_srs=None)

Get center coordinates of GDAL Dataset

pygeotools.lib.geolib.get_dem_mosaic_cmd(fn_list, o, fn_list_txt=None, tr=None, t_srs=None, t_projwin=None, georef_tile_size=None, threads=None, tile=None, stat=None)

Create ASP dem_mosaic command Useful for spawning many single-threaded mosaicing processes

pygeotools.lib.geolib.get_ds_srs(ds)

Get srs object for GDAL Datset

pygeotools.lib.geolib.get_geoid_offset(ds, geoid_srs=<Mock name='mock.SpatialReference()' id='140299428549904'>)

Return offset for center of ds

Offset is added to input (presumably WGS84 HAE) to get to geoid

Note: requires vertical offset grids in proj share dir - see earlier note

pygeotools.lib.geolib.get_geoid_offset_ll(lon, lat, geoid_srs=<Mock name='mock.SpatialReference()' id='140299428549904'>)

pygeotools.lib.geolib.get_outline(ds, t_srs=None, scale=1.0, simplify=False, convex=False)

Generate outline of unmasked values in input raster

get_outline is an attempt to reproduce the PostGIS Raster ST_MinConvexHull function

Could potentially do the following: Extract random pts from unmasked elements, get indices, Run scipy convex hull, Convert hull indices to mapped coords

See this: http://stackoverflow.com/questions/3654289/scipy-create-2d-polygon-mask

This generates a wkt polygon outline of valid data for the input raster

Want to limit the dimensions of a, as notmasked_edges is slow: a = iolib.ds_getma_sub(ds, scale=scale)

pygeotools.lib.geolib.get_proj(geom, proj_list=None)

Determine best projection for input geometry

pygeotools.lib.geolib.get_res(ds, t_srs=None, square=False)

Get GDAL Dataset raster resolution

pygeotools.lib.geolib.get_res_stats(ds_list, t_srs=None)

Return resolution stats for an input dataset list

pygeotools.lib.geolib.get_xy_1D(ds, stride=1, getval=False)

Return 1D arrays of x and y map coordinates for input GDAL Dataset

pygeotools.lib.geolib.get_xy_grids(ds, stride=1, getval=False)

Return 2D arrays of x and y map coordinates for input GDAL Dataset

pygeotools.lib.geolib.get_xy_ma(bma, gt, stride=1, origmask=True, newmask=None)

Return arrays of x and y map coordinates for input array and geotransform

pygeotools.lib.geolib.gt_corners(gt, nx, ny)

Get corner coordinates based on input geotransform and raster dimensions

pygeotools.lib.geolib.gt_extent(gt, nx, ny)

pygeotools.lib.geolib.heatmap(gdf, res=1000, out_fn='heatmap.tif', return_ma=False)

pygeotools.lib.geolib.heatmap_fn(fn, res=1000, return_ma=False)

pygeotools.lib.geolib.invertGeoTransform(geoTransform)

pygeotools.lib.geolib.itrf2ll(x, y, z)

pygeotools.lib.geolib.line2pts(geom, dl=None)

Given an input line geom, generate points at fixed interval

Useful for extracting profile data from raster

pygeotools.lib.geolib.ll2ecef(lon, lat, z=0.0)

pygeotools.lib.geolib.ll2itrf(lon, lat, z=0.0)

pygeotools.lib.geolib.ll2local(lon, lat, z=0, local_srs=None)

pygeotools.lib.geolib.ll2nps(lon, lat, z=0.0)

pygeotools.lib.geolib.ll2sps(lon, lat, z=0.0)

pygeotools.lib.geolib.lldist(pt1, pt2)

pygeotools.lib.geolib.localortho(lon, lat)

Create srs for local orthographic projection centered at lat, lon

pygeotools.lib.geolib.localortho_ds(ds)

pygeotools.lib.geolib.localtmerc(lon, lat)

pygeotools.lib.geolib.localtmerc_ds(ds)

pygeotools.lib.geolib.lon180to360(lon)

Convert longitude from (-180, 180) to (0, 360)

pygeotools.lib.geolib.lon360to180(lon)

Convert longitude from (0, 360) to (-180, 180)

pygeotools.lib.geolib.lyr2mask(lyr, r_ds)

pygeotools.lib.geolib.lyr_extent(lyr)

pygeotools.lib.geolib.lyr_proj(lyr, t_srs, preserve_fields=True)

Reproject an OGR layer

pygeotools.lib.geolib.ma_fitplane(bma, gt=None, perc=(2, 98), origmask=True)

Fit a plane to values in input array

pygeotools.lib.geolib.ma_fitpoly(bma, order=1, gt=None, perc=(2, 98), origmask=True)

Fit a plane to values in input array

pygeotools.lib.geolib.mapToPixel(mX, mY, geoTransform)

Convert map coordinates to pixel coordinates based on geotransform

Accepts float or NumPy arrays

GDAL model used here - upper left corner of upper left pixel for mX, mY (and in GeoTransform)

pygeotools.lib.geolib.map_interp(bma, gt, stride=1, full_array=True)

pygeotools.lib.geolib.mem_ds(res, extent, srs=None, dtype=<Mock name='mock.GDT_Float32' id='140299428550736'>)

Create a new GDAL Dataset in memory

Useful for various applications that require a Dataset

pygeotools.lib.geolib.mem_ds_copy(ds_orig)

pygeotools.lib.geolib.nps2geoid(x, y, z=0.0, geoid=<Mock name='mock.SpatialReference()' id='140299428549904'>)

pygeotools.lib.geolib.nps2ll(x, y, z=0.0)

pygeotools.lib.geolib.pad_extent(extent, perc=0.1, width=None, uniform=False)

pygeotools.lib.geolib.pixelToMap(pX, pY, geoTransform)

Convert pixel coordinates to map coordinates based on geotransform

Accepts float or NumPy arrays

GDAL model used here - upper left corner of upper left pixel for mX, mY (and in GeoTransform)

pygeotools.lib.geolib.polyfit2d(x, y, z, order=1)

pygeotools.lib.geolib.polyval2d(x, y, m)

pygeotools.lib.geolib.pt_within_extent(x, y, extent)

pygeotools.lib.geolib.raster_shpclip(r_fn, shp_fn, extent='raster', bbox=False, pad=None, invert=False, verbose=False)

Clip an input raster by input polygon shapefile for given extent

pygeotools.lib.geolib.res_compare(r1, r2, precision=0.001)

pygeotools.lib.geolib.round_nearest(x, a)

pygeotools.lib.geolib.sample(ds, mX, mY, xy_srs=None, bn=1, pad=0, min_samp_perc=50, circ=False, count=False)

Sample input dataset at map coordinates

This is a generic sampling function, and will return value derived from window (dimensions pad*2+1) around each point

By default, assumes input map coords are identical to ds srs. If different, specify xy_srs to enable conversion.

pygeotools.lib.geolib.scale_ps(lat)

This function calculates the scaling factor for a polar stereographic projection (ie. SSM/I grid) to correct area calculations. The scaling factor is defined (from Snyder, 1982, Map Projections used by the U.S. Geological Survey) as:

k = (mc/m)*(t/tc), where:

m = cos(lat)/sqrt(1 - e2*sin(lat)^2) t = tan(Pi/4 - lat/2)/((1 - e*sin(lat))/(1 + e*sin(lat)))^(e/2) e2 = 0.006693883 is the earth eccentricity (Hughes ellipsoid) e = sqrt(e2) mc = m at the reference latitude (70 degrees) tc = t at the reference latitude (70 degrees)

The ratio mc/tc is precalculated and stored in the variable m70_t70.

From Ben Smith PS scale m file (7/12/12)

pygeotools.lib.geolib.scale_ps_ds(ds)

pygeotools.lib.geolib.shp2array(shp_fn, r_ds=None, res=None, extent=None, t_srs=None)

Rasterize input shapefile to match existing raster Dataset (or specified res/extent/t_srs)

pygeotools.lib.geolib.shp2geom(shp_fn)

Extract geometries from input shapefile

Need to handle multi-part geom: http://osgeo-org.1560.x6.nabble.com/Multipart-to-singlepart-td3746767.html

pygeotools.lib.geolib.shp_dict(shp_fn, fields=None, geom=True)

Get a dictionary for all features in a shapefile

Optionally, specify fields

pygeotools.lib.geolib.shp_fieldnames(lyr)

pygeotools.lib.geolib.sps2geoid(x, y, z=0.0, geoid=<Mock name='mock.SpatialReference()' id='140299428549904'>)

pygeotools.lib.geolib.sps2ll(x, y, z=0.0)

pygeotools.lib.geolib.sps2local(x, y, z=0, local_srs=None)

pygeotools.lib.geolib.srs_check(ds)

Check validitiy of Dataset srs

Return True if srs is properly defined

pygeotools.lib.geolib.test_elev_api(lon, lat)

pygeotools.lib.geolib.tp2wgs(x, y, z)

pygeotools.lib.geolib.wgs2tp(x, y, z)

pygeotools.lib.geolib.wgs84_to_egm96(dem_ds, geoid_dir=None)

pygeotools.lib.geolib.wraplon(lon)

pygeotools.lib.geolib.xy2geom(x, y, t_srs=None)

Convert x and y point coordinates to geom