Tag Info

New answers tagged

2

Geopandas plot takes a colormap parameter GeoDataFrame.plot(column=None, colormap=None, alpha=0.5, categorical=False, legend=False, axes=None) The pandas documentation mentions colormaps such as colormap='cubehelix' colormap='Greens' colormap='winter'


0

Well, thank you so much @Aaron. The problem has solved. I have realized that the above code needs to be slightly modified for the case of Google Earth imagery regarding its conversion to grey-scale intensity image. Following is the output of my modified code. Impervious surfaces of building rooftops have been extracted. Nevertheless, the results can be much ...


1

You can simply use: my_sr = arcpy.SpatialReference(my_wkid) Then use my_sr however you wish. ESRI Documentation - ArcPy SpatialReference


1

You can use the OpenCV package in Python for image thresholding. This example shows not only how to perform the binary image thresholding, but also the limitations of this method. Here, I use a 1m spatial resolution NAIP image that shows a dirt road surrounded by arid vegetation. You can see that the road is extracted but there is also a significant amount ...


0

Just for fun, I noticed ESRI 10.3 only shows NULL as 0 so you have no idea what is in there. In QuantumGIS Brighton or later, the field shows as NULL, but you don't know which one. I modified this as the following (I added a zero after copying off Jamie): def RemoveNULL(x): if x is None: return '0' elif x == '': return '0' else: ...


1

Here is an Update Cursor approach: import arcpy shp = r'C:\path\to\your\shapefile.shp' with arcpy.da.UpdateCursor(shp, ["field1", "field2"]) as cursor: for row in cursor: row[1] = row[0][:4] cursor.updateRow(row)


4

kind of figured it out as soon as I posted. Figured I should post the answer... arcpy.CalculateField_management(myFeatureLayer, field2, "!field1![0:4]", "PYTHON") suppose I should add that one can apply any standard Python string functions to !field1! and that this syntax can be used in the Field Calculator if one toggles the Python radio button to the on ...


1

Indeed, you could create a python addin to do this for you! http://resources.arcgis.com/en/help/main/10.1/index.html#//014p00000025000000 Basically, you would be able to create a tool that could sit in your bosses tool bar. Using the tool your boss could click the map and get the data back for the area as you described. The tool would take care of all the ...


2

This works for me and I have labeling only with 2nd word: So maybe you have few atributres without " ", and you get "out of index". def label([SSA_NAME]): x = [SSA_NAME] x1 = x.split(" ") if len(x1)>1: return x1[0] + " " + x1[1] else: return x1[0]


4

Here's a script I put together. The script uses your point and your line feature class, and outputs a new line feature class of line features as desired. The basic steps: Iterate through points Create a buffer around each point Create East-West line from each point. This will be used to slice your buffer in half Create point north of each input point. This ...


0

With the OGR's function UnionCascaded() from osgeo import ogr shp= "yourfile" shapef = ogr.Open(shpor) lyr = shapef.GetLayer() unionc = ogr.Geometry(ogr.wkbMultiPolygon) for feat in xrange(lyr.GetFeatureCount()): fit= lyr.GetFeature(feat) geom= fit.GetGeometryRef() unionc.AddGeometry(geom) union= unionc.UnionCascaded()


1

If you want this as a python script, then you could do something with an update cursor... inputData = "" #Needs to be the full path codeField = "" updateField = "" with arcpy.da.UpdateCursor(inputData, [codeField, updateField]) as updateRows: for row in updateRows: if row[0] == "TestValue1": row[1] = "NewValue" elif #... etc ...


3

I hacked together a solution for this and wrote a blog article a while back on a very similar topic, which I will summarize here. The script is intended to extract a river from a 4-band NAIP image using an image segmentation and classification approach. Convert image to a numpy array Perform a quick shift segmentation (Image 2) Convert segments to raster ...


2

