qgis - How to access the road graph plugin in python?

I am developing a few scripts in Python for QGIS 1.8 to automate some geo-tasks and one of them would be to calculate the shortest distance from a hospital to each building in a town.

I would like to take advantage of the already existing Qgis plugin but, as it is a C++ plugin, I am not sure whether it is possible to access it from Python? Is there a binding or a trick for that?

If not, any further advice would be appreciated! As the base data is OpenStreetMap data stored in a PostGIS DB, I may also try with pgRouting. But the QGIS plugin would be my favourite choice if possible.

PS: this question is close to this one: Is there a way to access QGIS plugins in Python? but I believe it's still different enough to be worth a new thred.

Answer

gis-lab.info has a tutorial on using the network analysis library in QGIS using Python. It's not in English but Google Translate should be able to make it comprehensible. Good luck!

arcgis 10.2 - Prevent parameter defaults in service definition file

I have a Python toolbox with a single tool designed for 10.2. The tool is intended to be a geoprocessing service on a server. So I run the tool, and after it succeeds, I right click the result and choose Share As --> Geoprocessing Service. In the wizard, I choose "Save a service definition file" and move through the wizard until it switches to the dialog where I can set service properties. I Analyze, and there are no errors or warnings.

At this point, I notice that my parameters are being given default values. In particular, they get the values that I used when I ran the tool to create the results file. ArcGIS's dialog will not let me modify the default parameters; the text boxes are read only.

I want to prevent this. Ideally, my parameters would have no default values, but empty string default values would also be acceptable. All my input parameters happen to be strings. How can I prevent ArcGIS from using the result's parameters as the defaults?

Answer

As @AlexTereshenkov noted, there doesn't seem to be a way to prevent this.

