New answers tagged

0

you could create a buffer around your point, then intersect the buffer geometry with the line. Alternatively you could use linear referencing tools if you need a distance along the line (which could differ from distance to to point if your line is not a straight line near your point.


2

I know this is an old post, but it showed up on Google for something I was looking for, so I thought I'd post my solution. I don't see this as a 2D curve fitting exercise, but rather a 3D one. By considering the data as 3D we can ensure that the curves never cross one another, and can use information from other contours to improve our estimate for the ...


1

Perhaps you could use something like the following (untested) which uses the glob function to search for rasters in a given directory (you may need to import glob): def function(): myarray = np.array(src_ds.GetRasterBand(1).ReadAsArray()) rows = src_ds.RasterYSize cols = src_ds.RasterXSize BufferDistance = 15 IndexingSquare = ...


1

Use band.XSize instead of width e.g. fmt = "<" + ("h" * band.XSize) In the ReadRaster use row and column numbers, not the real coordinates e.g. for y in range(band.YSize): scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType) Check the data type with gdalinfo if it is really 2byte integer.


2

List pj_list gives you at least projections the proj command understands. The EPSG context in the proj.4 biotop is a simple file based "DATABASE" that connects ONLY EPSG number entries (KEYS) with initialization params of the tool proj. You will find the data in the file /usr/share/proj/epsg if you use proj.4 and one of it's bindings in LINUX. Each ...


0

Finally I found that , space in input feature class directory creates problem in topoToRaster tool execution. from this link :https://geonet.esri.com/thread/51293 ... and the "_inputFeature" is like this: _inputFeature = "D:\Temp\line.shp Max_depth CONTOUR"// "shape_file_path..space..Field..space..TYPE" thanks all....


0

Use the existing arcrest python library. Recommend using version 3.0.1 The sample for deleting rows is here: https://github.com/Esri/ArcREST/blob/master/samples/delete_rows_from_service.py


1

You are using parser() function to bring up a your GRASS module interface, but you are also setting up environment to run GRASS modules without starting GRASS explicitly. Although it is possible to use both, usually you do only one depending on what you want to achieve. In case you want to write a script which would run in GRASS environment as a module ...


0

I was able to use the randomForest package to answer my question which I explain in more detail here: http://stackoverflow.com/questions/34864182/predict-estimate-values-using-randomforest-in-r


3

You're close! The setGroupVisible function requires (int, boolean) as parameters (i.e. the integer you provided for your group as 0). Try the following: root = QgsProject.instance().layerTreeRoot() group1 = root.insertGroup(0, "PGB") legend = iface.legendInterface() legend.setGroupVisible(0, False)


-1

In case you need gdal+python on windows x64, you should compile it manually.


2

This is what fiona module would do, but in a sense you are then adding some GIS functionnalities to Python. Further GIS-like analysis on the geometries can then be performed using shapely.


0

The documentation says "Calling a Basemap class instance with the arguments lon, lat will convert lon/lat (in degrees) to x/y map projection coordinates (in meters)". So you are passing in points around 39.789808E, 3.3395N - those are definitely on land (in Kenya). If you were hoping for points near the coast in Spain, you've just reversed the coordinates. ...


1

The reason why it doesn't work is because the arguments for is_land(xin, yin) are referencing the grid, and aren't latitudinal and longitudinal coordinates, they're simply the index values. Since is_land(xin, yin) takes index values and not lat/long coordinates, you need to find out the index values on the map for the coordinates you're interested in. map ...


2

This is only half an answer as the following code can be used to import waypoints, routes and tracks but not route_points or track_points (these seem to be replaced by the tracks layer). import os path = "path/to/gpx/folder" names = ["waypoint", "route", "track", "route_point", "track_point"] for dirpath, subdirs, files in os.walk(path): for f in ...


1

This page has a good example of doing rule based renderer via python: Rule-based renderer The rule-based renderer styles features using a list of rules, defined by expressions. The rules are processed sequentially, with the first matching rule used to render a feature. The example below defines a new rule-based renderer for different ...


1

The custom function can be used to run a specific function on each feature in a layer. But this will be feature specific. For more uses for this please see the following posts specifically for this: http://www.lutraconsulting.co.uk/blog/2015/06/05/qgis-function-editor/ http://nathanw.net/2012/11/10/user-defined-expression-functions-for-qgis/ ...


1

