New answers tagged python
0
This is most easily solve using the equidistant azimuthal projection.
Guess some point for D, e.g., A. Transform A, B, C, to equidistant
azimuthal projection using D as a center. In projected space, solve for
D' (i.e., D' lies on BC and is a distance s from A). Project D' back to
lat,lon and update D with this position. Repeat. This converges
...
0
The ESRI example do indeed not work at all. This method works fine:
for row in rows:
feat = row.getValue(shapefieldname)
print "Feature %i: numpoints=%i" % (row.getValue(desc.OIDFieldName), feat.pointCount)
print 'Number of parts: ', feat.partCount
partnum = 0;
while partnum < feat.partCount:
print "Part %i:" % partnum
...
0
It looks like an indentation problem, plus had an extra ")" and was missing an "r" on a file path - try ...
import sys
import os
import atexit
import grass.script as grass
def cleanup():
pass
def main():
for i in range (2,4):
out = "reclasstest" + str(i)
rulesvar = r"C:\Grass_data\datareclassfile" + str(i)
...
1
Using ListLayers will help you in this task, especially if you want to apply this changes to all of the layers in the map document. In the question, arcpy is only seeing a list of text values, not layers. This will help get you to the FOR loop for the layers.
Set the symbology layers as you want it and save it. Reference it with a absolute path or set ...
4
While you can do it with the standalone installed copy of QGIS I would highly recommend installing QGIS and QGIS-dev using the OSGeo4W installer. OSGeo is easier to update and lets you get the latest stuff.
That will install Python into C:\OSGeo4W\app\Python27 and let you run Python from the OSGeo4W shell. python.exe lives in C:\OSGeo4W\bin.
I have some ...
2
You say "I'm stuck in iterating all the input folders" so that is the Question this Answer addresses.
Iterating through ArcGIS workspaces becomes much easier at ArcGIS 10.1 SP1.
However, it could be done using earlier versions of arcpy and arcgisscripting using os.walk as described here or glob and other techniques as described here.
2
You can do it by converting your bounding box to pixel offset and size:
from osgeo import gdal
def map_to_pixel(mx,my,gt):
#Convert from map to pixel coordinates.
#Only works for geotransforms with no rotation.
px = int((mx - gt[0]) / gt[1]) #x pixel
py = int((my - gt[3]) / gt[5]) #y pixel
return px,py
def ...
2
How about using the data access module? It looks like you can start an edit session with this module.
A few caveats:
I have not tried this module and the am not sure if it is 10.0
compatible. (New in 10.1?)
Example 1 shows the use of a with statement. This
is a great paradigm to implement as it handles potential exceptions
well.
You might be able to ...
0
Intersection of two paths given start points and bearings could be used. You would need to check the resulting point to make sure it is between B and C. Since there's no spheroid flattening, I'm not sure how close it would match up with Vincenty's formula.
Formula:
d12 = 2.asin( √(sin²(Δφ/2) + cos(φ1).cos(φ2).sin²(Δλ/2)) )
φ1 = acos( sin(φ2) − ...
0
You left off the search distance argument to SelectLayerByLocation. This is important -- for some reason, without it, a PointGeometry object will never intersect anything, at least in my limited testing.
Change:
arcpy.SelectLayerByLocation_management(lyr, "INTERSECT", pointGeom)
To:
arcpy.SelectLayerByLocation_management(lyr, "INTERSECT", pointGeom, "%d ...
1
I think your loop at the end is the problem. You will go though your list of matched flags, and for each one select the exp then switching then for the next item you would add the matched and switch again. You may want to simply move the switch selection out of the loop, that way you will select everything in matched then switch. Alternatively you could ...
0
To join line segments you can assign a from-node# and a to-node# to each line segment, where the node# is a hash of x and y rounded off. For example, a point at 123.45,567.8 might be computed to be 123568. This implies a tolerance of 1. You could generalize this to round off to the closest multiple of the tolerance.
Once you have From and To node#'s, you ...
0
You can use effectively Shapely, and Fiona to read a shapefile for example:
import fiona
# open a line shapefile
file = fiona.open('lines.shp')
# first element of the shapefile
first = file.next
print first
{'geometry': {'type': 'LineString', 'coordinates': [(203317.23, 90448.75), (203679.62, 90105.68), (203882.57, 89902.74), (204143.49, 89641.81), ...
1
It seems that Shapely can answer what you are expecting.
See the manual.
Create lines using native objects.
For what you need, only loop ;) to get segments of a line.
For getting intersection point, use operators like intersection
I can't say that is efficient but I suppose that because Shapely uses Geos behind, it's quite fast.
1
Seems like this is a "feature"...http://hub.qgis.org/issues/777
When the provider loops through the features, python stores the geometry, but as soon as the loop continues to the next iteration, the feature gets destroyed, and with it the geometry gets destroyed--passing invalid memory when trying to combine the geometries. So to save the geometries to ...
3
As far as I know, there is no way to pass objects between .net and python, but you can use a common data structure like JSON to pass data in string form back and forth and restructure using something like JSON.NET which has a method of deserializing the raw JSON string into a dictionary. It should be easy to find a library that does this on the python end ...
1
you need to add all the geometries before the addMapLayers
v_layer = QgsVectorLayer("LineString", "cable", "memory")
pr = v_layer.dataProvider()
# first element
seg = QgsFeature()
seg.setGeometry(QgsGeometry.fromPolyline([line_start, geomPoint]))
pr.addFeatures( [ seg ] )
v_layer.updateExtents()
# second element
seg = QgsFeature()
...
0
As @bananafish noted, no should be nos. You might also want to set i=0 before the loop.
nos = QgsFeature()
while nodedata.nextFeature(nos):
i = i + 1
print str(nos.gid())
4
Both Michael Markieta and gm70560 are correct. If you are running running large geoprocessing tasks, I would definitely do it via a stand-alone python script, preferably launched from the command line and not an IDE. For this sort of task, the overhead of importing ArcPy is well worth it.
However, a small task, especially one which requires user input, is ...
1
In the while loop, no doesn't refer to anything, you need to assign it to the current feature. I don't know how the QGIS API works, so here's a somewhat ugly approach:
while True:
try:
no = nodedata.nextFeature(nos)
except:
break
print str(no.gid())
10
In order to get the coordinates in decimal degrees, the data needs to be reprojected to WGS84.
import ogr, osr
driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open('testpoint.shp', 0)
lyr = shp.GetLayer()
feat = lyr.GetNextFeature()
geom = feat.GetGeometryRef()
# Transform from Web Mercator to WGS84
sourceSR = lyr.GetSpatialRef()
targetSR = ...
1
I can spot two things:
Change every instance of EPSG:900913 to EPSG:3857 -- EPSG:900913 (Or EPSG:GOOGLE if you squint hard enough) never existed in the EPSG database. It's a long story, but EPSG:3857 is the correct code.
It looks like you have one too many zeroes when you initialize your minLat and minLong variables. Compare it to your new ...
2
I think the crux of the question here is which tasks in your workflow are not really ArcGIS dependent? Obvious candidates include tabular and raster operations. If the data must start and end within a gdb or some other ESRI format, then you need to figure out how to minimize the cost of this reformat (i.e., minimize the number of round trips) or even justify ...
1
This functionality is covered in "Points2One" and "Points to Paths" plugins. You can take code from there.
0
So the answer to this took some time to find but ended up trivial but not straightforward to spot (to me at least). Damn you PATH...
When my plugin was running, I noted that if I ran os.getcwd() within my plugin, it did not actually output the plugin directory, but instead "C:\Program Files (x86)\QGIS" - the install directory of QGIS for my current machine. ...
11
you can use python for getting EXIF info:
from PIL import Image
from PIL.ExifTags import TAGS
from pprint import pprint
def getexif(im):
res = {}
try:
img = Image.open(im)
info = img._getexif()
for tag, val in info.items():
dec = TAGS.get(tag, tag)
res[dec] = val
except IOError:
print im
...
2
My first thoughts were also about lyr.getExtent() but what if you do it this way? Get the source of the layer and describe it.
lyrdesc = arcpy.Describe(lyr)
source = lyrdesc.catalogPath
shpdesc = arcpy.Describe(source)
df.extent = shpdesc.extent
If neither of these work there may be an issue with the extent of the new shp you are adding.
Edit: Actually ...
1
The proper way to set the -q flag in a GRASS-Python script inside a grass.run_command() is quiet = True. In the example given in the question (for grass64), that would be:
grass.run_command("r.stats", flags='l', input=<file>, output=<*.csv file>, fs=',', quiet = True)
Read also the GRASS-Wiki page GRASS_Python_Scripting_Library for examples.
1
I am still unclear on precisely what you are trying to do, but I think you will need to use ArcGIS Add-Ins for Python which only became available at ArcGIS 10.1.
1
If you had a shapefile contained in a folder (note a *.gdb), then Split Layer By Attributes would do the job. In any event, the code contained the scripts associated with the toolbox should provide sufficient background as to how accomplish this task.
1
I will assume you have many records per state that need to be geocoded, in which case you do not want to do this one record at a time, but in batches.
you can tweak your code a little to do this.
import arcpy
states = []
fc ="C:\GISProjects\ALLMembers\LEE_AllMembers.gdb\LeeMembers_Sort"
field = "State"
myList=set([row.getValue(field) for row in ...
1
Very easy, you can use the WKT or WKB formats see Geometry_Outputs:
in PostGIS:
WKT: ST_AsText without CRS or ST_AsEWKT with CRS
WKB:ST_AsEWKB with CRS
in QGIS , see Geometry Construction:
QgsGeometry.fromWkt(wkt)
or WKB
2
If you have not done so yet read up on route hatching. Then you should be able to author a layer file with appropriate route hatching that you can use to add/update a layer using ArcPy.
1
Look at the NME project
With it, you can have all metadata supported by GDAL referenced in an XML (also an SQL output options) using the python script gdalogr_catalogue.py with a command like below
python gdalogr_catalogue.py -d /home/my_directory
After, it's up to you to adapt the retrieved XML from this utility to be able to compare it with anzlic ...
2
Yes you can. Layer actions can call Python so you can do something like this:
from PyQt4 import uic
import os
uifile = r"{path to your UI file}\youruifile.ui"
uiinstance = uic.loadUi(uifile)
uiinstance.exec_()
You could even put that in a python file and call it rather then storing the code in the projet.
In yourcode.py
from PyQt4 import uic
import ...
0
I originally created the routes feature class without any geometry since it was not going to play any part in the network. For this reason arcpy.AddLocation_na could not locate the table where the attribute values were stored.
Once i created a new route from the streets line feature class that i had and included only the number of rows i needed (each row ...
1
You are running 10.1; have you installed any service packs? There is a list of bugs found for "tabulate area crashes" from esri's support page. Check to see if any apply to you, or check other searches about the tool.
0
Just to provide some further color on the solution to this. The first answer worked well to a certain point, but I still had problems with memory leak from ArcGIS.
The final working solution was to split my code into two scripts, where the first code split the files into chunks of 1000 and then fed this into ArcGIS via a subprocess.call-process. This way ...
1
It looks like you are attempting two cursors on the same feature class at the same time. I am fairly certain that will not work.
From Accessing data using cursors, 10.0:
Update and insert cursors cannot be created for a table or feature class if an exclusive lock exists for that dataset. The UpdateCursor or InsertCursor functions fail because of an ...
0
In the web only your solution works!!!! ;)
On Windows 7 & QGIS standalone (no installer OSGeo4W) my .bat is:
SET OSGEO4W_ROOT=C:\PROGRA~1\QUANTU~1
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%PYTHONPATH%
Set PATH=%OSGEO4W_ROOT%\apps\qgis\bin;%PATH%
set QGISHOME=%OSGEO4W_ROOT%\apps\qgis
echo off
call %OSGEO4W_ROOT%\bin\o4w_env.bat
call ...
0
I think the problem is knowing the path python looks in to find scripts. Within FME this can be confusing. A solutions for you might be to add the path to your script in the code like this:
import sys
sys.path.append("c:\\MyPath\")
import module2.py
You can see the python path (where it is looking for scripts) by adding this to your script
import sys
...
1
If you have a string that contains Python code you want to run, you can use the Python Built-In Function exec.
In your case what you could do is:
def CalculateFieldsBasedOnCondition(sourceTable, fieldNames_Condition):
'''
fieldNames_Condition = {'TestID': 'id = "LOC-" + str(row.getValue("OBJECTID"))',
'UPADATED':'Yes'}
usage: ...
1
GDAL doesn't have a ready-made utility for this that I know of. You would probably have to script something.
However, for filling sinks, GRASS has the r.fill.dir and r.terraflow functions. SAGA and some other FOSS raster GIS packages also have sink filling functions too.
Spikes pose a different problem and you'll need to consider a median filter and/or ...
2
You should be able to just add this line after the code you have in your question.
df.extent = lyr.getSelectedExtent()
If that does not work let me know and I will set up a test.
1
QGIS uses gdal_translate utility to do the clipping. It is trivial to run gdal_translate from a shell script over files in a directory to clip all the rasters.
See this post for an example
http://linfiniti.com/2010/11/batch-clipping-with-gdal-and-bash/
To make things easier, you can use the Raster -> Extraction -> Clipper tool to select your parameters. ...
6
Just to add on to the other answers, you can avoid having to explicitly call the file object's close() method if you use a with statement (which automatically calls close() even if an exception is raised):
with open(filename, 'w') as f:
f.write("open=" + jpg + " set to 100")
Also don't forget to write a line-ending if that is the intention...
3
Because you're asking in GIS, I'm including arcpy in my answer. ;)
import arcpy
from arcpy import env
env.workspace = r'C:\your\path\here'
jpgList = arcpy.ListFiles("*.jpg")
for jpg in jpgList:
jpgnum = jpg[:-4]
filename = jpgnum + '.fsv'
f = open(filename, 'w')
f.write("open=" + jpg + " set to 100")
f.close()
And that's ...
7
This is one way to go about it:
This python snippet will loop over a directory, get all of the filenames that end with a .jpg extension and the create a new text file with the new .fsv extension - this file will contain whatever text you choose to write to it.
import glob
import os
os.chdir("\dir") #the directory containing your .jpegs
for file in ...
2
First off, I don't think you want to be saving these as layer files (.lyr). A layer file is only a pointer to data. You need to save the data to a feature class or shapefile. The output from the Make XY event layer tool is a "in memory layer" and it is gone once the session is over. That needs to be converted to feature class to save it to your computer.
...
1
I think you could go either way with this code. Your first code block would run correctly.
I imagine there are some speed differences between looping through rows and doing a calculate values, so the latter is probably much faster.
Try this for the field calculator expression:
expression = "!Field_B! * constant"
arcpy.CalculateField_management(InTable, ...
Top 50 recent answers are included