So, I have come up with the following work around: I have revised my services so that bad inputs will generate empty output. (For me, this was an empty shapefile or geodatabase.) This is far, far from ideal. An error message would be more useful to users, but my particular scripts can live with this workaround since end users will never call them directly. All calls to this service come from a web application (It's not even publicly accessible.), and the web app can ensure the input is correct.

By generating empty results for bad inputs, I can run the tool on a desktop with nonsense inputs and generate an empty data set. Using this run to generate the service definition file, my services now have nonsense defaults and generate empty output if those defaults are used. This lowers the risk of someone getting the wrong output because of a programming error.

As I said, this isn't ideal, but it does do roughly what I need without modifying anything after publishing. I wouldn't recommend this solution in all cases, but I think it's workable in cases where some external service can filter the inputs.

javascript - Setting max zoom in mapbox map?

     mapObject = L.mapbox.map('mapDiv')
     .setView(defaultLocation,defaultZoom)
      .addControl(L.mapbox.geocoderControl('mapbox.places'));


             L.mapbox.styleLayer('mapbox://styles/gauravcatstech/civ82vos7008f2is5x92cp     vye').addTo(mapObject);

          L.control.scale({maxWidth: 100,metric: true, imperial: false}).addTo(mapObject);

How would i set max zoom in mapbox map..i have try object.maxZoome but sitt not able to set maxzoom

qgis - Joining parcel/building attributes based on coordinates

I downloaded a building footprint shapefile for a city and imported it as a vector layer in QGIS. I was hoping to see street address information in the attribute table, but there was none. I subsequently downloaded and imported a shapefile containing the land parcels for the same city (which DO have street addresses), and it overlays nicely such that buildings appear within parcels.

I would like to join the attribute tables such that my building footprint file contains addresses; however, there is no common key on which I can join (I see a Shape_STAr and Shape_STLe but can't figure out what these are).

Is there a way for me to (programmatically via a framework or otherwise) determine which buildings fall inside a parcel based on the coordinates such that I could join the data in the manner I desire?

I've done some Google searches but can't figure out where to start.

Answer

Seeing as though you have no matching attributes to join you will need to perform a spatial join. This is under Vector>Data Management Tools>Join Attributes by Location

geometry - How is ST_PointOnSurface calculated?

The PostGIS documentation states that ST_PointOnSurface returns "a POINT guaranteed to lie on the surface". It seems like this function could be trivially implemented to give results that satisfy the documentation but provide little real-world utility, though I'm certain that PostGIS provides a non-trivial implementation.

This introduction to PostGIS provides a nice comparison and contrast of ST_Centroid with ST_PointOnSurface and says that "[ST_PointOnSurface] is substantially more computationally expensive than the centroid operation".

Is there a more thorough explanation of how ST_PointOnSurface is calculated? I've been using ST_Centroid, but have encountered some edge cases in my data where the centroid is outside the geometry. I believe that ST_PointOnSurface is the correct substitute, but the function name and documentation leave room for uncertainty.

Further, is the computational expense of ST_PointOnSurface incurred even if the centroid does lie within the geometry already?

Answer

Based on a few experiments, I think ST_PointOnSurface() works roughly like this, if the geometry is a polygon:

1. Trace an east-west ray, lying half-way between the northern and southern extents of the polygon.


2. Find the longest segment of the ray that intersects the polygon.


3. Return the point that is half-way along said segment.

That may not make sense, so here's a sketch of a polygon with a ray dividing it into a northern and southern parts:

             _

            / \             <-- northern extent
           /   \
          /     \
         /       \
        /         \      __
       /           \    /  \
      /_ _ _ P _ _ _\  / _ _\  P = point-on-surface
     /               \/      \
    /                         \
   /            C              \   C = centroid

  /                             \
 /                              /
/______________________________/  <-- southern extent

Thus, ST_PointOnSurface() and ST_Centroid() are usually different points, even on convex polygons.

The only reason for the "surface" in the name, I think, is that if the geometry has 3D lines then the result will be simply be one of the vertexes.

I would agree that more explanation (and better naming) would have been useful and hope a GEOS programmer might shed some more light on the matter.

Is it possible to run an R script on a layer in QGIS?

I'd like to run an R script on a layer loaded in QGIS (to rotate a set of points created there using the Vector:Research Tools:Regular Points tool). I thought ManageR might work but cannot get it working on OSX ...

Of course, I could do it manually (create layer -> run script from within R -> find output and reload into QGIS) but surely there is a more elegant solution?

Answer

It is possible to get manageR working on OSX (Lion 10.7.2) at least.

It looks like the problem with manageR / rpy2 was actually problems with the previous R builds (2.13 at least) - symlinks were referring to the wrong location.

1. Update R to 2.14 (binary from r-project.org);


2. Reinstall QGIS (Kyngchaos installer);


3. Reinstall GDAL_complete (Kyngchaos again);


4. Reinstall rpy2 (latter via pip), rebooted ...

... and manageR now loads and works. Haven't tried this in Linux yet.

UPDATE: The SEXTANTE plugin can be used to write scripts that call R, with some modifications in OSX to make it work properly.

arcgis desktop - Removing spatially duplicate features using ModelBuilder?

I am creating a model to remove features which are duplicated within a polygon feature class. Most of the features are duplicated once, but about a quarter have up to three or four exact copies. I do not have ArcInfo, so I cannot use Delete Identical.

My model should contain the following:

1. Make Feature Layer on the input feature class specified by the user to output Output_FL


2. Get Count to count the number of rows within the feature class


3. Create a variable counter and initialise to 0


4. Add a While iterator where counter < row count


5. Select Layer by Attribute on Output_FL to create a NEW_SELECTION where OBJECTID = counter


6. Select Layer by Location on Output_FL where the relationship is ARE_IDENTICAL_TO


7. Select Layer by Attribute on Output_FL to REMOVE_FROM_SELECTION where OBJECTID = counter


8. Delete Features from Output_FL


9. Increment counter

I have finished all the steps except step 9. I am not sure where/how I should increment the counter. I'm also not sure if I should put a Copy Features tool somewhere before the while loop as a backup mechanism since this model will be modifying the input data, and not everyone will make a backup of their data before running the model.

Answer

You should try this, Add XY coordinates to your featureclass, then Dissolve on the X,Y fields. This should remove your stacked, duplicate features. You can then join the attributes from your source file. You may need to concatenate X,Y to form a single field for the join.

Drawing circles on map using OpenLayers 3?

I want to draw circles of certain radius on a map using openlayers . But the latitude,longitude & radius values are given in a excel sheet.

How do I do this work in OpenLayers?

I heard that OpenLayers doesn't take data from excel sheet so to which format I have to change my excel file?

And how to draw circles on the map?

Answer

Circles in Openlayers

http://geographika.co.uk/samples/geodesic/

http://geographika.co.uk/creating-a-geodesic-circle-in-openlayers you will need your data in json - Google Spreadsheets can do this

http://spreadsheets.google.com/feeds/feed/key/worksheet/public/basic?alt=json-in-script&callback=myFunc

http://code.google.com/apis/gdata/samples/spreadsheet_sample.html

spatial analyst - Adding two rasters using map algebra in ArcPy?

I want to simply add two rasters within a python script using map algebra, which from what I understand is the equivalent to adding two rasters in the raster calculator.

If I do this using raster calculator and simply type "Raster1 + Raster2", it works perfectly. However, if I run this in a python script, I get error000732: Input Raster: Dataset Raster1 does not exist or is not supported.

Here is my code:

import arcpy
from arcpy import env
from arcpy.sa import *
raster1 = ("C:/raster1")
raster2 = ("C:/raster2")
raster3 = Raster("raster1") + Raster("raster2")

I set up the syntax based on the Esri help page for Map Algebra http://desktop.arcgis.com/en/arcmap/latest/extensions/spatial-analyst/map-algebra/an-overview-of-the-rules-for-map-algebra.htm

I'm not sure why my rasters are recognized within the raster calculator but not map algebra.

Answer

What @Joseph says. The Raster class needs a string containing the path to your raster. You've correctly assigned this to a variable, so you then need to pass those variables to instantiate your Rasters. This code should work:

import arcpy
from arcpy import env
from arcpy.sa import *
raster1 = ("C:/raster1")
raster2 = ("C:/raster2")

#Note lack of quotes in following line
raster3 = Raster(raster1) + Raster(raster2)


#Don't forget to save eg:
raster3.save("C:/raster3")

The reason the documentation you link shows the rasters being passed to map algebra statements in quotes is because they assume you are typing in the python console in ArcMap. If you're raster datasets are open in ArcMap you can refer to them by their name in the table of contents, using quotes.

In which case the necessary code would be

import arcpy
from arcpy import env
from arcpy.sa import *


#No need to assign paths to variables.
#NB You still need to call `Raster()` in order to tell Arcpy that the '+' is for raster math. 
raster3 = Raster("raster1") + Raster("raster2")

#Don't forget to save eg:
raster3.save("C:/raster3")

Is it possible to put Geoserver behind a proxy and access the WMS/WFS using OpenLayers

I just start to work with Geoserver and Openlayers and one of the main concerns of the project now is do not expose the geoserver web interface outside of the internal network. No access to www.mydomain.com/geoserver/

However the application that is available in the Internet needs to get data (wms/wfs) from the GeoServer. This all through OpenLayers.

Is this possible?

Answer

It is possible, follow the instructions at http://ian01.geog.psu.edu/geoserver_docs/software/java.html but instead of /geoserver use /geoserver/wms and /geoserver/wfs

qgis - Why do area measurements from "measure area" tool and field calculator differ?

My problem is when using the "measure Area" tool on the toolbar I get the correct value but when using the field calculator in the Layer Attributes the area is very different including -ve areas. Both the project and layer are using a UTM CRS. Using version 2.0.1

arcgis 10.1 - Adding raster catalog to MXD using arcpy.mapping?

Using the bit of code below adds the raster catalog footprint polygon layer.

mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
addOrtho2011 = arcpy.mapping.Layer(r"P:\Mapping\Data\Orthos.gdb\Orthos_2011")
arcpy.mapping.AddLayer(df, addOrtho2011, "TOP")

I can't figure out what else I have to do to add the actual ortho part of the raster catalog. Any suggestions?

Answer

Your raster catalog is displaying as a wireframe due to a setting on the layer properties which is described here.

This setting does not appear to be accessible to ArcPy so as a workaround you could change the scale after adding it so that it is under the threshold set.

However, I think a better solution is to check the box on the layer properties so that wireframes do not show at all (which will be unwise if too many) and then Save As layer File so that you can add the layer from a layer file as per the code example here.

How to rotate (svg) symbols in QGIS map composer to match map rotation?

Because of my thesis layout requirements, I have to rotate 270° some of my maps in composer.
There are (many) svg smbols (icons), and they are rotated as well.
I would like to (non-)rotate them in order to keep them upright, so one doesn't have to continuosly rotate the page or the neck ;)
Is there a way to do that? Note that I would like to do it on the composer stage since some of my maps are normally placed north-south. Thanks for the attention!

Answer

Here's an approach which would work in QGIS 2.2 and above. You can set a data defined rotation on the symbol which matches whether the map is being rendered in a specific composer map or to the main map canvas - this is done using the '$map' variable. $map will be 'canvas' for the main window, or a map item's id if it's being rendered in a composition.

So, to begin, set the map item's id to something like "rotated_map_01":

Then, in your layer's symbol properties, click the "Data Defined Properties" button:

Check "Angle" and click the expression builder button. You want to enter an expression along this lines of

CASE WHEN $map LIKE 'rotated_map%' THEN 270 ELSE 0 END

This expression will return 270 if the map item's id begins with "rotated_map_", or 0 otherwise. So, in any composer maps with an id "rotated_map_01", "rotated_map_02", etc the symbol will be rotated 270 degrees!

d3 - Remove outlying territories/islands from Europe map

Using a shapefile from Natural earth, I am trying create a map of Europe with d3-geo’s new command-line interface. The following:

shp2json -n ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp \
| ndjson-filter '["ALB", "AND", "AUT", "BEL", "BGR", "BIH", "BLR", "CHE", "CYP", "CZE", "DEU", "DNK", "ESP", "EST", "FIN", "FRA", "FRO", "GBR", "GIB", "GRC", "HRV", "HUN", "IMN", "IRL", "ISL", "ITA", "LIE", "LTU", "LUX", "LVA", "MCO", "MDA", "MDA", "MKD", "MLT", "MNE", "NLD", "NOR", "POL", "PRT", "ROU", "SMR", "SRB", "SVK", "SVN", "SWE", "UKR", "VAT", "XKX"].includes(d.properties.adm0_a3)' \
| ndjson-reduce \
| ndjson-map '{type: "FeatureCollection", features: d}' \
| geoproject 'd3.geoConicEqualArea().parallels([43, 62]).rotate([0, 0]).fitSize([960, 960], d)' \
| geo2svg -w 960 -h 960 \
> europe.svg

...creates this output:

How can I remove all the overseas territories (e.g. French Guyana) and distant islands, so I'm left with a map of just 'mainland' Europe? Nothing west of Portugal, nothing south of mainland Spain etc. More like this:

Answer

I'd bring it into QGIS and use Vector > Geometry > Multipart to Singlepart.

In Natural Earth, France is a multipolygon so it is one feature with lots of individual polygons (France, French Guiana, Reunion etc).

Doing this breaks each of these polygons into its own feature - so your shapefile jumps in size from ~250 (roughly one feature per 'country') to several thousand features.

Then select the ones you want with the select tool, dragging your selection over the region of Europe you want, and save as a new shapefile (making sure you check "save selected features only")

You should now be able to plug in the new shapefile into your first line of code and it should work.

postgis - How do I get the closest point on a road to a point when I have imported OSM data using osm2pgsql?

I have been trying to use OSM data to find the closest point to a road for a given point. I imported the OSM data into postgis using osm2pgsql, and I am getting results back, but they are at the absolute edge of my imported data (and the point I'm trying to specify is near the middle) and distance reported is about 8000km. I'm pretty sure it has to be something to do with mixing different coordinate systems but I'm quite new to this so not sure.

The SQL query I'm trying to use is:

SELECT osm_id,

   highway,
   name,
   ref, 
   ST_X(ST_ClosestPoint(ST_Transform(r.way,4326),
       ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326))),
   ST_Y(ST_ClosestPoint(ST_Transform(r.way,4326),
      ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326))), 
   ST_Distance_Sphere(ST_ClosestPoint(ST_Transform(r.way,4326),
      ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326)),
      ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326))

FROM planet_osm_roads r 
ORDER BY 7 ASC
LIMIT 10;

Answer