SampleServer3 = ArcGIS 10.0. SampleServer6 = 10.3. As suggested, its the server which is having problems, not your code. The following works against 10.3 - note I modified your code to use built-ins instead of Requests. (Remember, Server underwent massive changes between 10.0 and 10.1. If you're stuck at 10.0 I'd suggest looking to upgrade as soon as ...


1

I would use a data access module da Update Cursor and the replace() method to accomplish this. The following is the preferred (new) method of using an Update Cursor. import arcpy shp = r'C:\path\to\shapefile.shp' with arcpy.da.UpdateCursor(shp, "your_field") as cursor: for row in cursor: row[0].replace("-", "") cursor.updateRow(row)


2

wait, it looks like it's expecting a file or a file-like object. when you say you have an json embedded in python, it's the whole json text, not the file path? a = "/home/george/boundaries.json" # file-path b = "{ 'type': 'FeatureCollection', 'crs': { 'type': 'name', 'properties': { 'name': 'EPSG:3857' } }, ...


1

As @polygeo suggests, POST requests do have a size limit, but it depends on server and browser. http://stackoverflow.com/a/2880777/3182496


0

Along a different line of thinking: import os ComputerName =os.getenv('COMPUTERNAME') txtFile.write("= Run By: "+str(ComputerName)),txtFile.write ('\n')


3

I suggest use parameter passing to the script while you call this script using task scheduler and grab this parameter using sys module- as below try: data = "= For ArcGIS 10.3.1: Date: "+str(timed)+ sys.argv[1]+'\n' #grab parameters passed except: data = "= For ArcGIS 10.3.1: Date: "+str(timed)+ "Arcgis tool used"+'\n' # else arcgis tool ...


2

You can use mapnik.save_map() import mapnik m = mapnik.Map(1440,800) # ... add styling and layers ... mapnik.save_map(m, "output.xml")


0

Would go for GDALwarp as well. Be sure to be consistent with "postings" and "cell" interpretations of rasters. http://www.remotesensing.org/geotiff/spec/geotiff2.5.html


1

A solution was given by Tim Sutton in A python layer action to open a wikipedia page in QGIS (in one single line): In my case, the field is "URL" and the complete path of one feature is for example http:/..../descriptions/4280614.htm (URL = "4280614.htm") The Python action is then (in one single line): from PyQt4.QtCore import QUrl; from ...


2

You could automate some tasks using the QGIS modeler which allows you to build and run models on selected layers like the following: There's a nice tutorial in the following link which goes into quite great detail on how to create a complex model: Automating Complex Workflows using Processing Modeler. You could also use the Python Console to call ...


0

Based on this new information, I would suggest the following setup to find the street: SELECT point.geom As pointgeom, COALESCE( (SELECT streets.geom FROM streets WHERE ST_DWihtin(points.geom, streets.geom, 30) ORDER BY streets.adt DESC LIMIT 1), (SELECT streets.geom FROM streets WHERE ST_DWihtin(points.geom, streets.geom, 50) ORDER BY streets.adt ...


2

If you look at QgsZonalStatistics.Mean it's part of an enumeration - it's actually the integer 4. You're calling it as if it's a function. Try removing the brackets:- zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath, 'pre-', 1, QgsZonalStatistics.Mean) zoneStat.calculateStatistics(None) Note that these values are done as powers of two; to get ...


4

QGIS has its own internal handling of layers and features. Different data "providers" are used to pull layers in from a variety of sources (eg PostGIS, OGR file types, WFS servers, etc) and translate them into QGIS' own internal formats. Once you have a project fully loaded, including a .shp file plus a .dbf, plus another CSV joined, are all the data in ...


2

You initiate newGeom as GEOMETRYCOLLECTION, but in the loop you combine geometries. combine is the geometric union operation and as such has a much greater job than simply adding a geometry to a collection. It must check topological relationship on a low level, and if objects are disjoint it aggregates them to MultiPolygon geometries. When I interpret your ...


1

The code looks fine, but I'd expect the performance to get worse over time. Each time round the inner loop, you're combining the current feature's geometry with the union of all the previous matching features' geometries. At the start you'll have one simple geometry. Each time round the loop, the geometry gets more complex, so each time round the loop it ...


3

Jibran, it seems that what you are looking for is a nearest neighbour selection. Since it wouldn't matter if the street is 30 or 50 meters away, you just want the closest and then snap to that one, you can do the whole thing from within postgis with the help of just two functionalities: Nearest neighbour selection Point on Line interpolation


5

As many users stated this question has no distinct answer. The results of an analysis may incomplete or fetch the wrong points, dependent on the definition of what should be considered as a corner. The simplest possibility would be the comparison of angles between adjacent points. As this question is familiar with this I will pick up this code fragment. ...


5

The question is about Shapely and Fiona in pure Python without QGIS ("using command line and/or shapely/fiona"). A solution is from shapely import shape, mapping import fiona # schema of the new shapefile schema = {'geometry': 'Polygon','properties': {'area': 'float:13.3','id_populat': 'int','id_crime': 'int'}} # creation of the new shapefile with the ...


0

This is due to a document error - you can't use lat/lon, you have to use llcrnrx, llcrnry, etc. –


3

You can do that in QGIS, without 'shapely' and 'fiona', by using PyQGIS. For a similar arrangement of shapefiles (see next image) from the answer in your link: How to calculate the size of a particular area below a buffer in QGIS This code: mapcanvas = iface.mapCanvas() layers = mapcanvas.layers() feats0 = [feat for feat in layers[0].getFeatures()] ...


2

I solve the similar problem this way: from qgis.core import * from qgis.gui import * layer = QgsVectorLayer(self.shpFilePath, "Track", "ogr") crs = layer.crs() crs.createFromId(4326) layer.setCrs(crs) QgsMapLayerRegistry.instance().addMapLayer(layer)


1

There are two issues here: You need to declare your matches variable before you use it, if you are going to use in the manner that you are. You have to get the unqualified names in order to compare them. Otherwise the name of the user or database etc. will make your names not match. I'm also not sure if you have your workspaces written correctly. I ...


5

I'll try to explain this the best I can so hopefully you can move forward. The short answer is, if your PYTHONPATH is pointing at an actual version of Python, when a product loads up and needs Python it's going to see this directory. In the case of ArcGIS Pro (64bit, using Python 3.4) and ArcMap (32bit, using Python 2.7)... well the use of PYTHONPATH and ...


3

Add an import arcpy statement as the first line of your script. EDIT: Because you can't import arcpy, you have an installation error.


5

There are 2 main reasons why a two raster with the same extent and the same pixel size have different size : 1) the pixel depth 2) the compression in your case, dividing by 2 might produce a float while your input is integer (e.g. unsigned 8bit). Furthermore, the input could be compressed while the output is not. For full control on the output type, you ...


