Skip to content

Commit

Permalink
Minor changes for Python3.x
Browse files Browse the repository at this point in the history
Revisions to code so it can be used in Python3.x
  • Loading branch information
elbeejay authored Aug 12, 2019
1 parent 823c4dd commit 10b6dd0
Show file tree
Hide file tree
Showing 3 changed files with 63 additions and 64 deletions.
4 changes: 2 additions & 2 deletions rivamap/delineate.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ def extractCenterlines(orient, psi):

# Find maxima along local orientation
nms = np.zeros(psi.shape)
for q, (di, dj) in zip(range(4), ((1, 0), (1, 1), (0, 1), (-1, 1))):
for q, (di, dj) in zip(list(range(4)), ((1, 0), (1, 1), (0, 1), (-1, 1))):
for i, j in zip(*np.nonzero(np.logical_and(Q == q, mask))):
if psi[i, j] > psi[i + di, j + dj] and psi[i, j] > psi[i - di, j - dj]:
nms[i, j] = psi[i,j]
Expand Down Expand Up @@ -71,7 +71,7 @@ def thresholdCenterlines(nms, tLow=0.012, tHigh=0.12, bimodal=True):
# Find connected components that has at least one strong centerline pixel
strel = np.ones((3, 3), dtype=bool)
cclabels, numcc = ndlabel(centerlineCandidate, strel)
sumstrong = ndsum(strongCenterline, cclabels, range(1, numcc+1))
sumstrong = ndsum(strongCenterline, cclabels, list(range(1, numcc+1)))
centerlines = np.hstack((0, sumstrong > 0)).astype('bool')
centerlines = centerlines[cclabels]

Expand Down
119 changes: 59 additions & 60 deletions rivamap/georef.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,51 +18,51 @@ def __init__(self):
self.projection = None
self.geotransform = None
self.rasterXY = None


def loadGeoMetadata(filepath):
""" Reads metadata from a geotiff file
Inputs:
filepath -- path to the file
Returns:
gm -- metadata
"""
ds = gdal.Open(filepath)

if ds is None:
raise ValueError('Cannot read file')

gm = GeoMetadata()
gm.projection = ds.GetProjection()
gm.geotransform = ds.GetGeoTransform()
gm.rasterXY = (ds.RasterXSize, ds.RasterYSize)

if gm.projection is None:
warnings.warn('No projection found in the metadata')

if gm.geotransform is None:
warnings.warn('No geotransform found in the metadata')

# Close file
ds = None

return gm


def saveAsGeoTiff(gm, I, filepath):
""" Saves a raster image as a geotiff file
Inputs:
gm -- georeferencing metadata
I -- raster image
filepath -- path to the file
filepath -- path to the file
"""

if (I.shape[1] != gm.rasterXY[0]) and (I.shape[0] != gm.rasterXY[1]):
raise ValueError('Image size does not match the metadata')

DATA_TYPE = {
"uint8": gdal.GDT_Byte,
"int8": gdal.GDT_Byte,
Expand All @@ -75,28 +75,28 @@ def saveAsGeoTiff(gm, I, filepath):
}

driver = gdal.GetDriverByName('GTiff')

ds = driver.Create(filepath, I.shape[1], I.shape[0], 1, DATA_TYPE[I.dtype.name])

ds.SetGeoTransform(gm.geotransform)
ds.SetProjection(gm.projection)
ds.GetRasterBand(1).WriteArray(I)
ds.FlushCache()

ds = None


def pix2lonlat(gm, x, y):
""" Convers pixel coordinates into longitude and latitude
Inputs:
gm -- georeferencing metadata
x, y -- pixel coordinates
Returns:
lon, lat -- longitude and latitude
"""

sr = osr.SpatialReference()
sr.ImportFromWkt(gm.projection)
ct = osr.CoordinateTransformation(sr,sr.CloneGeogCS())
Expand All @@ -105,21 +105,21 @@ def pix2lonlat(gm, x, y):
lat_p = y*gm.geotransform[5]+gm.geotransform[3]

lon, lat, _ = ct.TransformPoint(lon_p, lat_p)

return lon, lat


def lonlat2pix(gm, lon, lat):
""" Convers longitude and latitude into pixel coordinates
Inputs:
gm -- georeferencing metadata
lon, lat -- longitude and latitude
Returns:
x, y -- pixel coordinates
"""

sr = osr.SpatialReference()
sr.ImportFromWkt(gm.projection)
ct = osr.CoordinateTransformation(sr.CloneGeogCS(),sr)
Expand All @@ -133,110 +133,109 @@ def lonlat2pix(gm, lon, lat):

def exportCSVfile(centerlines, widthMap, gm, filepath):
""" Exports (coordinate, width) pairs to a comma separated text file
Inputs:
centerlines -- a binary matrix that indicates centerline locations
widthMap -- estimated width at each spatial location (x,y)
gm -- georeferencing metadata
filepath -- path to the file
"""