First of all, ST_Distance_Sphere returns in meters, so you are actually looking at 8km, which might be more reasonable. I suspect, also, that you have you lat/lon the wrong way round -- your point is somewhere off the coast of the Seychelles, not Carlisle, and there are not too many roads in the Indian Ocean.

Also, while the query planner will optimize this away, there is no need to repeatedly enter the same point. Instead, create it in a sub-query and reference it, eg:

SELECT 
 osm_id,
 highway,
 name,

 ref, 
 ST_X(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom),      
 ST_Y(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom),
 ST_Distance_Sphere(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom), point.geom)
FROM 
  planet_osm_roads r,
  (Select ST_SetSRID(ST_MakePoint(-3.198706, 55.938627), 4326) as geom) point
ORDER BY 7 ASC
LIMIT 10;

However, this is still very inefficient. To make better use of a indexes, look at using ST_DWithin (r.way, point, distance) ORDER BY ST_Distance(r.way, point). If you look at the docs for ST_DWithin you will get an idea how it is used.

ArcGIS median value from layer?

I've been hunting around ArcGIS 9.3 and having no success at finding a tool/function that can give me the median (not the mean/avg) tabular value of a spatial database. For example, if I have a data layer of mountain peaks, I would like to be able to determine the elevation of #35 of 70 mountain peaks that exist within a county.

Conceptually it should be pretty easy. Count the total records, sort by elevation ascending, get the value from the elevation field for record # (count/2). Would I need to script this out in Python, or am I missing something obvious that would allow me get this more quickly? I find it hard to believe no one else has had a need for this function and created a solution.

To add some complexity to the matter, what if I wanted the median elev value for peaks in every county within the state in a single table? Basically the same output format as Frequency or Summary Statistics tools when using the Case Field option.

---

I was able to adapt the script linked by @VietThanh Le to meet my needs. I have posted the updated script and ArcToolbox at: http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=A0CD1751-1422-2418-882E-001EE0DC0D35

Answer

I found a python script called Calculate Median Value in ArcScripts, but I have not tried it:

This script tool calculates the median value of one entire field and posts that single median value to every row in one other entire field.

add new SRID to sql server

Is it possible to add a new SRID to sql server 2012?

I've tried

insert into sys.spatial_reference_systems values (4938, 'GEOCCS["GDA94",DATUM["Geocentric Datum of Australia 1994",SPHEROID["GRS 1980",6378137.0,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0.0,0.0,0.0,0.0,0.0,0.0,0.0],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["m",1.0],AXIS["Geocentric X",OTHER],AXIS["Geocentric Y",EAST],AXIS["Geocentric Z",NORTH],AUTHORITY["EPSG","4938"]]', 'metre', 1)

But I get

Ad hoc updates to system catalogs are not allowed.

Haven't had much luck on google with any answers as well

Answer

Not answer per se ,but SRID in MS SQL server is allmoust useless because it dosen't support transformations. Also i think you can use SRID which isn't in spatial_reference_systems table for geometry. see : Thread about it

qgis - How to find the polygon inside which a point lies?

I have a layer with polygonal features. Each feature has attributes and values. I also have a list of co-ordinates and I would like to know which feature (or polygon) the co-ordinates lie in.

Could someone please guide me on how to go about this? Is there a function in the API that can help me achieve my goal or should I use some computational geometry algorithm to do it myself? I know how to do the latter but it would save me some time if there was a built in function already.

Thanks.

Answer

Few features

What you probably want to do is:

1. Create a list of QgsPoint from your coordinates


2. Iterate over all your layer features (polygons)


3. Loop over the list of points (within the iteration)


4. Call feature.geometry().contains( point ) to check if you've got a match

Of course you can now improve performance if you e.g. know, that a point can only be contained by one polygon, you can remove a point from the list, once the appropriate polygon has been found.

Lots of features (Using SpatialIndex)

As pointed out in the comments, a spatial index can be used to speed up the process significantly.

The steps here would be

1. Create a list of QgsPoint from your coordinates


2. Create a QgsSpatialIndex


3. Iterate over all your layer features (polygons)


4. Add the features to your index with insertFeature



5. Iterate over all your points


6. Call index.intersects( QgsRectangle( point, point ) ) to check if you've got a match

There is also a code example by NathanW

Fixing Script Error Invalid Pointer from ArcGIS Geoprocessing Tools?

I am just learning ArcGIS 10.2 and it was working well until I tried to use the Export to a Geodatabase (single) option for a dbs table I made. Since then, it seems I now get the same Script Error for import/export functions. However, I just tried to export to Geodatabase (multiple) and error did not pop up. I spoke with tech support and they could not help me so far. I really do not know much about GIS so it is hard for me to fix this problem and cannot proceed with some things.

latitude longitude - find latlngs that are some x distance away from a point to south to north to east to west

I have a latlng i want to find four points with equal distance from this latlangs in all directions east, west, south an north is there any formula .

I am trying to build polygons around higways so draw polygon around poly line the distance will be less than 1 km.

Answer

Background

For these purposes, the Earth is assumed to have the shape of an oblate ellipsoid of revolution around its axis. The major axis, a, is the radius of the Earth at its equator. The minor axis, b, is the radius of the Earth at either pole. The squared eccentricity (measuring the "squishedness") of this ellipse is e^2 = 1 - (b/a)^2.

For simpler but approximate formulas, we may ignore the eccentricity and take the earth to be a perfect sphere. In the formulas the follow, set e^2 = 0 and a = b = (2a+b)/3, which I will write as R. (The justification for this is found at How accurate is approximating the Earth as a sphere?.)

Formulas

Distances east-west lie along a parallel of latitude, which is a perfect circle. At the geodetic latitude f, the radius of this circle is

r = a * cos(f) / sqrt(1 - e^2 sin(f)^2).

(The spherical approximation is r = R * cos(f).) Therefore, a distance x in the east-west direction corresponds to an angle of

dl = x / r

radians. Convert that to degrees (or grads or whatever you prefer to measure longitude in) and add that value to the longitude (for an east direction) or subtract it (for a west direction).

Laying distances off in the north-south direction requires computing elliptic integrals. However, for the very short distances contemplated in this question (just a few kilometers, or less than 0.001 times a), we may assume the meridian is approximately circular. The best approximating circle ("osculating circle") has radius

s = a * (1 - e^2) / sqrt(1 - e^2 * sin(f)^2)^3.

(The spherical approximation is s = R.)

Therefore, a distance y in the north-south direction corresponds to an angle of

df = y / s

radians. Add that value to the latitude (for a north direction) or subtract it (for a south direction). If the result takes you over a pole, you will need then to subtract it from Pi radians (or 180 degrees or whatever) and add or subtract Pi radians (180 degrees) to the longitude.

Cautionary note

The east-west points are actually closer to the original point than you might think. This is because the desired distance is laid off along a circle of latitude rather than a straight line. For instance, a point 10,000 kilometers to the east of a location at 40 degrees latitude is located only 9100 kilometers away. These differences are inconsequential for short distances (less than about 1000 kilometers).

Data

In the World Geodetic System, the major axis has length a = 6378137 meters, minor axis b = 6356752.3142 meters, and inverse flattening f = 1/298.257223563. This is equivalent to a squared eccentricity e^2 = 0.00669438000426. (These values are connected by e^2 = 1 - (b/a)^2 = 2f - f^2.)

Worked Example

Consider the location at (39.9522° N, 75.1642° W). We will find four points exactly one kilometer away in the east, west, north, and south directions. First compute

sin(f) = sin(39.9522°) = 0.642148
sqrt(1 - e^2 sin(f)^2) = sqrt(1 - 0.00669438000426 * 0.642148^2) = 0.998619

Plugging these into the equations gives

r = 4,896,117.456 meters
s = 6,361,763.203 meters

(These values are over-precise for this application, whose coordinates are given only to six significant figures: the extra precision will help when checking a software implementation.) The offsets are

dl = 1000 / r = 0.0002042434662 radian = 0.0117023 degree
df = 1000 / s = 0.0001571891264 radian = 0.00900627 degree

That places the four points at

(39.9522000, 75.1759023), (39.9522000, 75.1524977)  (east and west)
(39.9612063, 75.1642000), (39.9431937, 75.1642000)  (north and south).

In the spherical approximation, the four points are