4

Saving this in a Python file and using it in the --code arg should do what you need: from PyQt4.QtCore import * from PyQt4.QtGui import * import os from qgis.core import * from qgis.gui import * from qgis.utils import iface for comp in iface.activeComposers(): print comp folder = r"C:\temp" title = "testing" printer = QPrinter() ...


0

The "optimal" size should be determined by experimenting on your machine, as different machines have different computing power. I've encountered a similar situation before. I recommend you make your chunks even smaller (eg 400*400*644). If it works, then you can consider increasing the size in the subsequent tests.


0

I managed to get it working by updating the .spec file to copy the missing qgis plugin folder into the .exe bundle. I also had to add a check in the code to see if i'm running bundled so I can choose the correct path to the bundled qgis_plugins folder. I assume a better long term fix would be to make this into a runtime hook but for now this is sufficient.


1

So you just want to mosaic all the tiles for a given day? That's a perfect Job for GDALs VRTs. gdalbuildvrt mosaic_049.vrt 049*.hdf or from Python for all days import subprocess import glob for day in range(0, 365): day = str('%0.3d' % day) cmd = ["gdalbuildvrt", "mosaic_"+day+".vrt", glob.glob(day+"*.hdf"] subprocess.call(cmd) edit: ...


2

Solved using Programmatic raster-vector calculation I had some troubles using directly gdal.RasterizeLayer() on layers like propose here in the cookbook but it seems that using "MEM" data source and finally writing it on the disk is maybe better. My solution (which is one of many possible solutions) def rasterizer(shapePath, rasterPath, attribute, ...


0

I couldn't get either to work fully. In the first example, in a directory with both geodatabases and shapefiles, I only got a listing of the feature classes in the geodatabase, but when I commented out the geodatabases portion of the script, I got a list of shapefiles. In the second example, the geodatabases portion didn't work at all, so I copied in the ...


1

Tidiest approach that I have seen is to use python to trigger the pgsql2shp command line utility which is part of the postgis package. Not very pythonic, but very simple and robust. e.g. See How can I get a shapefile from a postgis query? Once you have a shapefile it is straight forward to convert to most common formats using ogr2ogr, also over the ...


1

You can also use triple-quoted strings, this means you can leave your double quotes and single quotes (apostrophes) un-escaped. Can make things a little tricky to maintain, so if you can, leave one space padding at the start and end, like so expr = QgsExpression(""" "my_field" like '%something%' """)



Top 50 recent answers are included