POST
|
Has anyone heard if Arc will be supporting the OGC GeoPackage format: http://www.opengeospatial.org/standards/requests/115 Thanks.
... View more
11-26-2013
10:56 AM
|
0
|
1
|
3186
|
POST
|
Thanks, that is the one I created because I couldn't find one already up.
... View more
03-16-2012
01:31 PM
|
0
|
0
|
213
|
POST
|
Is there an idea up? I would like to vote on it. I couln't find one.
... View more
03-16-2012
12:52 PM
|
0
|
0
|
1742
|
POST
|
Any reason why the name of the library was changed from: libfgdblinuxrtl.so to libfgdbunixrtl.so ? Just curious, it messed with my builds.
... View more
03-16-2012
12:51 PM
|
0
|
0
|
472
|
POST
|
I managed to hack this together at one point, but found it was too slow for a web service. I rewrote the script using gdal for raster access and it was hundreds of times faster. The code for the arcpy script is below, with an example output graph. The commented code was for grabbing points from a feature set (for a toolbox or gp service). This ugly, but it works. disclaimer: This is hacked, old and ugly. sorry.
from polyline import *
import matplotlib.pyplot as plt
def fetch_cell_value(x, y):
xy_string = "{0} {1}".format(x, y)
band = "1"
arcpy.GetCellValue_management(r, xy_string, band)
results = arcpy.GetMessages()
lines = results.split('\n')
value = float(lines[2])
return value
#for arc script
#fs = arcpy.GetParameter(0)
#rows = arcpy.SearchCursor(fs)
# Assume feature set has only one feature and one shape.
#row = rows.next()
#shape = row.getValue("SHAPE")
#arcpoints = shape.getPart(0)
#points = list()
#for p in arcpoints:
# points.append([p.X, p.Y])
raster = "c:/data/someraster.tif"
r = arcpy.Raster(raster)
#test points
points = [[290626.474, 4875792.8594],
[282358.29,4863026.666],
[297836.3843, 4870699.62]]
p = poly_line(points)
n_points = 100
#if p.length() / n_points < cell_size:
# skip = cell_size
#else:
# skip = int(p.length() / n_points)
skip = p.length() / n_points
dists = list()
i = 0
while i * skip < p.length():
dists.append(i * skip)
i += 1
values = list()
for d in dists:
x,y = p.d_line(d)
v = fetch_cell_value(x,y)
values.append(v)
print("Time:{0}".format(time.clock() - start))
plt.plot(dists, values)
plt.xlabel("Distance (m)")
plt.ylabel("Elevation (m)")
plt.title("Elevation Profile")
plt.grid(True)
plt.show()
The poly_line is a little class I wrote to handle densifying lines. It assumes euclidean geometry (planar) so doesn't handle earth curvature. You could probably use some type of densify function for an arc type polyline and use the points in the shape. I didn't try it. poly_line:
import math
##
# Class representing a segment of a polyline
#
class segment:
p1 = [0,0]
p2 = [0,0]
##
# @brief Construct a segment
#
def __init__(self, p1, p2):
self.p1 = p1
self.p2 = p2
##
# @brief Get the length of the segment
#
# @todo Test for /0
#
def length(self):
a = self.p2[0] - self.p1[0]
b = self.p2[1] - self.p1[1]
l = math.hypot(a,b)
return l
##
# @brief Get the slope of the segment
#
# @todo Test for /0
#
def m(self):
a = self.p2[1] - self.p1[1]
b = self.p2[0] - self.p1[0]
m = a / b
return m
##
# @brief Get the y intercept of the segment
#
# @return y intercept
def b(self):
b = self.p2[1] - (self.m() * self.p2[0])
return b
##
# @brief Get the coordinate of a point some distance from point 1
#
# @param d distance from point 1
# @return x and y coordinate
def d_line(self, d):
if d == 0:
return self.p1[0], self.p1[1]
if d > self.length():
return None
m = self.m()
u = d / math.sqrt(m * m + 1)
if self.p2[0] < self.p1[0]:
u = -u
v = (m * d) / math.sqrt(m * m + 1)
if self.p2[1] < self.p1[1]:
v = -v
return (self.p1[0] + u, self.p1[1] + v)
##
# Class representing a polyline composed of segments
#
# @see segment
#
class poly_line:
segments = list()
def __init__(self, points):
for i in range(len(points) - 1):
s = segment(points, points[i+1])
self.segments.append(s)
##
# @brief Get the length of all the segments
#
# @return total length
#
def length(self):
l = 0.0
for s in self.segments:
l += s.length()
return l
##
# @brief Get the coordinate of a point some distance from the polyline
# origin
#
# @param d distance from the origin
# @return x and y coordinate for the point
#
def d_line(self, d):
if d == 0:
return self.segments[0].d_line(0)
l = 0.0
for s in self.segments:
if l >= d:
break
dd = d - l
l += s.length()
seg = s
return seg.d_line(dd)
I just ran it and wrote the attached output, but no guarentees it will work for anyone. I think that the service ESRI wrote used lower level raster access, because of it's speed. I got comparable speed with the gdal version. Sorry if I didn't help, but it is possible. I stopped working on it, but it may help others. I will watch the thread and try to answer questions. kss
... View more
02-29-2012
10:54 AM
|
1
|
0
|
816
|
POST
|
Next time, use gdalwarp, it does all of that together and the -multi flag uses parallel processing.
... View more
01-05-2012
01:52 PM
|
0
|
0
|
193
|
POST
|
GDAL has a command line tool, gdaldem. You can write a geotiff as the output.
gdaldem hillshade in.tif out.tif
... View more
12-19-2011
02:10 PM
|
0
|
0
|
279
|
POST
|
Are there actual cells with those values? Try computing statistics on them.
... View more
12-19-2011
02:03 PM
|
0
|
0
|
413
|
POST
|
On another quick note, the del statement doesn't actually delete objects, it dereferences the name from the object, which *may* result in deletion. The garbage collector decides that stuff. There are probably references to the object somewhere deep in arcpy or something. Again, hard tellin' not knowin' kss
... View more
12-12-2011
06:25 AM
|
0
|
0
|
970
|
POST
|
I am definitely not sure, so I am totally guessing. I suspect that Delete_management isn't made to interact with python/arcpy objects directly, same as many arcpy functions. It's possible that esri is slowly adding support to act on these, but if it were supported, I don't see any reason for an in_memory feature class/set. I tried this and got okay results:
line_shapefile = "C:\Users\jhook\Desktop\lines.shp"
output = "in_memory/output"
for i in range(0, 100, 1):
geom_obj = arcpy.CopyFeatures_management(line_shapefile, output)
arcpy.Delete_management(geom_obj)
I am guessing that most arcpy functions act on file-like objects, so by creating in_memory file like objects, people could use those functions. I don't know anything about future plans or anything. this is all guessing. I would stick with in_memory feature classes for now.
... View more
12-12-2011
06:18 AM
|
0
|
0
|
970
|
POST
|
What were your inputs? What is the spatial reference of the source raster?
... View more
12-02-2011
10:57 AM
|
0
|
0
|
140
|
POST
|
Can you post a code snippet? Are you running in arcgis python window?
... View more
12-02-2011
10:10 AM
|
0
|
0
|
970
|
POST
|
Both fgdb and esri's shape access suffer from size/performance correlation. It is strange, because when I use ogr to add fields, insertion is constant. Good example of 'Schlemiel the Painter's algorithm'[0]. Never thought I would see an example of that now a days. Attached is a graph of insertion times over the course of 1000 AddField calls. I also tested the copy of the mem to gdb, and that seems to be the best way to go in arc: Output of test over 1000 Add Field calls: C:\Users\kshannon\Desktop>c:\Python26\ArcGIS10.0\python.exe add_field_test.py
Creating feature classes...
Feature classes created.
GDB(Total:940.140323837,Max(924):1.72044671624,Min(60):0.582889467023,Average:0.940140323837
SHP(Total:352.133845124,Max(994):1.01701894889,Min(2):0.0293058087336,Average:0.352133845124
MEM(Total:52.2060125065,Max(994):1.01701894889,Min(2):0.0293058087336,Average:0.0522060125065
MEM to GDB:2.62204624813
OGR(Total:37.934927268,Max(994):1.01701894889,Min(1):0.0222360889709,Average:0.037934927268
Deleting feature classes...
Feature classes deleted. The values in parentheses are the iteration at which the min/max occurred. Thanks for the help. [0] http://en.wikipedia.org/wiki/Schlemiel_the_Painter's_algorithm
... View more
12-01-2011
06:34 AM
|
0
|
0
|
208
|
POST
|
I remember doing something like this for drive time polygons. I think I called CreateFeatureclass() and used the layer and sub layer as a template, then appended the data from the polygon layer.
arcpy.Solve_na(s_name)
arcpy.CreateFeatureclass_management(working_dir, l_name, "POLYGON",
"{0}\Polygons".format(s_name),
"", "", "", "", "", "", "")
arcpy.Append_management("{0}\Polygons".format(s_name),
out_name, "NO_TEST", "", "")
Not sure if that is what you are looking for, or if I did it right. It was in the middle of a multiprocess run where I split up a NA layer and ran it on 8 threads, so you could probably just use CopyFeatures().
... View more
12-01-2011
06:09 AM
|
0
|
0
|
184
|
Title | Kudos | Posted |
---|---|---|
1 | 02-29-2012 10:54 AM |
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:24 AM
|