(39.9522000, 75.1759316), (39.9522000, 75.1524684)  (east and west)
(39.9611601, 75.1642000), (39.9432399, 75.1642000)  (north and south).

The approximate locations are 2.5 meters too far in the east-west directions and 1.45 meters too close in the north-south directions. This is a typical level of error (it is approximately the size of the flattening, which is around one part per 300, or 3.3 meters out of 1000 meters). When this much error is acceptable, the simpler spherical formulas are fine; otherwise, you need the ellipsoidal formulas.

Limitations

The east-west results for the ellipsoid are completely accurate, but the north-south results depend on approximating an elliptical meridian by a circle. They will degrade in accuracy over longer distances, but will still be pretty good out to a few thousand kilometers. For instance--to continue the example--the point the formula places a thousand kilometers north is actually 1000.79 kilometers away, an error of less than one part in 1200.

arcpy - Why is FeatureClassToFeatureClass tool Failing on Long cell values (Arc 10.2)?

I have written a script to convert a CSV into a Feature Class. The actual line that does the conversion is the very last line in the script in which I call the FeatureClassToFeatureClass tool.

Here is my Code Snippet:

import arcpy
Source_CSV = r"\\int.ec\Data\Stations.csv"
Target_FC = r"Q:\Data\Stations.gdb\Stations"
Description = arcpy.Describe(Target_FC)
Target_FC_name = Description.name
Target_GDB = os.path.split(Description.catalogPath) [0]
Temp_Layer = "Temp_Layer"
SR = arcpy.Describe(Target_FC).spatialReference
arcpy.MakeXYEventLayer_management(Source_CSV, "Lng", "Lat", Temp_Layer, SR, "")

arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, "", "", "")

I am getting the following error:

File "Q:\scripts\CSVtoASMADPointFeatureClass.py", line 10, in 
    arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, "", "" , "")
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\conversion.py", line 1675, in FeatureClassToFeatureClass
    raise e

ExecuteError: ERROR 001156: Failed on input OID 1263820748, could not write value 'LS_ blue with white trim_ has overflow in Ruisseau Isabelle. Located between "Au P'tit Ruisseau" bar and Bas-Caraquet welcome sign_ just past Rue de Paul. A By-pass was reported at this LS on July 7_ 2012 and a temporary closed status zone was put in place.      ' to output field Comment
    Failed to execute (FeatureClassToFeatureClass).

The value that FeatureClassToFeatureClass could not write is from the 14th row in the CSV. I assume FeatureClassToFeatureClass proceeds sequentially through the CSV, which if correct tells me it had no problem with any of the cell values in the previous lines.

Here are the first 14 rows from my CSV's "Comment" field:

Comment
Foo Bar - overflow to vegetated surface      
overflow to vegetated landscape      
overflow to field      
overflow to field      

Flows into adjacent brook.  station is close to lagoon      this lift station serves a few homes - generally will pump 30 mins over two days.  have 2 pumps 

Lift station #6 (Rue de l'Ile). No standby power_ dual pumps_ total capacity 500 gal. per min. Overflow pipe to Caraquet Harbour. In case of failure_ alarm system sends signal to pager of LS operator.      
LS situated on Chemin St Simon_ No. #7. Located in a small baby barn building off road.      
Propane backup; next to lagoon;       
October 2015: Back up will be installed at this site this year.   Small grey shingled shack with RV dump station beside LS. Overflow_ near Ruisseau Isabelle. At small picnic park      
   October 2015 - lite and buzzer  small station - handles 3 houses and small campground.      Lift station on Rue François Gionet. UTM estimated. This LS pumps at LS Parc Fondateur.      less than 1 %
october 2015 - no lites just pager;  small flow station and no backup   White lift station with siding. Situated at 2264 rue Industrielle. No red warning light_ nor backup power observed.      
Blue & white trim lift station on "rue du Lac". New sewage system being installed at time of inspection. LS services 24 houses on Rue du Lac. LS does not have an apparent overflow pipe nor backup power.      
LS_ blue with white trim_ has overflow in Ruisseau Isabelle. Located between "Au P'tit Ruisseau" bar and Bas-Caraquet welcome sign_ just past Rue de Paul. A By-pass was reported at this LS on July 7_ 2012 and a temporary closed status zone was put in place. 

It has been suggested that the problem is with the single quotation mark in "P'tit", as shown in the error message.

The only problem is P'tit is part of a legitimate place name entered by someone else. I need every character in the comments field to be transferred.

Further, it has been suggested that the problem is not so much with FeatureClassToFeatureClass as it is the CSV converter not escaping single apostrophes.

If this is the problem what coding will force the CSV converter to escape single apostrophes?

Answer

It's nothing to do with the single quote in the 14th row comment field. Your 7th row also has a single quote and ArcGIS happily accepts that.

Your comment string that's failing is 259 characters long. It looks like ArcGIS 10.2 is creating the output FC Comments text field with a width of 256 (or 254 the shapefile limit).

When I run your code in ArcGIS 10.3, it works as expected as the output FC Comments text field is created with a width of 8000 characters by default.

Option 1. Upgrade ArcGIS.

Seriously! 10.2 is sooo old.

I understand if you can't though. We are stuck at 10.2 because of reliance on an ancient version of ArcSDE and Oracle. I had to get special permission to have 10.3 installed. Sigh.

Option 2. Use the field map parameter.

I managed to reproduce your error with this code:

arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, field_mapping="""Comment "Comment" true true false 256 Text 0 0 ,First,#,Temp_Layer,Comment,-1,-1""")

i.e.

c:\Python27\ArcGIS10.3\python test.py
Traceback (most recent call last):
  File "test.py", line 13, in 
    arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, field_mapping="""Comment "Comment" true true false 256 Text 0 0 ,First,#,Temp_Layer,Comment,-1,-1""")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\conversion.py", line 1789, in FeatureClassToFeatureClass

    raise e
arcgisscripting.ExecuteError: ERROR 001156: Failed on input OID -1, could not write value 'LS_ blue with white trim_ has overflow in Ruisseau Isabelle. Located between Au P'tit Ruisseau bar and Bas-Caraquet welcome sign_ just past Rue de Paul. A By-pass was reported at this LS on July 7_ 2012 and a temporary closed status zone was put in place.  ' to output field Comment
Failed to execute (FeatureClassToFeatureClass).

The following should work for you:

etc...
Description = arcpy.Describe(Target_FC)
Target_FC_name = Description.name
Target_GDB = os.path.split(Description.catalogPath)[0]
Temp_Layer = "Temp_Layer"

SR = Description.spatialReference
arcpy.MakeXYEventLayer_management(Source_CSV, "Lng", "Lat", Temp_Layer, SR, "")
# arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, "", "", "")
arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, field_mapping="""Comment "Comment" true true false 8000 Text 0 0 ,First,#,Temp_Layer,Comment,-1,-1""")

How does IP geocoding work?

I'm trying to understand how using each part of an IP address influences the accuracy of geocoding. For example, if

15.201.xxx.xxx

is used instead of

15.201.123.xxx

what is the approximate distance in KM? I'm using the IP2Location Lite database to examine who has visited a website. The README states the following:

1. IP2Location Lite Edition is free package with accuracy up to Class C (192.168.1.X) only. It is restricted for non-commercial use.

However, I'm having difficulty finding documentation on what the accuracy of this actually means. Does this accuracy vary over time also?

Answer

"Class C" is a historical reference meaning that your location data is summarized to the /24 level. For example, even if 203.0.113.1 and 203.0.113.254 happen to be in very different physical locations, the database will have only one lat/long for the entire 203.0.113.0/24 network.

Potentially useful generalizations about accuracy:

• 
  

  If you summarize to the /24 level (203.0.113.x), you can be reasonably confident that more than 80% of your data points will fall within the correct general metropolitan area.

  
  


• 
  

  If you summarize to the /16 level (203.0.x.x), you can expect that more than 80% of your location data points will be on the correct continent and probably even in the correct country.

  
  


• 
  

  If you summarize to the /8 level (203.x.x.x), your location data will be almost completely meaningless.

  
  

IP addresses that are likely to geocode reasonably well (within 25 miles or less) at the /24 summary level include:

• College/university campuses


• Corporate headquarters


• Cable modem connections serviced by major carriers like Comcast