You could look at clustering in scikit-learn. You will need to read the data into numpy arrays (I'd suggest rasterio) and from there you can manipulate the data so that each band is a variable for classification. For example, assuming you have the three bands read into python as red, green, and blue numpy arrays: import numpy as np import sklearn.cluster ...


0

The difference is most likely due to outdated statistics in your GeoTiff. band.GetMaximum() will read the value from the precomputed statistics, which are included in the GeoTiff. If you want to update the statistics you can either use gdalinfo -stats or do it with the Python bindings, provided you have write access to the file. from osgeo import gdal ds = ...


1

The solution is to use from PyQt4.QtGui import QRegExpValidator and from PyQt4.QtCore import QRegExp then use the code: rx = QRegExp('^[1-9]\d{3}$') validator = QRegExpValidator(rx) self.boundwidgets['flatenr'].widget.setValidator(validator) See https://github.com/DMS-Aus/Roam/issues/239 for more details about the solution.


2

Several things to check here: Firstly, make sure that %QGIS_PATH%\apps\qgis\python\plugins\processing is added to the PYTHON_PATH environment variable. Secondly, instead of: import processing from processing import * Use: import processing from processing.core.Processing import Processing Finally, before calling any algorithms using processing, do ...


1

When you created your script tool you would have gone through that wizard to wire up your script to a tool interface. At the stage where you point it to the script python file (.py) there is a check box which is usually un-ticked which is Show command window when executing script, try ticking that on?


0

You set a workspace, so any outputs from the script will default to that location unless you explicitly put them somewhere else. Create a folder for the outputs, then write your output rasters to that folder. import arcpy import os workspace = "D:\Jiawei default download" feature_classes = [] outputFolder = r'c:\outputfolder' walk = ...


3

It's just a widget so you sure can. Normally you would add it to a layout however that will resize your canvas when it pops up which is a bit annoying IMO. What you can do is just some logic to resize the bar to match the parent size when it is shown. Here is some code that I use to do just that: class MessageBar(QgsMessageBar): def __init__(self, ...


2

There's several methods but I will mention the ones I normally use, both using the Python Console. For SINGLE shapefiles: Method 1 - Using a method utilising the QGIS algorithm - Delete Column: output_0=processing.runalg('qgis:deletecolumn', myshapefile.shp, "myFirstField", None) output_1=processing.runalg('qgis:deletecolumn', output_0['OUTPUT'], ...


0

Not directly related to the question, but since I landed here I'm guessing others might as well. I was also having problems with numpy on Ubuntu 14.04 (when using GRASS 7.0.0). The issue in my case was the numpy version, which was 1.9 whereas the subpackage numpy.oldnumeric (needed by GRASS) is only shipped with numpy until version 1.8. The solution was ...


1

You should specify more your question. It's unclear what mean by doing a mean and stdev. Here's a sample code you can run in QGIS Python console (still without further iteration): """ RPG - Random Point Generator Version: 2015-06-30 Geometa Lab HSR, http://twitter.com/geometalab License: MIT license, http://opensource.org/licenses/MIT To run this: 0. ...


0

Just solved it. Uninstalled and intalled SciPy module (http://www.kyngchaos.com/software/python) again and it worked.


2

You're script above probably still has a lock onto the data. Specifically this line: rows = arcpy.SearchCursor(path, fields="alpha; beta") You've got two ways to get around the lock. Delete the rows variable after your loop rows = arcpy.SearchCursor(path, fields="alpha; beta") for row in rows: print("{0} => {1}").format(row.getValue("alpha"), ...


1

I was able to figure this out myself using the Shapely library, like so: from shapely.geometry import LineString, Point, Polygon epsilon = 0.001 def is_between(a, b, c): crossproduct = (c.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y) if abs(crossproduct) > epsilon : return False # (or != 0 if using integers) dotproduct = (c.x - a.x) * ...


2

You can activate Python plugin "ftools" and do "Random Points..." (from menu Vector > Research Tools).


1

There is a much simpler way using a Dictionary object. To be accessible/updateable for each iteration it must be outside the scope of the function (global) otherwise it's reset for each calculation. preDict = {} def CalcField3(Field1,Field2): global preDict # tell python this isn't a new variable but the other one if preDict.has_key(Field2): # if ...


2

You are on the right track. The thing that bit you are the backslashes in your prefix path \ which are escape sequences in python. Just do it the same way as you do it with the layer path and use slashes / instead. qgishome = "C:/OSGeo4W64/apps/qgis" A useful snippet that you can run after QgsApplication.initQgis() that makes sure that the providers are ...


2

The folder for the PYTHONPATH must be the compiled one and not the source one. Your build path should be something like /home/user/QGIS-master/<builddir>/output/. I am guessing you are running from build your QGIS version so you just export the environment var like: export PYTHONPATH=/home/user/QGIS-master/<builddir>/output/python/ maybe it ...


1

The "browse" feature may have gotten "stuck" on an empty or illegal file path and name, perhaps involving a directory in the path. Rather than use the browse button, try typing in a fully qualified file name, including the path. Start at the "beginning" with the top level device or "/" directory. If all works well, as you type you should receive file ...


1

Assuming "Nearby:" is at the beginning of the line (and never in the middle), you can simplify this logic with a one-liner, without the need of a codeblock: str(!Street1!).replace('Nearby:', '').lstrip()


2

as written, the function returns a string with 'Nearby:' removed only if the 'Nearby:' test is true (if 'nearby:' exists within the field string. If the test is false, the function is set to 'pass', which returns no value (deleting the original value in the process) to remedy, a value should be passed for both cases - true and false - something like def ...


0

You need to install cx_oracle module so python will be able to import it. Use this https://pypi.python.org/pypi/cx_Oracle After you install it then open python console and try import cx_Oracle. Make sure you download the cx_Oracle to match you python version. By the way what is your OS?


1

assuming that the goal is to determine an elevation with a user-supplied set of coordinates, there should be a way to do this with gdal - i use a function like this, which came almost verbatim from this answer (assumes raster is not rotated). rast = 'myIMGfile.img' mx = XCOORDINATE my = YCOORDINATE def readRastPix(rast,mx,my): src_ds=gdal.Open(rast) ...


1

Since the given Spatial References in ArcGIS have got their own parameters for the vertical part, as far as I know there is no way of introducing vertical part by using the given factory codes or like sr = arcpy.SpatialReference(horizontal_wkid,vertical_wkid) that one may think it should work, but it is not. However, you can create your own spatial ...


0

For a pure python solution, you can use os.path.isfile to check if the shapefile exists and is indeed a file. import os myfile = r'C:\temp\myshp.shp' if os.path.isfile(myfile): print "The file exists" # Do something


1

As a cavaet, be very careful about doing any geometry comparisons particularly geometry objects. One should be aware that the geometric operations are 2D whether you have 3D objects. You will have to compare the individual values to check for equality, preferably in Python. On a second, but related point, be careful what the geometry means. Consider the ...


2

Use the ModelBuilders' Calculate Field to insert a python code like this: Expression: calculate(!Field1!, !Field2!) And a code block of: previousValue1 = "" previousValue2 = "" previousUnchangedValue1 = "" def calculate(field1, field2): global previousValue1 global previousValue2 global previousUnchangedValue1 #First item in list if not ...


1

If we print out your formatted string we get this: '"OB_ID" =\'id_LCC\'' This is not what you want. Format considerations for file GDBS are as follows: Field names are not qouted If the value is a string, use single quotes If the value is a number, do not use quotes Review the help docs for more details Use .format to simplify the variable substitution ...


0

@Vince is right, this would be trivial in GIS software. But if I can assume you are always working in the same area and you have access to GIS to get started: A raster data file is essentially an array of data with a few lines of header on top. If you build a "land mask" you will end up with an array of presence/absence of land (or water, depending on which ...


1

I find a way of overlay vector files (if there is a vector file already open) and a raster file. orto = '...image.tif' fileCV = '.../layer1.shp' layers = [] # display file in canvas fileInfo = QFileInfo(fileCV) layer = QgsVectorLayer(fileCV, fileInfo.baseName(), "ogr") # raster layer rfileInfo = QFileInfo(orto) rlayer = QgsRasterLayer(orto, ...


2

areas = [] for line_feature in line_layer.getFeatures(): cands = area_layer.getFeatures(QgsFeatureRequest().setFilterRect(line_feature.geometry().boundingBox())) for area_feature in cands: if line_feature.geometry().intersects(area_feature.geometry()): areas.append(area_feature.id()) area_layer.select(areas)


2

Yes, it is possible. ArcPad comes with a toolbox and geoprocessing tools. There are two tools that might be of interest, the ArcPad Checkout and the ArcPad check in. If you export them to py you will learn their arcpy parameters and variables. arcpy.ArcPadCheckout_ArcPad("", "", "", "", "", "", "", The_full_path_and_name_of_the_AXF)


3

Use the geometry object's equals method. With polygons and polylines, if their symmetric difference is empty they are considered equal. So for example, in the below, g2 is the same as g1, it just has an extra vertex in the middle, and g3 is opposite to g1, so always equal. >>> g1 = arcpy.Polyline(arcpy.Array([arcpy.Point(1, 1), arcpy.Point(3, ...


0

Supposing that geom is a shapely geometry object, you can use from shapely.geometry import * import json geojson = json.dumps(mapping(geom))


4

You can use Find Identical, using the geometry of the features. Alternately, you can spatially join a feature class to itself, and use the "ARE_IDENTICAL_TO" as your match option. This can be useful if you have a bunch of features stacked on top of each other and you want to deduplicate (features with a Join_Count > 1). I recently used this to clean up some ...


3

Or you could use a describe and .baseName fcs = arcpy.ListFeatureClasses("","","") for fc in fcs: desc = arcpy.Describe(fc) print fc print desc.baseName ActualName = desc.baseName


2

This is really a Python rather than ArcPy question, so would be best researched at Stack Overflow. However, you can try changing: arcpy.CalculateField_management(fc, "Policy_Ref", "'" + fc + "'") to: arcpy.CalculateField_management(fc, "Policy_Ref", "'" + fc.replace(".shp","") + "'")


0

(answering own question, hoisted from comment to question) It seems that I may have reversed lat,lng. Based on transforming (180,0), it seems that this epsg:3857 coordinate system is centered at the equator and scales by meters at the equator, using a circumference of the earth of 40075016.68557849m. So I guess that's how to scale the result.



Top 50 recent answers are included