centerlineWidth = widthMap[centerlines]
[row,col] = np.where(centerlines)
with open(filepath, 'wb') as csvfile:

with open(filepath, 'w') as csvfile:
writer = csv.writer(csvfile, delimiter=',')
writer.writerow(["width","lat","lon"])

for i in range(0, len(centerlineWidth)):
lon, lat = pix2lonlat(gm, col[i], row[i])
writer.writerow([centerlineWidth[i], lat, lon])

def exportShapeFile(centerlines, widthMap, gm, filepath):
""" Exports line segments to a ShapeFile
Inputs:
centerlines -- a binary matrix that indicates centerline locations
widthMap -- estimated width at each spatial location (x,y)
gm -- georeferencing metadata
filepath -- path to the file
"""

# initiate the shp writer
w = shapefile.Writer()
w = shapefile.Writer(filepath)
w.field('width', 'N', decimal=2)

R, C = centerlines.shape

# make a copy of the centerlines matrix since we'll be modifying it
c_copy = centerlines.copy()

# scan the centerline matrix skipping the boundary pixels
# convert neighbor centerline pixels into line segments
for r in range(1, R-1):
for c in range(1, C-1):
if c_copy[r,c]:
c_copy[r,c] = 0
if c_copy[r,c]:
c_copy[r,c] = 0
lon_orig, lat_orig = pix2lonlat(gm, c, r)

# connect neighbors (check diagonals first)
if c_copy[r-1, c-1]:
c_copy[r-1, c] = 0
c_copy[r, c-1] = 0
lon, lat = pix2lonlat(gm, c-1, r-1)
width = ( + widthMap[r-1, c-1]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])
w.line([[[lon_orig, lat_orig], [lon, lat]]])

if c_copy[r-1, c+1]:
c_copy[r-1, c] = 0
c_copy[r, c+1] = 0
lon, lat = pix2lonlat(gm, c+1, r-1)
width = (widthMap[r, c] + widthMap[r-1, c+1]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])
w.line([[[lon_orig, lat_orig], [lon, lat]]])

if c_copy[r+1, c-1]:
c_copy[r+1, c] = 0
c_copy[r, c-1] = 0
lon, lat = pix2lonlat(gm, c-1, r+1)
width = (widthMap[r, c] + widthMap[r+1, c-1]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])
w.line([[[lon_orig, lat_orig], [lon, lat]]])

if c_copy[r+1, c+1]:
c_copy[r+1, c] = 0
c_copy[r, c+1] = 0
lon, lat = pix2lonlat(gm, c+1, r+1)
width = (widthMap[r, c] + widthMap[r+1, c+1]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])
w.line([[[lon_orig, lat_orig], [lon, lat]]])

if c_copy[r-1, c]:
lon, lat = pix2lonlat(gm, c, r-1)
width = (widthMap[r, c] + widthMap[r-1, c]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])
w.line([[[lon_orig, lat_orig], [lon, lat]]])

if c_copy[r+1, c]:
lon, lat = pix2lonlat(gm, c, r+1)
width = (widthMap[r, c] + widthMap[r+1, c]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])
w.line([[[lon_orig, lat_orig], [lon, lat]]])

if c_copy[r, c-1]:
lon, lat = pix2lonlat(gm, c-1, r)
width = (widthMap[r, c] + widthMap[r, c-1]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])
w.line([[[lon_orig, lat_orig], [lon, lat]]])

if c_copy[r, c+1]:
lon, lat = pix2lonlat(gm, c+1, r)
width = (widthMap[r, c] + widthMap[r, c+1]) / 2
w.record(width)
w.line(parts=[[[lon_orig, lat_orig], [lon, lat]]])

w.save(filepath)

w.line([[[lon_orig, lat_orig], [lon, lat]]])

w.close()
4 changes: 2 additions & 2 deletions rivamap/singularity_index.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ def _createFilters(self):
theta3 = 2*np.pi/3

# Create a meshgrid for second order derivatives
X, Y = np.meshgrid(range(-ksize2,ksize2+1), range(-ksize2,ksize2+1))
X, Y = np.meshgrid(list(range(-ksize2,ksize2+1)), list(range(-ksize2,ksize2+1)))
u1 = X*np.cos(theta1) - Y*np.sin(theta1)
u2 = X*np.cos(theta2) - Y*np.sin(theta2)
u3 = X*np.cos(theta3) - Y*np.sin(theta3)
Expand Down Expand Up @@ -107,7 +107,7 @@ def applyMMSI(I1, filters, togglePolarity=False, narrow_rivers=True):

# Compute the multiscale singularity index
for s in range(0, filters.nrScales):
print "Processing scale: " + str(s)
print("Processing scale: " + str(s))

# Downscale the image to the current scale. We use a pyramid instead of
# increasing the sigma and size of the kernels for efficiency
Expand Down

0 comments on commit 10b6dd0

Please sign in to comment.