• Fiber-to-the-premise connections serviced by major carriers like Verizon and Frontier

IP addresses that are reasonably likely to geocode quite poorly (off by 200 miles or more) at the /24 summary level include:

• 3G and 4G wireless phones/hotspots


• Corporate branch offices


• DSL connections serviced by CLECs

clustering - Grouping points in QGIS?

I need to group 400 points in 10 groups with 40 points each. The goal is to get the 40 points geographicly together like this:

Any idea how to do this with QGIS?

raster - A tool in QGIS like region group

I am trying to find a tool in QGIS that does the same function as Region group in ArcGIS.

What I am trying to do is to identify uniquely each patch of a raster that has a certain value (e.g "Value" =1) using neighbouring cells.

From the region group docs:

For each cell in the output, the identity of the connected region to which that cell belongs is recorded. A unique number is assigned to each region.

line - Remove polyline vertex outside polygon in Arcpy

I am trying to remove polyline vertex which are outside the buffer polygon and not creates new vertex on intersection.

I found a script which remove vertex but inside polygon. What is more it creates new vertex on intersection. It may help with this task.

import itertools, arcpy
arcpy.env.overwriteOutput = True
line_FC = arcpy.GetParameterAsText(0)
poly_FC = arcpy.GetParameterAsText(1)
output = arcpy.GetParameterAsText(2)

#clone input lines to memory for fast processing
temp = arcpy.CopyFeatures_management(line_FC, "in_memory/temp")


polygons = arcpy.CopyFeatures_management(poly_FC, arcpy.Geometry())

for poly in polygons:
    with arcpy.da.UpdateCursor(temp, ["SHAPE@"]) as uCursor:
        for line in uCursor:
            line = line[0]
            diff = line.difference(poly)
            #if the two lines are not equal, that means it intersected the polygon
            if not line.equals(diff):
                #the result of geometry.difference() is a multipart line of only those

                #parts that lie outside the polyon
                parts = diff.getPart()

                #if parts is empty that means the line is completely within the polygon
                #i.e., no difference
                if parts:
                    #We'll need to "join" the end of part1 to the beginning of part2
                    #so we'll just flatten the list of lists
                    joined = list(itertools.chain.from_iterable(parts))


                    #and create a new polyline object to update the shape
                    poly_trimmed = arcpy.Polyline(arcpy.Array(joined))
                    uCursor.updateRow([poly_trimmed])

arcpy.CopyFeatures_management(temp, output)

Is there a way to do it in arcpy? I am using ArcGIS Pro 1.3.

Fixing geometry validity errors in QGIS 2?

I am working on ownership data at a US county level. I have one shapefile per county that contains thousands of tax parcels. I need to dissolve the polygons representing parcels owned by the same individual. When the shapefile is not clean, the "Dissolve" tool in QGIS doesn't seem to work, it freezes. I have to fix the geometry validity issues first using the "Check Geometry Validity" tool. However, I have sometimes hundreds of errors. It takes a lot of time to fix them manually. Is there any other strategies that could save me some time?

I am using QGIS 2.4.0.

shapefile - Why OSM data has different attribute values from the same element and same time frame?

Here is my question: Does anyone know why OSM has different attribute values from the same element and same frame time?

The example explains itself, when I try the Shapefile from different mirrors like Planet.osm.org , bbbike.org and geofabrik.de; I find that they have a mistake with the Direction of Travel on its attribute value.

When I check in the website: here the link , I see on its attribute value that there is no mistake at all and has its right value.

So, I thought that those mirrors extract the exact information from osm servers, but I think that there is a lack of info. Is that right?

Thanks for your time,

Answer

Metadata is data about data - the coordinate system of a layer, what the allowable attribute values are, or what a given attribute is/means. In a shapefile (which is a database) what you're asking about is an attribute. However in the OSM data this information is stored as tags because it is xml (Extensible Markup Language) document and not a database. See this related question: shapefile terminology for key-value pair

As for your question, it appears that something in the extraction process that converts the xml to a database is dropping the negative sign. This could be by design or a bug / error. You would have to contact one of the organizations doing the conversion and go through their help/support system to report or determine what is going on (unless somewhere here knows something - you may want to edit your question to use the correct terminology).

The alternative would be to process the OSM data to a shapefile yourself, so you know exactly what conversions are being made.

arcpy - How to update attributes of a shapefile using list values?

I am using an update cursor to update the attributes of a shapefile using values from a list. It seems very easy, but I couldn't make it run correctly. When I run the following code, it updates all rows based on the final value in the list. Please suggest how to fix it.

cur1 = arcpy.UpdateCursor(fc)
for v in cur1:
    for s in list:
        val = s

        v.Val = val
        cur1.updateRow(v)

Answer

This looks like a simple logic error. I'm assuming that you want the first feature in the cursor to be updated with the first item in the list, then the second feature with the second item, etc. As written, you're updating the same row over and over again with each element of the list. (Then repeating the process for each row.)

Maybe your loop should look something like this instead:

listIndex = 0
cur1 = arcpy.UpdateCursor(fc)
for c in cur1:
    if listIndex >= len(list):
        arcpy.AddMessage("More rows than items in the list.  Quitting early.")

        break
    val = list[listIndex]
    v.Val = val
    cur1.updateRow(v)
    listIndex = listIndex + 1

EDIT: Corrected logic error (needed to check array length before accessing it, not after). Also fixed missing : on if statement.

EDIT: Following the suggestion in Paul's comment is more "Pythonic" and reduces the lines of code:

cur1 = arcpy.UpdateCursor(fc)
for tuple in enumerate(cur1):

    index = tuple[0]
    row = tuple[1]
    if index >= len(list):
        arcpy.AddMessage("More rows than items in the list.  Quitting early.")
        break
    row.Val = list[index]
    cur1.updateRow(row)

fields attributes - Writing vector layer to CSV file with geometry using PyQGIS?

I use Python try to write vector layer file to csv file with this code:

QgsVectorFileWriter.writeAsVectorFormat(mylayer, r'c:\temp\xyz.csv', "utf-8", None, "CSV")

It can export to xyz.csv but only attributes show in the csv, not the geometry column.

How can I export both attribute and spatial data into csv file?

How do I remove a join (between a standalone table and a feature layer) with ArcObjects?

I already know how to do a table join with ArcObjects, via the IDisplayTable, IMemoryRelationshipClassFactory etc. interfaces. Here are some links to resources that have example code for this:

What I need to know is, how do I undo, i.e. remove such a join?

---

Once I've joined a standalone table to a feature layer, I end up with references to:

• 
  

  an IRelationshipClass (resulting from the join operation);

  
  


• 
  
  

  an ITable/IStandaloneTable/IDisplayTable (the table that was joined to the feature layer); and

  
  


• 
  

  an IFeatureLayer (the feature layer to which the table was joined).

  
  

Do these interfaces, or rather the objects behind them, allow unjoining at all? I've seen that IRelationshipClass has various DeleteRelationship… methods, but I can't see how they could work towards that end.

Answer

Code below works for me for a featurelayer. Similar logic could be used for a standalone table.

private void RemoveAllJoins(IFeatureLayer fLayer)

{
    var dispTable = fLayer as IDisplayTable;
    var rqt = dispTable.DisplayTable as IRelQueryTable;
    if (rqt != null)
    {
        Debug.Print("source: {0}", ((IDataset)rqt.DestinationTable).Name);
        Debug.Print("dest: {0}", ((IDataset)rqt.SourceTable).Name);
        fLayer.FeatureClass = (IFeatureClass)rqt.DestinationTable;
    }
    else

        Debug.Print("there are no joins");
}

The documentation for IRelQueryTable at 10 seems to be missing this important graphic, that appears in 8.3 help doc:

qgis - Multiple layers auto-updated when one of them is modified

