In ArcGIS, there is a tool called "Identity." Say you have a bunch of points. Some points fall in polygons of interest, some don't. Each polygon has a unique ID. This tool will add a field to the points and assign the unique ID of the polygon that contains it. Is there a tool that does this in QGIS?
This question does not ask the same question I am.
I'm looking for a free high resolution map that clearly indicates the land and sea in the world. I need it to calculate distances over sea between ports. Additionally it should also include inland waterways and canals.
Does anybody know whether this data is freely available?
Answer
After some trial and error the most practical solution turned out to be:
Use a shape file of all world countries, add the Suez and Panama canal manually using software like QGIS, add the ice of the arctic and Antarctica using shapefiles as well. This gives you a map of all the places a ship can go.
I left out the inland waterways in the end. However these could be included using the link provided by Aaron.
I've got this csv file that I'm trying to import to QGIS 3.2. It's just a data table with no geometry attached.
When I just drag it from the Browser on the left-hand side, all the columns get incorrectly read in as string. However, when I use the "Open Data Source Manager" and import the file through there, the column types are all correct (first one is integer and the others are double). But I really didn't specify that the columns should be double instead of string.
So, finally, my question is: how/where can I try to force certain column type at import for each of my columns?
Also, extra brownie points to whoever knows what's going on in the background that makes dragging a CSV from the Browser different from using the "Open Data Source Manager".
Answer
QGIS uses OGR in the background, and OGR interprets all columns as strings.
The OGR CSV driver returns all attribute columns as string data types if no field type information file (with .csvt extension) is available.
Using CSVT file with the same name as the CSV file allows you to specify types of columns in a CSV file.
Limited type recognition can be done for Integer, Real, String, Date (YYYY-MM-DD), Time (HH:MM:SS+nn), DateTime (YYYY-MM-DD HH:MM:SS+nn) columns through a descriptive file with the same name as the CSV file, but a .csvt extension. In a single line the types for each column have to be listed with double quotes and be comma separated (e.g., "Integer","String"). It is also possible to specify explicitly the width and precision of each column, e.g. "Integer(5)","Real(10.7)","String(15)".
Difference between dragging a CSV from the Browser and using the "Open Data Source Manager" is that "Open Data Source Manager (ODSM)" guesses data types of columns.
Sample .csv and .csvt file:
sample.csv
1,high,2.3
2,low,5
3,low,7.2
sample.csvt
"Integer(5)","String(10)","Real(3.5)"
If you import the file to QGIS using Browser or ODSM, the attribute table looks like this:
Further information, please look at Comma Separated Value (.csv)
I am new to R and using the raster package. I have a problem extracting polygons from an existing raster file. If I use
extract(raster, poly_shape)
function on the raster it always creates a list with the data. What I really want is to extract another raster file that I can load with ArcGIS again. After reading a bit more I think the crop function is what I really need. But when I try to use this function
crop(raster, poly_shape)
I get this Error:
Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap
The files raster and poly_shape are the same for both functions. Can you tell me what could be wrong here? Is it even right that the crop function creates another raster and not a list?
EDIT: The extent() function does not work for me. I still get the same error. But I am sure the 2 datasets overlap! With the
extract(raster, poly_shape)
I get the right data from it. Just as a list and not as a raster like I want to have it. I just loaded the datasets in ArcGIS before and they fit very well so I did not check the projection. Now I tried
projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"
and you can see that the projections do not fit. The extract function seems to be able to automatically transform the files in the right way. I know that because I did the following:
The 2 have the exact same result so I guess the conclusion should be that the extract function did work correct. Now I have 2 options I guess:
What would you suggest to do here?
Answer
What I actually searched for was the mask() function.
mask(raster, poly_shape)
works without errors and returns what I searched for.
I have a string field 'Oneway' with three values '+' '-' and Null showing a street accessibility, I would like to replace them with values 'F' 'T' and null as the script for Network dataset expect these values.However it seems that it is not possible to replace values if Null values are present. I tried also changing the field type to double or text but the field properties are not editable somehow (I tried this in arc catalog). I am only able to replace value '+' and '-' with numbers using the field calculator but need to replace them with letters.
I plan to run a scheduled server job on a feature class in an Oracle 18c/10.7.1 geodatabase.
The job would:
Note: I would like to do this without needing to temporarily store any sort of intermediate results in a temporary feature class/GDB.
I've gotten part-way there:
import arcpy
conn = "Database Connections\\DEV.sde\\"
input_lines = "OWNER.ROAD"
existing_fc = "OWNER.SCHED_ROAD_MIDPOINT"
arcpy.TruncateTable_management(conn + existing_fc)
arcpy.FeatureVerticesToPoints_management(conn + input_lines, conn + existing_fc, "MID")
#arcpy.Append_management(, , "TEST", "", "")
Question:
I'm not sure how to use FeatureVerticesToPoints_management as the input in the append process.
Is there a way to append/insert geoprocessing results into an existing feature class using ArcPy?
Answer
I've built on @user2856's answer:
Caution! The target FCs will get truncated! All records in the target FCs will be deleted.
import arcpy
conn = "Database Connections\\\\"
#This script only works for polyline feature classes (not polygons or points).
#Caution! The target FCs will get truncated! All records in the target FCs will be deleted.
# INPUT FC TARGET FC ID FIELD (for both the input FC and the target FC)
lstFC = [\
[".", ".S_SIDEWLK_MIDPOINT", "SDW_ID" ],\
["." , ".S_SEWER_MIDPOINT" , "SEWER_ID" ],\
["." , ".S_ROAD_MIDPOINT" , "ROAD_ID" ],\
]
def ReloadMidpoints():
for fc in lstFC:
inputFC = conn + (fc)[0]
targetFC = conn + (fc)[1]
idField = (fc)[2]
arcpy.TruncateTable_management(targetFC)
print "<" + targetFC + "> has been truncated."
with arcpy.da.SearchCursor(inputFC, ['SHAPE@', idField]) as s_rows:
with arcpy.da.InsertCursor(targetFC, ['SHAPE@', idField, "MIDPOINT_X","MIDPOINT_Y"]) as i_rows:
for row in s_rows:
Midpoint = row[0].positionAlongLine(0.50,True).firstPoint
rowVals = [(Midpoint.X,Midpoint.Y),s_rows[1],Midpoint.X,Midpoint.Y]
i_rows.insertRow(rowVals)
print "<" + targetFC + "> has been reloaded with midpoints that were generated from <" + (inputFC) + ">."
print ""
ReloadMidpoints()
print "Complete."
I am struggling with the DescribeFeatureType query with WFS KVP encoding.
I have this WFS: https://inspire-wfs.maanmittauslaitos.fi/inspire-wfs/au?request=GetCapabilities&service=wfs&version=2.0.0 and I can see in the GetCapabilities that there is a FeatureType called au:AdministrativeUnit.
I could of course use DescribeFeatureType to give me the description of all Feature Types of the service, but I would like to use &typeName=XXXX in the kvp to choose only this FeatureType au:AdministrativeUnit
So I basically do the following request: https://inspire-wfs.maanmittauslaitos.fi/inspire-wfs/au?&request=DescribeFeatureType&service=WFS&version=2.0.0&typeName=au:AdministrativeUnit
But the result is not filtering at all to give me only the Feature Type with the name au:AdministrativeUnit.
I saw that in the OGC specifications, there is sometimes written typename or typenames. I tried both of them with changing the version and I don´t get it!
Am I missing something important here?
With an other WFS, I can do the query with filtering typenames without issue!
I am relatively new to Linux, so can you indicate the exact commands (step-by-step) that I need to type to get the latest QGIS installed on Ubuntu 14.04.
I have tried to follow the instructions provided here https://www.qgis.org/en/site/forusers/alldownloads.html but I haven't managed to install anything.
Can someone explain to me step by step process of doing it. Right from adding the QGIS repository on the sources list.
I have Ubuntu 14.04.1 LTS installed on a 32-bit Dell Latitude E4310 Machine
I am trying to make a report with the QGIS composer. This report contains multiple maps. I don't have any trouble to make it, I just have to lock the maps with a right click to avoid them being refreshed while I change the layers in QGis.
However, when I decide to export the report as a pdf, the result is a pdf with copies of the same map. It looks like the "export as pdf" tool refreshed all maps.
How can I avoid that and have a pdf that matches my composer?
Thank you!
Best Windows 7 - QGis 1.8.0 Lisboa
Answer
Bap, right clicking on a QGIS composer map item only locks its position, not the layers displayed within the map item. To lock the layers displayed, you must enable the "Lock layers for map item" check box in the map item properties panel:
I've been let loose in the workplace to learn python to do things in Arcmap 10. so, I am learning python as I go and trying to remember the programming I have done.
Where I am in this project is converting a dbf, or csv, xls in a simple fashion. from there, all the files will be copied together into one file. I've got the all-in-one xls working, but I can't find an easy, simple dbf to xls solution.
I condensed code found here: http://blog.gmane.org/gmane.comp.python.education/page=12
into:
from xlwt import Workbook
import dbfpy.dbf
input0 = '...file.dbf'
output = '...file.xls'
def test1():
dbf = dbfpy.dbf.Dbf(input0)
book = Workbook()
sheet1 = book.add_sheet('Sheet 1')
for row in range(len(dbf)):
for col in range(2):#chop to two for my purposes, gm
sheet1.row(row).write(col, dbf[row][col])
book.save(output)
test1()
This works, minus the lack of field names.
Answer
Like whuber says, you have to write out the headers explicitly. I loaded up dbfpy and xlwt in a virtualenv and ran this:
from xlwt import Workbook, easyxf
import dbfpy.dbf
from time import time
def test1():
dbf = dbfpy.dbf.Dbf("pipelines.dbf", readOnly = True)
header_style = easyxf('font: name Arial, bold True, height 200;')
book = Workbook()
sheet1 = book.add_sheet('Sheet 1')
for (i, name) in enumerate(dbf.fieldNames):
sheet1.write(0, i, name, header_style)
for (i, thecol) in enumerate(dbf.fieldDefs):
name, thetype, thelen, thedec = str(thecol).split()
colwidth = max(len(name), int(thelen))
sheet1.col(i).width = colwidth * 310
for row in range(1,len(dbf)):
for col in range(len(dbf.fieldNames)):
sheet1.row(row).write(col, dbf[row][col])
book.save("pipelines-xls.xls")
if __name__ == "__main__":
start = time()
test1()
end = time()
print end - start
This gives me headers in my xls:
I'm a fairly novice GIS user and know just enough to get in trouble in a lot of cases. I am trying to look at the connection between several measured parameters at a few hundred stations along several rivers and the landuse patterns in the areas drained in ArcGIS 10. I have created basic watershed analysis layers (flow directions, flow accumulation, catchments, etc) for the entire area of interest. Is there a way to easily sum up the land usage for the area drained by each station (i.e. get the number of cells of each type of land use along with the total cells drained) so that I could use that for further analysis? If so, how would be the easiest way to implement that so I don't have to do it individually for each station? I would assume there would be a way to script it, but I don't know anything about scripting in ArcGIS 10. Thanks so much!
I'm trying to create a Leaflet/Mapbox webmap with line layer (GeoJSON or KML) with some IDs and data coming from the .csv (after joining through this ID) file located on the same path as the index.html. Joined attributes need to show on lines popup windows.
I've found THIS article with exactly the same assumptions as mine:
But honestly I've fallen during recreation this example for my needs. It seems to be more complicated than just joining the .csv with lines.
Here is my zipped example website (clean, without any joining implemented, just .csv and kml variables)
The working example of map from mentioned article is HERE.
My clean code (with no joining yet) looks like this:
Projekty Budżetu Partycypacyjnego Dzielnicy Bielany
Join method from other article is described here: Joining CSV to GeoJSON in Leaflet?
For one of my agricultural projects, I really need reasonably accurate field boundaries for fields in the US, specifically Iowa right now but later, other similar states.
The USDA/NRSC/FSA Common Land Unit (CLU) data is exactly what I need, but apparently since 2008 it cannot be accessed by the public.
Does anyone know of a similar source, commercial or otherwise, of this type of data? I don't need much for attributes, just field boundaries that will let me aggregate analytic data within a field, accurately.
I used to be able to iterate over selected features using the following script in v2.0.
import processing
features = processing.getfeatures(layer)
for feature in features:
#Do whatever you need with the feature
However, this doesn't seem to work when I upgraded my QGIS to v2.2. Is there any new way on how to iterate selected vectors?
Answer
You can use the solutions given in Using processing algorithms from the console:
But if you look at what's in the module (version 2.2):
import processing
dir(processing)
no more getfeatures() or getFeatures()
You can control this with a little function adapted from Script de Python para filtrar por patrón de texto los métodos de Clases en PyQGIS de José Guerrero
import re
def get_patt(keyword, L):
return [item for item in dir(L) if re.search(keyword,item)]
get_patt("getfeatures",processing)
[]
get_patt("getFeatures",processing)
[]
but there is:
get_patt("features",processing)
['features']
So the command is:
features = processing.features(layer)
But for me, it is easiest with:
without processing:
layer = qgis.utils.iface.activeLayer()
with processing:
layer = processing.getObject("name_of_the_layer")
All the features:
for feat in layer.getFeatures():
geom= feat.geometry()
attr =feat.attributes()
Selected features only:
for feat in layer.selectedFeatures():
geom= feat.geometry()
attr =feat.attributes()
and, for example, an "universal solution":
if layer.selectedFeatureCount():
geom = [feat.geometry() for feat in layer.selectedFeatures()]
else:
geom = [feat.geometry() for feat in layer.getFeatures()]
I have over 100 features in a shapefile and I want to create individual viewsheds for each and then run processes on each of those resulting viewsheds in ArcGIS using Python.
Is it possible to write a script using a cursor to run through each "row" and use that individual observer point to run the viewshed?
Then update the cursor to the next row and rerun it until all features have been exhausted?
Or would I have to first Split the shapefile into over 100 files and run each individually?
Also, is it possible to designate the filenames for the resulting viewsheds (or shapefiles) using a specific field?
Answer
I've never used the viewshed tool so I can't speak to the specifics of using that tool, however in regards to your suggestion in this question I will add an answer here to build on KHibma's solution.
The first thing we need to address is the format of your SQL query: "[OID] = count". As it stands, this query will fail since "count", as represented here, is a string and strings need to be surrounded by single quotes: "[OID] = 'count'". Now, the second reason this query will fail is that OID is an integer field and as such there will be no OID value equal to "count". This is where string formatting comes into play: '"OID" = {}'.format(i) or '"OID" = {0}'.format(i) for those using Python 2.6.x (Use quotation marks to enclose the fieldname when using a shapefile. See ArcGIS help.).
To avoid certain pitfalls of making queries on the OBJECTID field, we can nest the loop into a cursor to make sure that we are only querying OID values that exist:
arcpy.MakeFeatureLayer_management ("C:/data/pts.shp", "pts")
with arcpy.da.SearchCursor('pts',['OID']) as cursor:
for row in cursor:
oid = row.OID # or whatever your OBJECTID field is called
arcpy.SelectLayerByAttribute_management ("pts", "NEW_SELECTION", '"OID" = {}'.format(oid))
outViewshed = Viewshed("elevation","pts",2,"CURVED_EARTH",0.15)
outViewshed.save("C:/sapyexamples/output/outvwshd"+oid)
I have a .NET application dependent on ArcGIS 9.3.1 which I am upgrading to ArcGIS 10.
I have an AxMapControl sat in a panel in a split container on a tab. The Dockstyle is set to Fill.
Since upgrading to ArcGIS 10 from 9.3.1, the painting of the map inside the map control has become inconsistent.
The initial size of the map as painted is 'wrong' and the map ends up only filling, say, half the panel, and a second vertical scrollbar appears beside the map. This second scrollbar is a glitch; you cannot use it, and it is eventually painted over. The issue with the the map seems to be painting only: all other information presented (for example in the legend) is as if the map were painted correctly across the whole panel.
After a few operations (such as zooming in/out etc.) the map painting seems to 'catch up' and the map is painted correctly, filling the whole panel. This 'catching up' seems to get faster each time I run the application which leads me to suspect some sort of timing or caching issue.
There is no consistency about when the problem manifests. Sometimes adjusting the screen resolution affects it, other times not.
This seems to be a problem internal to the control. Does anyone have a workaround (other than reverting to the 9.3.1 DLLS)?
UPDATE: I have tried the following without success:
Answer
The IScreenDisplay interface has a ScaleContents property; set this to true
As this question reveals, ESRI changed the default behaviour of the map control at ArcGIS 10.
Less redrawing in data view
In previous versions of ArcMap, if you changed the size of the ArcMap display while working in data view, either by resizing the ArcMap window or by docking/undocking/resizing a dockable window, by default your map was completely redrawn to fit inside the available display area. So the scale changed and the extent stayed the same (although you may have seen some extra geographic coverage based on how well the extent fit inside the new shape of the display area). At version 10, the default drawing behavior in data view has been changed so that when the display size is changed, your map is no longer completely redrawn to fit inside the display. Instead, the scale stays the same and the extent will change. If you make the display bigger, you'll see a larger geographic extent, and vice versa.
This has the performance advantage that the portion of the display unaffected by resizing doesn't need to be redrawn at all. For example, if you close a docked window, only the portion of the display that was obscured by the window needs to be redrawn. It is also easier to work with the display because geographic features on your map don't move around as you dock and undock windows. Features on your map remain in the same location in the display until you manually pan or zoom the map.
You can reverse this change by using the IScreenDisplay interface to set the ScaleContents property to true. This can be set once the map is loaded into the control.
For me, the line of code would read:
mapControl.ActiveView.ScreenDisplay.ScaleContents = True;
This restores the old map redrawing behaviour - and as a side-effect, cures the redrawing problem explained here!
I have not yet pulled this down to a code snippet that has just enough in it to reproduce the error below that I am getting when I run my code in IDLE.
pythonw.exe has stopped working
Windows is checking for a solution to the problem
However, I am seeing it when my script repeatedly runs Trace Geometric Network which I am doing to try and identify isolation valves for almost 1,000 locations on a water network.
These are a few lines from the code:
print_message("Trying to trace")
arcpy.TraceGeometricNetwork_management(
r"C:\Users\101791\test.gdb\WaterNetwork\WaterNetwork",
r"in_memory\WaterDistribution_Net",
r"C:\Users\101791\test.gdb\Flag",
"FIND_CONNECTED", "", "", "", "", "wSystemValve",
"TRACE_ENDS", "", "", "", "AS_IS", "", "", "", "AS_IS")
print_message("Traced OK")
and this is the ouput I get. The error in the image above occurs after it has run Trace Geometric Network eight times successfully immediately after printing "Trying to trace" the ninth time.
1
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'RV71818']
2
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'RV71818']
3
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'RV130266']
4
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'RV12053', u'RV415348']
5
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'RV12052', u'RV12051']
6
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'RV12060', u'RV12059']
7
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'RV12060', u'RV12061']
8
2074
Trying to trace
Traced OK
Isolation Valves for KCA 2074: [u'WSF405337']
9
2086
Trying to trace
If I try to start the process on the ninth flag pythonw.exe stops immediately. If I set my code to skip trying the 9-10th and 12-14th flags then all the others work until it gets to the 28th. The first failure occurs when the KCA value changes for the first time, but the 11th flag (which works) has the same KCA value of 2091 as the 12th-14th (which work) so particular KCA values do not appear to be the cause.
Does anyone have any thoughts on how I may be able to troubleshoot this?
I have a suspicion there may be some interplay between how much memory is needed to perform a geometric network trace and the way the Trace Geometric Network tool writes its result as an in_memory group layer.
I am using ArcGIS Desktop 10.2.1.01 on Windows 7 SP1.
I created a layer that contains points and I added it to a map as WMS layer in OpenLayers. I did a CQL filter and it works but I'd like to zoom to the selected point when I worked with getextent or getdataextent it doesn't zoom to the right place. It missed with a few Km. Here is my the code of filter:
function filterNom(){
var filterType ="cql";
var filter = document.getElementById('batima').value;
// by default, reset all filters
var filterParams = {cql_filter: null};
if (OpenLayers.String.trim(filter) != "") {
if (filterType == "cql")
filterParams["cql_filter"] = "batiment='"+filter+"'";
}
// merge the new filter definitions
wfs_layer.mergeNewParams(filterParams);
var bounds = new OpenLayers.Bounds
bounds = wfs_layer.getExtent()
map.setCenter(wfs_layer.getExtent(),18);
var center = map.getCenter();
center.transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));
alert(center.toString());
}
and here is the code to add the WMS layer:
wfs_layer = new OpenLayers.Layer.WMS("Patrimoines",
"http:// localhost/geoserver/cite/wms",
{layers: 'cite:patrimoine' , transparent: true },
{isBaseLayer: false, opacity: 1, singleTile: true, visibility: true,
projection: 'EPSG:900913'});
Answer
Your wms layer is a image, so the zooming to it will not work. I set up a vector layer with a protocol and filter and then am able to zoom to the particular feature. I did some cutting and pasting but I think I got the important parts. Let me know if this is sufficient
filter = new OpenLayers.Filter.Comparison({
type: OpenLayers.Filter.Comparison.EQUAL_TO,
property: "field name",
value: 123
});
protocol = {
version: "1.1.0",
srsName: "EPSG:900913",
url: your url,
featurePrefix: "prefix",
featureNS: namespace,
featureType: featuretype,
geometryName: "geom"
};
vector_LAYER = new OpenLayers.Layer.Vector("vector", {
strategies: [new OpenLayers.Strategy.Fixed()],
styleMap: styleMap,
filter: filter,
protocol: protocol
});
vector_LAYER.events.on({
featureadded: function() {
map.zoomToExtent(vector_LAYER.getDataExtent());
}
});
Then somewhere down down in my code I set the filter and refresh the layer.
filter.value = filter;
vector_LAYER.refresh();
How can I convert an excel file with x, y columns to a point shapefile?
There are some somewhat optional requirements in addition to the correct creation of a shapefile:
When using ArcGIS Desktop, I often use the Identify tool by clicking and dragging the tool over my data and identifying the features within it.
Is there a way to do this in QGIS or is setting the search radius within the map tools options the best that can be achieved for searching across multiple datasets?
In the ogr2ogr documentation it says that when using -a_srs srs_def:
Srs_def can be a full WKT definition (hard to escape properly), or a well known definition (ie. EPSG:4326) or a file with a WKT definition.
Where are the default well-known definitions listed? Do they include ESRI as well as EPSG? If not, why not, and can I add all the ESRI definitions listed at spatialreference.org? Is there an easy way to load an srs definition using an SRID from a PostGIS spatial_ref_sys table, aside from just running some sql to get the text?
Thanks!
Answer
Well known definitions are in the gdal/data directory. You can browse the current directory source online: http://trac.osgeo.org/gdal/browser/trunk/gdal/data. Look at the gcs.csv and pcs.csv files.
In the QGIS (2.99) Python Console, the name/variable qgis seems to be automatically defined.
However, in writing a plugin, there is an error whenever I refer to qgis. For example, in I have qgis.iface.mapCanvas() in my Plugin's `initGui(self)' function, I get an error
NameError: name `qgis` is not defined
How can fix this error?
By the way, I am trying to put the following code from this GIS.SE answer(https://gis.stackexchange.com/a/45105) into my plugin
tool = PointTool(qgis.iface.mapCanvas())
qgis.iface.mapCanvas().setMapTool(tool)
The code works fine inside QGIS main window (in Python Console), but generates the above error when moving it into the Plugin. I can't figure out how to let the PointTool respond inside a plugin. The code for PointTool class from the linked answer is:
class PointTool(QgsMapTool):
def __init__(self, canvas):
QgsMapTool.__init__(self, canvas)
self.canvas = canvas
def canvasPressEvent(self, event):
pass
...
I try to extract coordinates from gpx file and it returns empty list. I used xml.etree and have any idea why it is empty, because d.attrib return proper values.
import xml.etree.ElementTree as ET
f = ('St_Louis_Zoo_sample.gpx')
p = ET.parse(f)
root = p.getroot()
lats=[]
lons=[]
for d in root:
if d.tag == 'wpt':
y = d.attrib['lat']
x = d.attrib['lon']
lats.append(y)
lons.append(x)
Here it is a gpx file
St Louis Zoo sample
This self guided,GPS enabled tour of the world famous St. Louis Zoo, has 85 points of interest. Narratives in english,explaining each exhibit and provides guidance to zoo facilities.This audio tour guide can enhance your next visit.
wizardone, using GeoTours
St Louis Zoo sample
St Louis Zoo sample
Audio tour guide
St.Louis Mo.
Zoo
Forest Park
Animals
Asian Elephant
elephant
Waypoint
15.24
SymbolAndName
Bactrian Camel
camel
Waypoint
15.24
SymbolAndName
Black Rhinoceros
black rhino
Waypoint
15.24
SymbolAndName
Cafe Restroom
cafe restroom
Waypoint
15.24
SymbolAndName
Polar Bear
polar bear
Waypoint
15.24
SymbolAndName
California Sea Lion
sea lion
Waypoint
15.24
SymbolAndName
Cheetah
cheetah
Waypoint
15.24
SymbolAndName
Grizzley Bear
grizzly bear
Waypoint
15.24
SymbolAndName
Jaguar
jaguar
Waypoint
15.24
SymbolAndName
Lion
African lions live in a number of different habitats: grassy plains, open woodlands, semi-desert areas, even high mountains.
Waypoint
15.24
SymbolAndName
Answer
GPX files use XML namespaces (look at Parsing XML with namespace in Python via 'ElementTree' or Read GPX using Python ElementTree.register_namespace? for example).
The namespace is {http://www.topografix.com/GPX/1/1}
import xml.etree.ElementTree as ET
tree = ET.parse("St_Louis_Zoo_sample.gpx")
for elem in tree.findall("{http://www.topografix.com/GPX/1/1}wpt"):
print elem, elem.attrib['lon'], elem.attrib['lat']
-90.29408 38.63473
-90.28679 38.63368
-90.29323 38.63408
-90.29019 38.63533
-90.28976 38.63677
-90.28948 38.63496
-90.29458 38.63421
-90.29083 38.63633
-90.28715 38.63395
-90.28769 38.63347
With ElementTree (xml.etree), you can define the namespace before
namespace = {"gpx": "http://www.topografix.com/GPX/1/1"}
for elem in tree.findall('gpx:wpt', namespace):
print elem.attrib['lon'], elem.attrib['lat']
-90.29408 38.63473
-90.28679 38.63368
-90.29323 38.63408
-90.29019 38.63533
-90.28976 38.63677
-90.28948 38.63496
-90.29458 38.63421
-90.29083 38.63633
-90.28715 38.63395
-90.28769 38.63347
As with lxml
from lxml import etree
NSMAP = {"gpx": "http://www.topografix.com/GPX/1/1"}
tree = etree.parse("St_Louis_Zoo_sample.gpx")
for elem in tree.findall("gpx:wpt", namespaces=NSMAP):
print elem.attrib['lon'], elem.attrib['lat']
-90.29408 38.63473
-90.28679 38.63368
-90.29323 38.63408
-90.29019 38.63533
....
I'm trying to make a supervised classification in Grass 6.4.3 under Linuxmint 16.
After digitizing the training areas, I transformed them from vector to raster. Then I've done the following: just in case: g.region rast=train_area_2013_03_05_raster@PERMANENT,2013_03_05_toar4@PERMANENT
then
i.group group=group subgroup=sub_group input=2013_03_05_toar4@PERMANENT,2013_03_05_toar5@PERMANENT,2013_03_05_toar3@PERMANENT
and
i.gensig trainingmap=train_area_2013_03_05_raster@PERMANENT group=group@PERMANENT subgroup=sub_group signaturefile=superv_class
but the resulting signature file, only contains one character: "#".
The output of i.gensig in the console is the following:
Finding training classes...
3 classes found
Calculating class means...
Calculating class covariance matrices...
Signature 1 not invertible
Signature 2 not invertible
Signature 3 not invertible
i.gensig complete.
Obviously if I try: i.maxlik group=group@PERMANENT subgroup=sub_group sigfile=superv_class class=results_class
I obtain the following error message:
ERROR: Unable to read signature file
I've been searching, but I couldn't find the solution.
I've tried: this Problem running i.gensigset in Grass. Any ideas?
My problem is similar to this: i.maxlik cannot read i.class output signature file. But here, there's no solution suggested. As in this case, I also used i.cluster to generate an unsupervised signature file and i.maxlik reads this perfectly well.
Here there's another similar question, but again, with no solution...
What does Signature X not invertible means?
Any idea? Thanks for your help!
If I save layer files in ArcGIS 10, can they be opened and modified in ArcGIS 9.3?
I am a beginner using QGIS (1.8.0 Lisboa) on a Mac (OS X 10.6.8 Snow Leopard). I work on a file with two layers in the following sequence (top to bottom):
Sometimes when changing the scale (scrolling in and out) the Google layer updates a little wrong in size or position compared to the polygon layer. But as soon as I "pan" the map everything reconstructs right. Now the problem I have concerns the print composer.
The result of any kind of export - pdf, or image or svg - of my map (vector layer + Google layer in the background) differs from what is shown in the preview of the print composer: the vector layer and the Google layer positions are mismatching and in a different scale. I can send some pictures showing screenshots of the print composer preview and the final results if required.
I tried almost everything:
The problem appears only using an open layer map (Google, OpenStreetMap...). If I try the same using two vector layer (shapefile) - one with the buildings, another one with the city in background - everything works out great.
...exasperating! Since it is a real important research project I would be very grateful for some help and tips.
I have a shapefile with a feature that is non-contiguous. I want to split the areas on the left from the area on the right. However, QGIS's Split Features tool requires you to draw a line over the feature. When I draw a line between the two areas of this feature, it doesn't do anything.
How should I split this feature?
Here's a screen shot. The feature I want to split is in yellow.
Answer
You could try Multipart to Single Parts in the Vector toolset, which should split what you have there, a multipolyon, into its constituent polygons, from where you should be able to delete the polygons you need to.
You can always recombine them later if you need to.
I have a polygon layer which has a number attribute and a text attribute. When I make a legend in the composer, it automatically chooses the number to display while I would like to have both the number and the text.
How do I display two attributes in a legend? If not possible, how do I select which attribute is displayed in the legend?
I am not looking to have a multi-attribute label in the main window, I wish to have a multi-attribute label in the legend of the composer, or choose which attribute is displayed in the legend.
QGIS 2.8 Win8 64 bit.
Answer
You could simply concatenate the three classification fields in the layer's properties.
"Unit_Num" || ' ' || "Age" || ' ' || "Name"
Then it will also appear in the legend in the Print composer:
Unfortunately, the concatenated text will also appear in the main window.
Edit: The legend is derived from your classification, but you can of course label your features with another field:
I am in search of an algorithm to solve the first/direct geodetic problem precisely in 3D. The input is:
The output should be the geodetic position (again, WGS84) of point A.
I have read Algorithm for offsetting a latitude/longitude by some amount of meters and this is how I've computed the solution in the past. It is obviously computationally efficient, but not precise enough anymore.
Currently, my idea is to
Of course, I still would like things to be fast, but precision has a definite priority here.
My questions:
What do you think of this approach? Does it work, is it precise? (I fear that compensating for pitch/roll not being relative to the center of the earth would probably imply a LOT of additional complexity).
Is PROJ.4's pj_transform() capable of computing the conversions in steps 1 and 4? If not, what C/C++ library can you recommend? This must be a common routine in some library?!
Note to @whuber who referred to 30448! I think I'm slowly starting to resolve misunderstandings, getting an idea of what you do. Questions about your answer to that question:
When you say "geocentric", you mean cartesian coordinates from earth center (=ECEF), right? (I first interpreted "geocentric coordinates" just like "geodetic coordinates", except the latitude being determined using earth-center instead of ellipsoid-normal. Now I know that "geocentric latitude" != "geocentic coordinates". Whew.)
In the first image of 30448, you show "the geocentric frame at the North Pole". The origin of that coordinate system really is at the center of the earth, right? But if your first image started at earth's center, simply rotating wouldn't get you to the origin of the local coordinate system. So why does it start at (0,0,6370) in the first picture? You could have also chosen (6370,0,0), but then rotations would be different. Is it that each initial offset corresponds to an Euler convention (Z-Y'-Z'')?
If so, it then seems to me that my drafted approach is very much similar to your solution at 30448. Are there other approaches that are more precise? Given I compute using doubles, what accuracy do you expect?
Slightly OT: At the end, you add a local Z displacement motivated by differences against MSL - but I thought the WGS84 datum you referenced earlier is not a tidal/MSL datum, but an ellipsoidal. Assuming you sticked to WGS84, you'd add the local Z offset based on the distance to the ellipsoid, right?
@whuber, @martin f provided quick replies and patience - I read the whole day, trying to understand more.
I have a layer (SQLlite) of city lots and would like to accomplish the following:
For each city lot,
I am stumped, as spatial join will only aggregate data.
Also, I don't have administrator privileges on my computer so cannot install plugins.
(Using QGIS 2.6)
Answer
Concerning your first request (concatenate the list of unique lot id's which are bound within a radius of 30m from the given lot), I guess you can do that with a SQL query in QGIS :
SELECT b.id as ID_ref, group_concat(a.id) AS ID_within_30m
FROM city_lots a, city_lots b
WHERE b.id = 1 AND b.id != a.id
AND ST_Distance(a.GEOM, b.GEOM) < 30;
Which should output something like:
ID_ref | ID_within_30m | -------|---------------| 1 | 5,6,8 |
(There is probably more efficient ways to do this if you are handling a large dataset)
Concerning your second request, if it is about doing the same analysis for each lot, I guess you can do that with a JOIN (and a GROUP BY) in SQL (and using a spatial index):
CREATE TABLE t (startrange INTEGER not null, endrange INTEGER not null);
INSERT INTO t VALUES (1, (SELECT count(rowid) FROM city_lots));
SELECT b.id as ID_ref, group_concat(a.id) AS ID_within_30m
FROM city_lots a,
t JOIN city_lots b on (b.rowid >= t.startrange and b.rowid <= t.endrange)
WHERE a.rowid IN (
SELECT rowid FROM SpatialIndex
WHERE f_table_name = 'city_lots'
AND search_frame = ST_Buffer(b.GEOM, 30))
-- AND b.id != a.id
AND ST_Distance(b.GEOM, a.GEOM) < 30
GROUP BY b.ID;
Which should output the entire result :ID_ref | ID_within_30m | -------|-----------------| 1 | 1,5,6,8 | 2 | 2,7,9 | 3 | 3 | .....
(the clause b.id != a.id is commented in order to get a row in ID_ref even if there isn't any other lots in the 30m radius)
I want to check whether there is a way to update the vertices x, y value and remove all values from the 3rd decimal. I only want 2 decimals (no matter the 3rd decimal is 0 or 9, just remove it).
The original data and the desired data is like the following:
25642.882 38471.416 (need to change to 25642.88 38471.41)
25648.078 38446.104 (need to change to 25648.07 38446.10)
25655.240 38406.615 (need to change to 25655.24 38406.61)
Any advice?
I have a postgis table of simple linestrings, each one is a straight line with just two points. I want to select every linestring, but extended from both sides by 1 meter. So i have lines, that are longer by 2 meters. How can this be done?
Answer
I seems very similar to post "postgis, extrapolate a line". If I avoid repetition of cited post, I think you just need to extrapolate beyond your extreme points. In a query you get something like this should work:
SELECT ST_MakeLine(ST_TRANSLATE(a, sin(az1) * len, cos(az1) *
len),ST_TRANSLATE(b,sin(az2) * len, cos(az2) * len))
FROM (
SELECT a, b, ST_Azimuth(a,b) AS az1, ST_Azimuth(b, a) AS az2, ST_Distance(a,b) + 1 AS len
FROM (
SELECT ST_StartPoint(the_geom) AS a, ST_EndPoint(the_geom) AS b
FROM ST_MakeLine(ST_MakePoint(1,2), ST_MakePoint(3,4)) AS the_geom
) AS sub
) AS sub2
This question is very similar to How to define layer's CRS and avoid the CRS dialog in PyQGIS?, but it's doesn't work for me.
This is my try:
vlayer.setCrs(QgsCoordinateReferenceSystem(25830))
iface.mapCanvas().mapRenderer().setDestinationCrs(QgsCoordinateReferenceSystem(25830))
But the dialog pop up asking for CRS to choose, and when I ask for the CRS of the layer, it show a diferent CRS:
print "CRS: " + vlayer.crs().geographicCRSAuthId()
It show CRS: EPSG:4258 and It must be 25830
I am trying to add an interactive map HERE to my ArcMap? I am using ArcGIS 10 and the path I have been following so far is through "Add WMS Server" (Catalog -> GIS Servers -> Add WMS Server). Unfortunately, this is not working for me. The way I see it, this can have two reasons:
In case of the former, what is that I am doing wrong?
I am using the AppID and AppCode from my HERE developer account. I have tried different things for the URL. The last one I tried was simply https://www.here.com/?map=51.54916,5.64706,10,normal (I was running out of ideas). Usually it tells me there is a connection error.
In case of option 2, that this isn't the way to go, what should I do to display HERE maps in ArcMap? It surely isn't offered through Add Basemap/ArcGIS Online.
In QGIS it is possible to add HERE maps in multiple ways. See this post: HERE background maps in QGIS possible? I hope there is a way around to get it fix in ArcMap as well.
I think that there is a third option, and that it isn't possible at all, period. But I sure hope somebody can tell me how I can get it done.
I'm looking at putting together a light Python install on a thumb drive that I can take to client offices and do some basic GIS processing. I'm thinking of installing the following:
Does anyone have any suggestions of any libraries / software that might be a useful addition? Typically I'll be working with raster datasets, but any and all suggestions would be appreciated.
(I'll cw the question if people think that its appropriate).
Answer
Portable Python http://www.portablepython.com/
Interesting project you have will it be available to the GIS community?
Maybe of interest Portable GIS http://www.archaeogeek.com/blog/portable-gis/
In QGIS, I have a polygon layer with square polygons in it. The polyogons are of same size but different orientation. I want to create a set of maps using the atlas feature. My aim is to rotate the map accordingly to the orientation of the polyogns.
The "item properties" offer the possibility to insert an expression for the map rotation. If I had a rotation-field in the property table of the poygons layer I'd be able to automate map rotation accordingly, I suppose. So my question is: is it possible to calculate the rotation of quadratic polyon features relative to the N-S (or W-E) axis?
Here I found a similar question; however, not being familiar with PostGIS, I want to accomplish this in QGIS.
Answer
You can use $x_at(n) to get the n-th x coordinate. Then you can either use atan2 with the coordinate differences:
(180/pi())* atan2($y_at(1)-$y_at(0), $x_at(1)-$x_at(0))
or azimuth by making points:
180/pi() * azimuth(make_point($x_at(0), $y_at(0)),
make_point($x_at(1),$y_at(1)))
The 180/pi() gets you degrees. I think the only difference between these two is the angle for zero, and maybe the direction (clockwise or anti-clockwise...).
I am running ArcGIS 10.2.1 on Win7 64bit. I have finished doing an NDVI, and after defining a vegetation presence threshold at 0.05, opened the Raster Calculator (via Map Algebra in the Arc Toolbox) to generate a new raster with a 0 value for no vegetation and a 1 value for vegetation. I'm getting the following error:
Executing: RasterCalculator ""NDVI_ik_06.tif" > 0.05"
ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
File "", line 1, in
File "", line 2, in rcexec
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\__init__.py", line 24, in
from arcpy.toolbox import *
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\toolbox.py", line 356, in
from management import Graph, GraphTemplate
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 22, in
import _management
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\_management.py", line 14, in
import _graph
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\_graph.py", line 27, in
import numpy
ImportError: No module named numpy
Failed to execute (RasterCalculator).
Failed at Tue Apr 07 11:24:48 2015 (Elapsed Time: 0.09 seconds)
Any ideas?
How do i create a linestring from points sort by date in PostgreSQL?
I tried:
SELECT observations.id,
ST_MakeLine(observations.geom ORDER BY observations.date) AS newgeom
FROM observations
GROUP BY
observations.id;
This is going wrong. It returns me a linestring with same coordinates. So, instead of a line from a to b, it gives me a to a.
What am I doing wrong?
Answer
The point is, if you GROUP BY the row id column, you will get one result row per input row (this is equal to grouping by the actual date column)! And since a Linestring is only valid with a minimum of two points, ST_MakeLine adds the same point twice.
Either run
SELECT ST_MakeLine(geom ORDER BY date) AS geom
FROM observations
;
to get one line for all points, or
SELECT ROW_NUMBER() OVER() AS id,
*
FROM (
SELECT ST_MakeLine(geom, LEAD(geom) OVER(ORDER BY date)) AS geom
FROM observations
) q
WHERE geom IS NOT NULL
;
to get a linestring between each pair of consecutive points.
How can I add an external online map in ArcMap by inserting the http address, to use as a basemap in order to georeference multiple aerial photos?
I don't really care about the precision of the georeferencing. I just want to use the template as a quick-to-use index for the aerial photos.
I have to do a multi-criteria analysis to answer the question : "which is the best lot to develop".
A few of the criterias are :
I thought I'd give it a try with QGIS, and here's how I've done :
add the following columns in my lots layer attributes table :
Convert my lots layer to points using "polygons to centroids"
Run the "distance matrix" tool
Open the CSV to run an operation in excel (bus stop grade is 1.0 if nearer than 200m, and 0.0 if more than 750m, but I coulnt find the MIN() function in QGIS)
Join the resulting CSV back in QGIS
Repeat the same for shops
Run the "point in polygon" tool to select all the points in nature-protection area
Set 0.0 to all selected points
Repeat for other "in ... area" criterias
Run the "spatial join" tool to merge flood danger area information
Run a calculation using the column calculator to have the mean grade (using determined factors for each criteria)
Once all that done, re add the BUILDING LOTS shapefile once for each criteria
For each criteria, join the converted layer (that one with the centroids) on the LOT id
Set the display to a gradient from red to green according to the corresponding criteria attribute and the mean grade attribute
Now, after a good 2 days of work, i now have all my criterias displaying in green if good choice for building, and red if bad choice, and I have my synthesis which aggregates all my criterias in one beautiful red-green map. (and I also have a huge mess in my "shapefiles" folder)
Now the problem.
What if :
Do I have to start all over again ?
Am I correctly using the wrong tool, or am I using the correct tool wrong ?
Would it be easier with a commercial GIS software ?
I see what answerers/commenters mean, and i didn't really think of using rasters.
However, the main question was more about the ability to try different scenarios or update the base data without having to restart all the process from scratch.
It seems that your suggestions are not much more flexible than what I suggested (even maybe more complex) since you have news steps : - (for each criteria) rasterization. - (in the end) sampling (quite complex if you want to include partial overlaps)
That Sextante Model builder seems awesome; in fact I was exactly thinking at something like that when posting my last comment.
I've used Grasshopper3D quite alot (it has nothing to do with GIS software) which is a great plugin for the Rhino3D modeler and which uses the same concept of node graph workflow construction. (example : http://designreform.net/2009/07/rhino-grasshopper-parametric-truss )
This seems so well adapted to a lot of GIS data analysis that I'd love to see a GIS software really built around such a node graph tool.
I'm looking forward to try Sextante Modeler and let you know how it worked out. I wish i had found about it by myself by googling it, but i didn't know the keyword "model builder".
Using standalone script I loaded a layer, and now I need to export it as an image. See code below:
from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
QgsApplication.setPrefixPath("C:\OSGeo4W\apps\qgis", True)
qgs = QgsApplication([], True)
qgs.initQgis()
rasterpath = "E:/MODIS DATA/2016/15-10-2016/TIRUNELVELI.tif"
lyr = QgsRasterLayer(rasterpath, "TIRUNELVELI")
QgsMapLayerRegistry.instance().addMapLayer(lyr)
I checked the post QGis Save Raster as Rendered Image, but the code is executed using Python console. If I use the same code it shows an error:
NoneType object has no attribute 'clone'.
Updated Code after Comment
from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
rasterpath = "E:/MODIS/TIRUNELVELI.tif"
layer = QgsRasterLayer(rasterpath, "TIRUNELVELI")
QgsMapLayerRegistry.instance().addMapLayer(layer)
uri = "E:/MODIS/newStyleqgis.qml"
layer.loadNamedStyle(uri)
extent = layer.extent()
width, height = layer.width(), layer.height()
renderer = layer.renderer()
provider=layer.dataProvider()
crs = layer.crs().toWkt()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
file_writer = QgsRasterFileWriter('E:/MODIS/abcd.tif')
file_writer.writeRaster(pipe,
width,
height,
extent,
layer.crs())
updated but showing error as shown : 
I'm new to GIS.
My understanding of shapefiles is that they are (contain?) featureclasses, that contain features (points, polylines, polygons). Is there a clear hierarchy for these files. Would the next level after features be individual points? Are polygons and polylines merely point collections? What implication does this have for working with the data - e.g. can all geometries be worked with at the "point" level?
Is there a resource that would help me to understand the data structure and/or object model of these files? What I am essentially looking for is an explanation as to what exactly shapefiles are and what I should know when working with them, especially when using ArcObjects. In addition to this - how does this relate to other software like QGIS? Do they operate in the same way?
I am relatively new to ArcPy 10. and I want to return the list of feature class from a dataset. I can see them through the message, but cannot see them from the result window. Actually there are only inputs, environments and messages I can view from the results window. Am I missing something in the script? Thanks.
import arcpy
import json
from arcpy import env
env.workspace=r"C:\sde\pan18.sde"
featureDatasets = arcpy.ListDatasets("SDE.POINT_OF_INTEREST")
fd=featureDatasets[0];
fcList=arcpy.ListFeatureClasses("*", "ALL", fd)
pf=[]
for fc in fcList:
pf.append(fc)
arcpy.AddMessage("--------Encode the list to a json string-\n")
layerlistEncode=json.dumps({'resource':pf})
print layerlistEncode
arcpy.AddMessage("--------Encode the json string to a list\n")
layerlistDecode=json.loads(layerlistEncode)
print layerlistDecode
print layerlistDecode["resource"]
arcpy.SetParameterAsText(0, layerlistEncode)
arcpy.AddMessage(layerlistEncode)
arcpy.AddMessage("finished1")
Edit / Delete Edit Post Reply With Quote Reply With Quote Multi-Quote This Message Top Bottom
Answer
You need to set up a derived, output parameter for SetParameterAsText to have any effect. See the help on Setting script tool parameters.
In ArcMap, it is possible to automatically calculate the unique values for a raster and then apply a different style to each unique value (see 1st image).
However, in QGIS, I have to manually add values for styling when using the "Singleband pseudocolor" style (see 2nd image). Is there a way of auto-populating the unique values like ArcMap does?
ArcMap:
QGIS (how to auto-populate the value list with unique values?):
Answer
For those still looking for this. Unique raster values have been added to QGIS 3. "Added by Nyall Dawson about 1 year ago
[FEATURE] Allow classifying paletted renderer using unique values from a raster layer
Adds an easy way to style discrete rasters such as landuse classes using the Paletted renderer. Just select the Paletted renderer, pick a band, then hit the "Add Unique Values" button. The unique pixel values will be fetched from the layer and a color assigned to each using the currently selected color ramp."
I am very interested in the HTML5 canvas element for maps; work like Cartagen from an MIT Media Lab researcher looks very promising, for example. There is some interesting SVG based work at carto.net. WebSocket seems like a fantastic API for live geospatial data. I have done a few HTML5 experiments, for example here, with Flickr and Picasa data. What are people getting up to with these new technologies, or what have you tried?
I have just downloaded the newest version of QGIS (3.2) and the processing toolbox is nowhere to be found. It says that it is now apart of the core package, but i cannot find it on anywhere on the program.
I am trying to calculate the centroids of a polygon and add the x- and y-coordinate to their respective fields:
arcpy.env.workspace = r"C:/..."
input_fc = "shapefiles/precincts_v.shp"
arcpy.AddField_management(input_fc, "X", "DOUBLE")
arcpy.AddField_management(input_fc, "Y", "DOUBLE")
arcpy.CalculateField_management(input_fc, "X", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[0]", "PYTHON")
arcpy.CalculateField_management(input_fc, "Y", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[1]", "PYTHON")
but I an invalid syntax error. Where is the mistake?
Edit: exact error message
Traceback (most recent call last):
File "C:\Users\mquentel\Dropbox\PLSS\CA voting data\Python\CONSTRUCT_merge_voting_and_census.py", line 17, in
arcpy.CalculateField_management(input_fc, "X", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[0]", "PYTHON")
File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 3360, in CalculateField
raise e
ExecuteError: ERROR 000539: SyntaxError: invalid syntax (, line 1)
Failed to execute (CalculateField).
Answer
Replace your Field Calculator lines with
arcpy.CalculateField_management(input_fc, "X", "!SHAPE.CENTROID.X!", "PYTHON_9.3")
arcpy.CalculateField_management(input_fc, "Y", "!SHAPE.CENTROID.Y!", "PYTHON_9.3")
Works for me with Decimal Degrees.
I have a line in ArcMap 10.2, and want to create lines to the right of it at 100 m, 200 m, 300 m etc. up until 1.7 km. Obviously typing in these values under 'multiple ring buffers' will take too much time. Is there a way for ArcMap to compute this automatically across the equal intervals?
If you could answer through basic ArcMap methods that would be great. I do not understand the coding language that some of you use!
I have a feature class of 20k records. trying to import it into postgis using arcpy and psycogp2. I have a pretty basic script
cur.execute('''drop table if exists developed_lands.daysemtric_test''')
cur.execute('''CREATE TABLE developed_lands.daysemtric_test
(pid serial primary key,
pop_dns float,
shape geometry);''')
with arcpy.da.SearchCursor(shp, ["POP_DNS","SHAPE@WKT"]) as cur2:
for row in cur2:
sql = "insert into developed_lands.daysemtric_test(pop_dns,shape)values(%s,ST_GeomFromText(%s)"
cur.execute(sql,(row[0],row[1]))
x+=1
print x
conn.commit()
error:
Traceback (most recent call last):
File "C:\Users\rizagha\Desktop\module1.py", line 16, in
cur.execute(sql,(row[0],row[1]))
InternalError: geometry requires more points
HINT: "...6391.65285167098 737288.47085659206)," <-- parse error at position 13049795 within geometry
^
This script works because If I skip the first feature it inserts the records up to row 2,000 then throws the same error. I put the function st_makevalid() around ST_GeomFromText(%s) and it did nothing to help. running repair geometry on the feature class now, but that usually takes hours...
the first row where the error is happening, the geometry is extremely complex. is it possible the WkT is too long?
UPDATE
repair geometry finished and I reran the script and I got the same error
Answer
TL;DR;: PostGIS restricts some kinds of invalid WKT geometries that are otherwise allowed by some formats and software, such as shapefiles and ArcGIS. You can sneak these in via WKB.
Let's take an invalid geometry, described by WKT:
POLYGON((5 5, 20 30, 30 5, 5 5), (22 13, 16 13, 22 13))
This shape can be loaded into Esri and GDAL/OGR software.
# Esri
import arcpy
e = arcpy.FromWKT('POLYGON((5 5, 20 30, 30 5, 5 5), (22 13, 16 13, 22 13))')
print(e.JSON)
# {"rings":[[[30,5],[20,30],[5,5],[30,5]],[[16,13],[22,13],[16,13]]],"spatialReference":{"wkid":null}}
# GDAL/OGR
from osgeo import ogr
g = ogr.CreateGeometryFromWkt('POLYGON((5 5, 20 30, 30 5, 5 5), (22 13, 16 13, 22 13))')
print(g.IsValid()) # False
print(g.ExportToJson())
# { "type": "Polygon", "coordinates": [ [ [ 5.0, 5.0 ], [ 20.0, 30.0 ], [ 30.0, 5.0 ], [ 5.0, 5.0 ] ], [ [ 22.0, 13.0 ], [ 16.0, 13.0 ], [ 22.0, 13.0 ] ] ] }
But JTS, GEOS and PostGIS do not allow this kind of invalid geometry in WKT form:
SELECT 'POLYGON((5 5, 20 30, 30 5, 5 5), (22 13, 16 13, 22 13))'::geometry;
ERROR: geometry requires more points
LINE 1: SELECT 'POLYGON((5 5, 20 30, 30 5, 5 5), (22 13, 16 13, 22 1...
^
HINT: "...0, 30 5, 5 5), (22 13, 16 13, 22 13))" <-- parse error at position 55 within geometry
However, you can load the WKB equivalent into PostGIS without issue:
SELECT '010300000002000000040000000000000000001440000000000000144000000000000034400000000000003E400000000000003E400000000000001440000000000000144000000000000014400300000000000000000036400000000000002A4000000000000030400000000000002A4000000000000036400000000000002A40'::geometry AS g
INTO TEMP geoms;
SELECT ST_AsEWKT(g) FROM geoms;
st_asewkt
---------------------------------------------------
POLYGON((5 5,20 30,30 5,5 5),(22 13,16 13,22 13))
(1 row)
So you can modify your INSERT to use hex-encoded WKB instead of WKT. I'm not sure what Esri's e.WKB is (it's not ISO or EWKB), but you can use OGR for this: g.ExportToWkb().encode('hex')
And you should repair these geometries to avoid further complications:
SELECT ST_AsText(ST_CollectionExtract(ST_MakeValid(g), 3)) FROM geoms;
st_astext
--------------------------------------
MULTIPOLYGON(((5 5,20 30,30 5,5 5)))
(1 row)
(but use UPDATE instead).
Could someone write a brief albeit dumbed down idiot's guide to installing postgis and loading national Tiger data on ubuntu? I've tried a few guides, namely, http://wiki.bitnami.com/@api/deki/pages/302/pdf, but I'm not having much luck. I apologize for the open ended nature of this question.
Answer
Since you have PostGIS 2.1.1 you're ahead of the game. Make sure you have wget installed, it is what will download the data from the Census FTP site.
Create a gisdata directory with:
sudo mkdir /gisdata
Use chown and chgrp commands to change the ownership and group of /gisdata so that your normal user can read and write to /gisdata.
Start psql and connect to your database. Once in psql use
\a
and
\t
so that the results of the query are formatted correctly.
I forgot this part initially! Before you can use the loader script you need to do some house cleaning. First thing is to make sure the tiger schema is in your search path. Next up check the values in tiger.loader_platform and tiger.loader_variables. These two tables control variables for the loader script like your username and password. I usually just edit them in PGAdmin. Next up you'll need to run a script that populates lookup tables and other bits of background goodness the geocoder will need. First set an output file:
\o nation_generator.sh
then run:
SELECT loader_generate_nation_script('sh');
Then exit psql and run the file:
sh ./nation_generator.sh
Then hop back into psql and type:
\o loader_script.sh
to output the results of the query to a textfile called loader_script.sql. Then execute the function that generates the loader script:
SELECT loader_generate_script(ARRAY['DC','RI'], 'sh');
This is the query whose output will be redirected to loader_script.sql. Replace 'DC' and 'RI' with the two letter abbreviations of the states you want to download.
Exit out of psql and run the script with this command:
sh ./loader_script.sh
This will download the files for the state(s) you select, unzip them, and import the data into your PostGIS database.