I just discovered this post and I found what I searched here. But I was wondering if there was a way to have multiple layers (let's say l1, l2, l3, l4) being updated automatically when you change something on one of them.

For example if I change the outline style of l1 from solid to dot I would like to see l2, l3 and l4 updated the same way without copy/pasting the style like mentioned in the link

(I'm using qgis 2.18.28 on windows)

Making circular polygon in QGIS?

I am trying to make a floor plan for my office, have managed to make all of the rectangular components (e.g. the rooms) however some of the desk layouts are circular and I haven't been able to work out how to make exact and measured circles. They need to have a radius of 1.45m.

Answer

Take a look at the plug in CADDigitize. It may be helpful. edit : This tool allows the user to draw circles, squares, resctangels, ellipses from many ways.

Selecting column like Primary Key in View from PostGIS using PyQGIS?

Postgresql doesn't support Primary Key in views. When I load a postGIS view in QGIS "manually", I can choose the column that will be the Primary Key, but I would like to know if I can select a column as Primary Key when I load it using pyQGIS script.

Answer

I have found the answer reading the QGIS API

This is the method:

void QgsDataSourceURI::setKeyColumn (QString column )

this is my example:

uri = QgsDataSourceURI()

uri.setConnection('localhost', '5432', '-mydatabase-', '-myowner-', '-mypasswd-')
uri.setDataSource('public', 'parcelas', 'the_geom')
#public: schema, parcelas: the view, the_geom: spatial column
uri.setKeyColumn('id') #name of my primary key

vlayer = QgsVectorLayer(uri.uri(), '-name_of_the_layer-', 'postgres')
if vlayer:
  print "All Ok!"

QgsMapLayerRegistry.instance().addMapLayer(vlayer)

Disabling caching in geowebcache.xml using GeoServer?

I need to publish a geometry table (in PostGIS) as a layer of vector tiles (in PBF format) in GeoServer. As the table is dynamically written by a production web app, i need to set the layer's caching as disable. I got a hint that this could be done by editing the geowebcache.xml, specifically with and tags, using the special value -1 (that means no-cache).

What i can not find anywhere is a a real example of configuring this in the geowebcache.xml. The following is my default file (of GeoServer 2.13.2 on Apache Tomcat 8.5):

  xmlns="http://geowebcache.org/schema/1.8.0"
  xsi:schemaLocation="http://geowebcache.org/schema/1.8.0 http://geowebcache.org/schema/1.8.0/geowebcache.xsd">
  1.8.0
  120
  
  
    

    
  

  
    

  

Do i:

1. 
   

   simply put the and out of comment?

   
   


2. 
   

   Or do i need to make a whole block for that specific layer (from , , , , all the way down to ), that specifies and ?

   
   
   

qgis - When I finish drawing a polygon, it disappears again (although appearing in my attribute table)

I'm working with Quantum GIS at this moment. And as the title says, when I finish drawing a polygon, it disappears again. The new object does appear in my attribute table and in the legend/'filling of the polygon' should not be the problem. Before this never happend, so I supposed it was some change I made in the options. But even after reinstalling the program the problem stays. After searching for a while I discovered that on a new project it works correctly till the moment I import a vector-file I made before with line objects. How can importing a file have such an influence?

Thanks in advance!

Need Basic Viewer template for ArcGIS JSAPI 3.2/3.3

I am looking for basic viewer for ArcGIS JSAPI but not getting proper sample/example.

I am going through ArcGIS resource also but looks like its not working. Going through this question also.

Checking out layout example and getting confused that what is difference between Viewer template and layout?

So any help regarding this will be great!

Answer

There is no viewer template for the Javascript API, at least not one offered by Esri. But... Esri developer and software engineer David Spriggs has been working on an open source project to produce one: https://github.com/DavidSpriggs/ConfigurableViewerJSAPI

qgis - Why can't I save my custom CRS?

I am unable to create custom CRSes in QGIS. The 'save' button doesn't seem to have any effect. My best guess is that I don't have permission to save in the folder or database where custom CRSes go. I have looked at the screenshots of other users who are having trouble with this, and I usually see "1 of 1" in the custom CRS browser. On my system, I see "* of 0" instead

Any assistance appreciated.

QGIS version: 1.8.0

Ubuntu 12.04 LTS

Answer

I was able to solve this problem on my own. After some investigation, I found that the ~/.qgis folder was owned by root and had permissions 644. With these permissions, the qgis.db cannot be created when QGIS runs, and settings cannot be saved. There are two ways to fix this problem.

1. 
   

   Delete the ~/.qgis folder. It will automatically be recreated when you run QGIS, but with you as owner and permissions 776.

   
   


2. 
   

   Manually modify the ~/.qgis folder using the chown, chgrp, and chmod commands so that you are the owner and permissions are 776.

   
   

Based on issue ticket #914, this problem has been fixed in the source: http://trac.osgeo.org/osgeo/ticket/914

Installing postgis on postgresql 9.1 not working

I'm trying to install Postgis 2.0 on a perfectly functionig postgresql 9.1.10 database.

I installed the package via:

richard@srv1:/$ sudo apt-get install postgresql-9.1-postgis postgis
Reading package lists... Done
Building dependency tree       
Reading state information... Done
postgis is already the newest version.
postgresql-9.1-postgis is already the newest version.
0 upgraded, 0 newly installed, 0 to remove and 2 not upgraded.

With out any errors. Then I try to add postgis as an extension to an existing (empty) database, with:

 dev=# CREATE EXTENSION postgis;
 ERROR:  could not open extension control file "/usr/share/postgresql/9.1/extension/postgis.control": No such file or directory

If I look at the location the file does not exist.

What am I missing? I can't figure it out. If I google I only get ancient results....

Where can I find the version of the installed postgis? As I'm not able to do a SELECT PostGIS_version();

Environment: Postgresql 9.1.10, Ubuntu 12.04

Update

I'm willing to install postgis2.0; apparently I have 1.5 installed. After purging that version I added the repository

deb http://apt.postgresql.org/pub/repos/apt/ precise-pgdg main

After trying to install the package via

sudo apt-get update
sudo apt-get install postgis

it ends with:

Errors were encountered while processing:
/var/cache/apt/archives/libgdal1_1.9.0-3.1~pgdg12.4+1_amd64.deb
E: Sub-process /usr/bin/dpkg returned an error code (1)

Update 2

It appears because of the naming of the packages, it installed the wrong architecture type package. But I ended up using this website and these steps after making sure I had gdal>= 1.9.0 and GEOS >=3.3.3:

wget http://postgis.refractions.net/download/postgis-2.0.0.tar.gz
tar xfvz postgis-2.0.0.tar.gz
cd postgis-2.0.0
./configure --with-gui

make
sudo make install

sudo ldconfig
sudo make comments-install


sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/shp2pgsql
sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/pgsql2shp
sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/raster2pgsql

TADA

durham=#  select PostGIS_full_version();

NOTICE:  Function postgis_topology_scripts_installed() not found. Is topology support enabled and topology.sql installed?
                                                                     postgis_full_version                                                                     
--------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER
(1 row)

Answer

By default, you don't have PostGIS 2.0 but 1.5 (see http://packages.ubuntu.com/fr/precise/postgresql-9.1-postgis)

By doing CREATE EXTENSION postgis; you're trying to initialize from PostgreSQL with the new loading mechanism a PostGIS version (1.5) that do not support it and in this case, you will need to follow the PostGIS 1.5 official doc

If you need 2.0, follow various instructions already available in previous gis.stack.exchange.com answers e.g Install PostGIS on Ubuntu 12.04 or Can POSTGIS 2.0 be install with POSTGRES 9.2 on Ubuntu 12.04 OS at the moment?

Edit due to changes in question:

When I look in the repository (http://apt.postgresql.org/pub/repos/apt/ precise-pgdg main), I see you can get different PostGIS versions for different PostgreSQL versions, so sudo apt-get install postgis is not enough clear for the package manager

Did you try?

sudo apt-get install postgresql-9.1-postgis-2.0

modelbuilder - arcpy Tool works in debug mode but not in normal model

When I try to run the following simple script from arcmap it fails but when run using right--click debug and then F9 (debug) in python scripter it works fine. I don't understand why???

# Modified by: Georgec Corea (ATGIS) coreagc@gmail.com

# Author:  ESRI
# Date:    July 5, 2010
# Version: ArcGIS 10.0
# Purpose: This script will print one or more map documents to a local printer.
#          The script is intended to run within a script tool.  There are two

#          parameters:
#               1) Select Map Documents to Print,
#               2) Select Output Printer (auto populated using a validation script)
#
#Notes: The print order of the MXDs is based on how they are entered.  The MXD
#       at the top of the list is first followed by those below it.

import arcpy, string, datetime, shutil, os
import arcpy.mapping as MAP


#Read input parameters from script tool
MXDList = string.split(arcpy.GetParameterAsText(0), ";")
#printer = arcpy.GetParameterAsText(1)

count=0
#Loop through each MXD and print
for MXDPath in MXDList:
    count=count+1
    arcpy.AddMessage(str(count)+'. Working on '+str(MXDPath)+' @ '+str(datetime.datetime.now())+"\n")
##    try:

    mxd=arcpy.mapping.MapDocument(MXDPath)
    outPDFpath = mxd.filePath[:mxd.filePath.rfind('.')]+'.pdf'
    ddp = mxd.dataDrivenPages
    #MXD = MAP.MapDocument(mxd)
        #arcpy.AddMessage(ddp, MXD)
    if str(os.path.isfile(outPDFpath))=='False':
            outputPDF = arcpy.mapping.PDFDocumentCreate(outPDFpath)
            mxd.dataDrivenPages.exportToPDF(outPDFpath)
        #arcpy.mapping.ExportToPDF(mxd, outPDFpath)
            name='c:\\temp\\'+outPDFpath[outPDFpath.rfind('\\'):]

            shutil.copy2(outPDFpath,name)
            arcpy.AddMessage('Printed '+str(outPDFpath)+' @ '+str(datetime.datetime.now())+"\n")
    else:
            arcpy.AddMessage('File Exists '+str(outPDFpath)+' skipped pdf creation @ '+str(datetime.datetime.now())+"\n")
##    except:
##        errorm=arcpy.GetMessages()
##        arcpy.AddMessage(' Error creating file...' +str(outPDFpath)+'\n'+str(errorm)+'\n'+' continuing...')

    #MAP.PrintMap(MXD, printer)
    #Remove variable reference to file

del mxd

When run in debug

Answer

I changed the tool to run only in the foreground and after that it works fine. It would be nice to understand what part of the script can not be run in background mode in a tool.

Thanks Get Spatial and Nicksan for your input.

routing - Add time restriction in ArcGIS Network Analyst

I'm trying to use train tracks as a network in ArcGIS Network Analyst.

We are trying to create a model for security agents on board our trains and we'd like to optimize their daily "route". For now, it's all made by hand.

My approach is to take our gtfs data and create a network using stop times to determine when the train segment is "traversable". When there's a train on the segment, it's open, if not, it's restricted.

The problem is I can't set a restriction based on time in Network Analyst (or wasn't able to). I would like to be able to say by a formula or script than you can take that segment only from 7:03 AM to 7:10 AM for example.

Is it possible?

Thanks!

Answer

What you want to do is possible, but not in the manner you're thinking.

A network could generally be considered static in that the rules and restrictions for the edges are governed by single, static values. Accounting for change and time of day is another level known as Traffic.

You may want to read through this ArcGIS help page for a detailed explanation, but basically either storing a variety of values for a single edge or multiple networks with varying values would take up too much space. Instead, the travel time properties (or a profile) for an edge are referenced from a lookup table based on time of day.

What you need to do is create historical traffic data for your network based on the train schedule. There are several tutorials and help pages that can assist you with this, but covering the whole process would be too long for a SE question/answer.

qgis - How to remove duplicate point features?

I have created many .gpx datas with my Garmin device, but due to the amount of data, I have imported the mapping data from my device several times (also because I didn't want it to be lost)...

As a result I have now many double or even 6 times the same point (including tracks, which I displayed as points also, so that I have the elevation) after I saved all the gpx-files as shapefiles and then merged them.

Is there a way to automatically erase the duplicates with QGIS, or if not with any other free software working on Ubuntu 12.4?

I very much appreciate your help!

Answer

Try the 'mmqgis'-plugin. There, go to Modify -> delete duplicate geometries and select your layer.

QGIS 1.8 will not start on Windows 7

I have downloaded QGIS 1.8 and installed it on my windows 7 system. When I click on the desktop icon, nothing happens.

I am currently running QGIS 1.7.4 with no problems and the installs of previous versions have always worked with no problems. I'm not sure what's going on here. Has anyone else experienced this and have a solution?

Change vector layer symbology PyQGIS 3?

I'd like to change the symbology of my vector layer and can't figure it out how to read the documentation. I would like to have an polygone with an orange outline, a dashed fill style and an orange stroke color. However that would be nice to have an answer that explains how to access most of the symbology change! (if possible ;-) )

Here is what I managed so far:

import os
import processing
from qgis.core import *
from qgis.gui import *
import qgis.utils


# initialize project
project = QgsProject.instance()
project.setCrs(QgsCoordinateReferenceSystem(21781))
project.setFileName("C:\\test_symbology.qgs")

my_shp = "Y:\\shp_collection\\A300.shp"
my_shp_name = "A300"

# load polygon
shplayer = iface.addVectorLayer(my_shp , my_shp_name , "ogr")

if not shplayer.isValid():
    print(".shp failed to load!")
shplayer.setCrs( QgsCoordinateReferenceSystem(21781,QgsCoordinateReferenceSystem.EpsgCrsId))

layer = iface.activeLayer()

# make my polygone orange
single_symbol_renderer = layer.renderer()
symbol = single_symbol_renderer.symbol()
symbol.setColor(QColor.fromRgb(255,128,0))

How can I set the other symbology properties?

Answer

The easiest for me would be to check the properties of the symbology layer which contains symbol properties affecting the whole layer:

dir(symbol)

>>>['DynamicRotation', 'Fill', 'Hybrid', 'Line', 'Marker', 'RenderHint', 'RenderHints', 'ScaleArea', 'ScaleDiameter', 'ScaleMethod', 'SymbolType', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_getLineString', '_getPoint', '_getPolygon', '_getPolygonRing', 'appendSymbolLayer', 'asImage', 'bigSymbolPreviewImage', 'changeSymbolLayer', 'clipFeaturesToExtent', 'clone', 'cloneLayers', 'color', 'createSimple', 'defaultSymbol', 'deleteSymbolLayer', 'drawPreviewIcon', 'dump', 'exportImage', 'forceRHR', 'hasDataDefinedProperties', 'insertSymbolLayer', 'layer', 'mapUnitScale', 'opacity', 'outputUnit', 'renderFeature', 'renderHints', 'renderPolygon', 'renderUsingLayer', 'renderVertexMarker', 'setAngle', 'setClipFeaturesToExtent', 'setColor', 'setForceRHR', 'setLayer', 'setMapUnitScale', 'setOpacity', 'setOutputUnit', 'setRenderHints', 'startRender', 'stopRender', 'symbolLayer', 'symbolLayerCount', 'symbolLayers', 'symbolRenderContext', 'takeSymbolLayer', 'toSld', 'type', 'usedAttributes']

To access individual symbol layers inside, you could use:

dir(symbol.symbolLayer(0))

>>>['Property', 'PropertyAngle', 'PropertyArrowHeadLength', 'PropertyArrowHeadThickness', 'PropertyArrowHeadType', 'PropertyArrowStartWidth', 'PropertyArrowType', 'PropertyArrowWidth', 'PropertyAverageAngleLength', 'PropertyBlurRadius', 'PropertyCapStyle', 'PropertyCharacter', 'PropertyCoordinateMode', 'PropertyCustomDash', 'PropertyDisplacementX', 'PropertyDisplacementY', 'PropertyDistanceX', 'PropertyDistanceY', 'PropertyFile', 'PropertyFillColor', 'PropertyFillStyle', 'PropertyGradientReference1IsCentroid', 'PropertyGradientReference1X', 'PropertyGradientReference1Y', 'PropertyGradientReference2IsCentroid', 'PropertyGradientReference2X', 'PropertyGradientReference2Y', 'PropertyGradientSpread', 'PropertyGradientType', 'PropertyHeight', 'PropertyHorizontalAnchor', 'PropertyInterval', 'PropertyJoinStyle', 'PropertyLayerEnabled', 'PropertyLineAngle', 'PropertyLineDistance', 'PropertyName', 'PropertyOffset', 'PropertyOffsetAlongLine', 'PropertyOffsetX', 'PropertyOffsetY', 'PropertyOpacity', 'PropertyPlacement', 'PropertyPreserveAspectRatio', 'PropertySecondaryColor', 'PropertyShapeburstIgnoreRings', 'PropertyShapeburstMaxDistance', 'PropertyShapeburstUseWholeShape', 'PropertySize', 'PropertyStrokeColor', 'PropertyStrokeStyle', 'PropertyStrokeWidth', 'PropertyVerticalAnchor', 'PropertyWidth', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_renderPolygon', 'angle', 'brushStyle', 'clone', 'color', 'copyDataDefinedProperties', 'copyPaintEffect', 'create', 'createFromSld', 'dataDefinedProperties', 'drawPreviewIcon', 'dxfAngle', 'dxfBrushColor', 'dxfBrushStyle', 'dxfColor', 'dxfCustomDashPattern', 'dxfOffset', 'dxfPenStyle', 'dxfWidth', 'enabled', 'estimateMaxBleed', 'fillColor', 'hasDataDefinedProperties', 'isCompatibleWithSymbol', 'isLocked', 'layerType', 'mapUnitScale', 'offset', 'offsetMapUnitScale', 'offsetUnit', 'ogrFeatureStyle', 'outputUnit', 'paintEffect', 'penJoinStyle', 'prepareExpressions', 'properties', 'propertyDefinitions', 'renderPolygon', 'renderingPass', 'restoreOldDataDefinedProperties', 'setAngle', 'setBrushStyle', 'setColor', 'setDataDefinedProperties', 'setDataDefinedProperty', 'setEnabled', 'setFillColor', 'setLocked', 'setMapUnitScale', 'setOffset', 'setOffsetMapUnitScale', 'setOffsetUnit', 'setOutputUnit', 'setPaintEffect', 'setPenJoinStyle', 'setRenderingPass', 'setStrokeColor', 'setStrokeStyle', 'setStrokeWidth', 'setStrokeWidthMapUnitScale', 'setStrokeWidthUnit', 'setSubSymbol', 'startRender', 'stopRender', 'strokeColor', 'strokeStyle', 'strokeWidth', 'strokeWidthMapUnitScale', 'strokeWidthUnit', 'subSymbol', 'toSld', 'type', 'usedAttributes', 'writeDxf']

Near the end of this list of properties, we can see there's settings for setBrushStyle and setStrokeColor. The setBrushStyle requires Qt.BrushStyle to set a brush style, we can check the properties inside Qt using:

dir(Qt)

So doing a little search, we can find brush patterns such as BDiagPattern or FDiagPattern

---

So now we can apply the relevant settings:

layer = iface.activeLayer()
single_symbol_renderer = layer.renderer()
symbol = single_symbol_renderer.symbol()
#Set fill colour
symbol.setColor(QColor.fromRgb(255,128,0))
#Set fill style
symbol.symbolLayer(0).setBrushStyle(Qt.BrushStyle(Qt.FDiagPattern))
#Set stroke colour
symbol.symbolLayer(0).setStrokeColor(QColor(255,128,0))

#Refresh
layer.triggerRepaint()
iface.layerTreeView().refreshLayerSymbology(layer.id())

tilecache - WMS tile cache that works with time dimension

Is there a good WMS tile cache that understands the WMS time dimension? And what configuration is necessary? Open source preferred.

I am serving some WMS layers with the time dimension, so the map requests from the client include dates. For example http://…?…request=GETMAP&…&TIME=2010-12-10 returns the map for 10 Dec 2010. In my undestanding the tile cache will need to treat the TIME parameter as an identifier for the maps, just like the layer name.

I'm aware of Geowebcache, Tilecache, and MapProxy but I can't find any reference to the time dimension in their documentation. The time dimension has existed since WMS 1.1.0 so I would hope it's been considered by someone? Our layers are served by ncWMS and MapServer, if that's relevant. We would like to use a tile cache to increase the speed for our most frequently requested content. Can anyone recommend a tile cache, and preferably point us to the relevant documentation too?

arcgis desktop - Select "interior" lines of a polygon-to-line layer

I have a polygon featureclass that has adjacent features. Using Arc10, I am turning those polygons into lines using the "Featureclass to Lines" tool.

Now I would like to select ONLY THE INTERIOR line(s) (color=grey) so that I can run some ADDITIONAL processes on them. The key is that some (but not all) feature may touch the boundary line (color=blue). I am trying to use "Select By Location", but I cannot come up with a way to accomplish this task.

I might just be overthinking this one, but I am stuck. Any suggestions?

Answer

In doing some brainstorming with a colleague, we finally came up with a solution to this problem.

1. We exported all polygons to lines


2. Dissolved the original polygons and preserved no attributes and not allowing multi-part features (giving the fewest # of single-part polygons possible).


3. Select by Location, using the lines as the input layer, the dissolved polygons as the selection layer, and "SHARE_A_LINE_SEGMENT_WITH" as the selection method. This selected all features that were on the boundary.



4. Switch the Selection. This gave us all interior lines that did not match the outer dissolved boundary.

While I have not tested it for 100% accuracy, it seems to have accomplished the task that we were after very well. Of course, to automate this task, I just dropped all of these functions into a ModelBuilder routine, which makes the task much less tedious for future operations.

arcgis desktop - Storing points from multiple GDBs as object

I often download data that comes in 30 minutes to 1-hour windows and then have to compile that data manually into one feature class, from multiple GDBs. I am trying to iterate through the separate GDBs and store each point feature class in an object. Here's what I have so far (Python 2.7.14):

import arcpy
import os
from arcpy import env
outputOverwrite = True

print('imported')


arcpy.env.workspace = r'C:\Vector\20190902'
print('ws set')

workspaces = arcpy.ListWorkspaces(workspace_type = 'FileGDB')

for workspace in workspaces:
   datasets = arcpy.ListDatasets(feature_type = 'All')

for ds in datasets:
   arcpy.env.workspace = ds

   fc = arcpy.ListFeatureClasses(feature_type = 'Point')
   print(fc)

I am unable to post pictures or copy and paste my code, as the system is on a totally different network.

Edited to reflect changes made by @BERA.

This code outputs

imported
ws set
>>>

It doesn't print the feature classes like i expected it to.

I traced this back to the line of code

for workspace in workspaces:
    datasets = arcpy.ListDatasets(feature_type='All')

When I tried to do

print(datasets)

It gave me blank lists

[]

[]
[]
[]
...

Hope this clears up the ambiguity on my end.

Answer

ListFeatureClasses will search in current workspace, which is r'C:\Vector\20190902'. You need to change it to point at each gdb when listing feature classes.

I dont know why you are listing datasets, are you sure that is what you want?

If you have a structure like this, code below will work:

import arcpy, os

arcpy.env.workspace = r'C:\GIS\ArcMap_default_folder'

fclist = []
for ws in arcpy.ListWorkspaces(workspace_type='FileGDB'):
    arcpy.env.workspace = ws
    for fc in arcpy.ListFeatureClasses(feature_type='POINT'):
        fclist.append(os.path.join(ws, fc))


print (fclist)
#arcpy.Merge_management(inputs=fclist,...

['C:\\GIS\\ArcMap_default_folder\\Default.gdb\\point1', 'C:\\GIS\\ArcMap_default_folder\\test.gdb\\point2']

pyqgis - How can I trigger a QGIS plugin Event from map refresh?

I have a combobox in my plugin (python Qt dockable widget) that needs to contain the currently open rasters.

How do I trigger an update of the plugin when something new gets added to the map?

Answer

There's a couple of methods which I use (note that it may not be the most efficient solution):

• 
  

  (using Combo Box) - Add the following line before you populate your combo box:

  
  

  self.dockwidget.comboBox.clear()
  

  
  

  Then when you add a new layer, close your plugin and open it again, the combo box should be re-populated with the new additions.

  
  

---

• 
  

  (using QgsMapLayerComboBox) - This is automatically updated when any changes are made to the layer list (i.e. layers added/removed, names changed etc) without you having to refresh the plugin. I would recommend this method. You can set it so that it only recognises raster layers by using something like:

  
  

  self.dockwidget.mapLayerComboBox.setFilters( QgsMapLayerProxyModel.RasterLayer )
  

---

Note: When using QGIS custom widgets, you may need to edit your dockwidget_base.ui file, find the QgsMapLayerComboBox class (usually at the end of the file) and replace the

to qgis.gui.h.