Using ArcMap 9.3, I want to georeference a scanned map of Sudan (Africa) as .jpg file format,coordinate points are marked on scanned map in 4 degree intervals. I am not familiar to georeferencing in ArcMap.
So looking for a detail guide line for georeferencing in ArcMap (using 4 control points) and also looking for a guide line to edit the spatial reference in base map with the suitable PCS information for the AOI.
I would like to know how to get all the raster values within a polygon using gdal or pygeoprocessing, without reading the entire grid as an array.
pygeoprocessing and gdal can do zonal statistics but only the min, max, mean, stdev or count are available from such a function. Since zonal statistics need to access the values, would it be easy to extract values the same way ?
I found a very similar question here : (Getting pixel value of GDAL raster under OGR point without NumPy?) but only for a particular "point".
In QGIS 2.4, i am unable to set the default value as NULL if i have not updated any value on date fields. Its automatically taking the current system date. I am expecting to set the default values(NULL) for date fields in properties of layer, since that default values should be displayed while capturing the new feature itself. It may be edited if required. Please refer the attached screen shots.
I am programming for an web editing application with ArcGIS API for JavaScript. I've published a version of ArcSDE with Feature Acces and I need to know how can I change in the web application from this version to another through the user that is connected to the application using the security manager of ArcGIS Server.
I've seen that exist a parameter gdbVersion and SetGDBVersion of featurelayer, but I don't know how to use. I need to see some example of this because the information in the documentation is not enough.
Can anyone help me? Any sample, please? Thanks
Answer
i definitely disagree with Sunil's comment above.
take a look at the following sample https://developers.arcgis.com/en/javascript/jssamples/widget_identitymanager_client_side.html
within the storeCredentials function you can see that the username of the person that signed in can be dug out using a statement like this
idObject.credentials[0].userId;
>>>"rick"
you would have to write your own logic to use information like this to determine how to set an appropriate corresponding SDE version on the featureLayer for editing.
i have shapefile that contains a lot of parcels "parcel has like fields (n° parcelle,proprietaire, cin) see picture :
and i want for all parcels to generate this plan :
take a parcel and generate this. (red text is clarification)
software or language programming not important, i'm open to any idea or answer that can solve this problem(product presentation format for all parcels with the format in picture 2 ).
I would like to know what is the sense of using st_buffer operation. In PostGIS reference guide there is only the syntax and a short description of this operation.
Could you give some examples to understand the advantages of this operation?
I have a geotiff file and I use geotools library for parsing it.
This is my code:
File f = new File("myFile.tif");
ParameterValue policy = AbstractGridFormat.OVERVIEW_POLICY.createValue();
policy.setValue(OverviewPolicy.IGNORE);
ParameterValue gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();
ParameterValue useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
useJaiRead.setValue(true);
GridCoverage2D image = new GeoTiffReader(f).read(new GeneralParameterValue[]{policy, gridsize, useJaiRead});
Rectangle2D bounds2D = image.getEnvelope2D().getBounds2D();
GridGeometry2D geometry = image.getGridGeometry();
System.out.println(geometry);
System.out.println(bounds2D);
The output of geometry is:
GridGeometry2D[
GeneralGridEnvelope[0..1199, 0..1399],
PARAM_MT[
"Affine",
PARAMETER["num_row", 3],
PARAMETER["num_col", 3],
PARAMETER["elt_0_0", 996.4799194335938],
PARAMETER["elt_0_2", -598386.1975402832],
PARAMETER["elt_1_1", -999.6777954101562],
PARAMETER["elt_1_2", 650290.4111022949]
]
]
The output of bounds2D is:
java.awt.geom.Rectangle2D$Double[x=-598884.4375, y=-748758.6635742188, w=1195775.9033203125, h=1399548.9135742188]
What does this information mean?
Two rasters
I want to erase the red raster from the slope raster 
I found this thread
https://community.esri.com/thread/190381-how-to-perform-reverse-clip
and they suggested something like this (for my case block is the red raster and buff is the slope)
Con(IsNull("block.tif"),"buff.tif")
this produces the output in the multi colors. But it is not the full erase, part of the buff/slope is missing
I do not want to go the vector conversion route.
Answer needs to use Raster Calculator syntax
Answer
Use the Reclassify tool to turn your red areas into NoData values and the red area NoData values to zeros. Before you do this ensure use the geoprocessing settings to ensure the output from the Reclassify has the same extent at both data sets. Then add the resulting reclassified surface with your slope using the Plus tool.
Con(IsNull("block.tif"), "buff.tif")
I use QGIS, version 2.18.9. I´ve several polygon-shapes with topological errors. Before I tried to apply the "intersection"-tool and got the message: Input layer A contains invalid geometries (feature 77). Unable to complete intersection algorithm. Firstly, I used the "check-validity"-tool, which detected some line-intersection and double nodes.
Then, I tried the "topology checker" with "must not have invalid geometries".
I´ve had the impression, that all errors where found then. So I applied the LWGEOM algorithms "Make valid" but got a error another messages ([error 87] falscher Parameter see log for more details).
I´m an absolute QGIS-beginner. So, I have not a single idea how to solve this problem. It has to be a simple automatically way, ´cause I´ve a big dataset here.
I installed some python libraries using pip (pandas, scipy, simplekml) and ever since every time I try to load QGIS I get this error:
Warning: loading of qgis translation failed
[/usr/share/qgis/i18n//qgis_en_US]
Warning: loading of qt translation failed
[/usr/share/qt4/translations/qt_en_US]
Warning: QCss::Parser - Failed to load file "/style.qss"
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
QH6248 qh_lib_check: Incorrect qhull library called. Caller uses reentrant Qhull while library is non-reentrant
QH6249 qh_lib_check: Incorrect qhull library called. Size of qhT for caller is 8184, but for library is 2896.
QH6255 qh_lib_check: Cannot continue. Library 'qhull 7.2.0 (2015.2 2016/01/18)' uses a dynamic qhT via qh_QHpointer (e.g., qhull_p.so)
What is wrong? How can I fix it ? I have version 2.14. Should I remove it? Or maybe remove qhull library ?
Answer
Got the same problem on Linux Mint and finally fixed it.
The problem is that you've installed a few python libraries and one of them is a dependency of QGIS.
QGIS ~2.18 works fine with Scipy 0.17, and after upgrade this library to the newest version you have 0.19.
QGIS won't load with 0.19 because libqhull error.
All you need to do is downgrade, simply:
pip install scipy==0.17
And everything should be ok. You can also try with newer version but I'm not sure which one is good. I downgraded to 0.11 and it works form me again.
My QGIS 2.8.2 standalone app provides several decorating layers (e.g. cartographic grid, background layers, reference points) which should be not changed by the user (neither editing nor hiding), as well as layers, the user can modify and query.
The app has a legend widget which shows all layers, since all layers are registered with map registry. The widget is populated with the layer tree as described in How to add a legend to a canvas in a standalone PyQGIS application?.
Instead of exposing all layers in the tree view, I want to restrict the tree view to the layer group, whose children represent the user mutable layers. Thus the user can change visibility and ordering of only these layers, and must keep the decoration as it is.
What I tried so far is to initiate QgsLayerTreeMapCanvasBridge with this group node (instead of layerTreeRoot), pass the group node to QgsLayerTreeModel(), or changed the root group of the model with QgsLayerTreeModel.setRootGroup(). In each case Python crashes, with an error message like this:
QObject::connect: Cannot connect (null)::willAddChildren( QgsLayerTreeNode*, int, int ) to QgsLayerTreeModel::nodeWillAddChildren( QgsLayerTreeNode*, int, int )
Is it possible to expose only some group oder layer nodes in the treeview, and how can I implement this?
I'm trying to get my head around the theoretical side of XML and GML. The OGC web page for GML states that:
The Geography Markup Language (GML) is an XML grammar for expressing geographical features.
But what does it mean by a "grammar"? I can't see that phrase (as a noun) used anywhere else. Is it a unique GML meaning?
I see the terms "language", "schema", "format" and "standard". Are they the same thing?
So can I say that GML is an XML schema? Or would that mean something different? If so, is there a single XSD file that defines the entire GML specification?
Answer
When someone designs a class of XML documents for representing information in a particular domain, they will sometimes call this an XML grammar, or a vocabulary, or a schema, or a document type, or even a language. The terminology isn't consistent. There's perhaps a different emphasis: calling it a schema implies that an XML Schema is the primary way in which the grammar/vocabulary is specified; but they all mean essentially the same thing.
I have some points and some features. These points fall outside the features. I want to connect each point with a line to the closest feature. How can this be done using ArcGIS 10.2?
I have checked around toolbox but didn't manage to find the right tool for it.
I am adding Georgia data to my GIS database and, unlike the California data I've added, have decided to use my existing EPSG:3361 (South Carolina StatePlane) coordinate reference system for calculations.
I've done this with reasonable success in North Carolina (South Carolina's northerly neighbor) but never did the research to determine the error suffered by using a linear reference system from a non-intended, albeit adjacent area.
I'm guessing the error shouldn't be more than 50 feet, am I correct in assuming this?
More technical information:
Answer
The northern two-thirds of Georgia should be pretty good, because the South Carolina coordinate reference system is Lambert conformal conic-based. Thus, the standard parallels extend through Georgia too. I ran a point at 31N 85W through the National Geodetic Survey's SPC program to see what the distortion would be. Note: South Carolina's zone is 3900. It returned:
=======================================================
Latitude Longitude Datum Zone
INPUT = N310000.0 W0850000.0 NAD83 3900
=======================================================
NORTH(Y) EAST(X) AREA CONVERGENCE SCALE
METERS METERS DD MM SS.ss
----------- ----------- ---- ----------- ----------
-85051.112 227348.537 SC -2 13 3.35 1.00086349
So the distortion is 1.00086349 or 86 parts per 100000 (or 8.6 parts per 10000). A State Plane zone was originally designed to be no worse than 1 part per 10000.
I am trying to repeat the behavior of this code from the ESRI website.
I want to replicate the output example in ALL_Type field, so that it concatenates a field but it deals (ignores) any that have NULL values.
Unfortunately this example no longer works at 10.2.2 because all example of NULL return a "None" text. Instead I just want a blank to show.
And here is the code.
"*args" allows this routine to accept any number of field values.
# the values are passed as a Python tuple, essentially a
# non-editable list
#
def concat(*args):
# Initialize the return value to an empty string,
# then set the separator character
#
retval = ""
sep = "_"
# For each value passed in...
#
for t in args:
# Convert to a string (this catches any numbers),
# then remove leading and trailing blanks
#
s = str(t).strip()
# Add the field value to the return value, using the separator
# defined above
#
if s <> '':
retval += sep + s
# Strip of any leading separators before returning the value
#
return retval.lstrip(sep)
Answer
I think this should do what you want:
def concat(*args):
sep = "_"
nonnull_args = [str(arg).strip() for arg in args if arg] # Filter NULLs
good_args = [arg for arg in nonnull_args if arg] # Filter blanks
retval = sep.join(good_args)
return retval
The args just get passed through a couple filters and joined by sep at the end.
I am trying to use a cgi script calling mapserver to display a raster layer that is stored in a PostGIS database. However, only a blank white image is returned and I don't know where the problem is.
The request seems to go through, but here is an error in the postgres log that I don't understand and the result is a white image.
My questions are:
Here is some more detailed info for this problem.
The mapserver request is:
http://someurl.com/mapserv.cgi?map=../../usr/maps/test.map&mode=map&layers=kahoo&
The error from the postgres server log:
ERROR: column "o_column" does not exist at character 39
STATEMENT: select o_table_name, overview_factor, o_column, o_table_schema from raster_overviews where r_table_schema = 'public' and r_table_name = 'kahoo' and r_column = 'rast'
From the mapserver log:
.632473 msDrawMap(): rendering using outputformat named png (AGG/PNG).
.632541 msDrawMap(): WMS/WFS set-up and query, 0.000s
.633883 msDrawRasterLayerLow(kahoo): entering.
.859561 msDrawMap(): Layer 0 (kahoo), 0.227s
.859588 msDrawMap(): Drawing Label Cache, 0.000s
.859594 msDrawMap() total time: 0.227s
.862701 msSaveImage(stdout) total time: 0.003s
.862751 mapserv request processing time (loadmap not incl.): 0.231s
.862759 msFreeMap(): freeing map at 0x266bc60.
The mapfile used is:
MAP
CONFIG "MS_ERRORFILE" "../../tmp/ms_error.txt"
DEBUG 5
NAME "test"
SIZE 218 152
EXTENT 739335 2268075 7761335 2284075
UNITS METERS
STATUS ON
PROJECTION
"init=epsg:26704"
END
WEB
METADATA
"wms_title" "WMS Demo Server"
"wms_onlineresource" "http://someurl.com/mapserv.cgi?map=../../usr/maps/test.map&"
"wms_enable_reqest" "*"
END
END
LAYER
NAME kahoo
TYPE RASTER
STATUS ON
DATA "PG:host=*** port=*** dbname='***' user='***' password='***' schema='public' table='kahoo' mode='2'"
PROCESSING "NODATA=0"
PROJECTION
"init=epsg:26704"
END
END
END
Here are the versions I am using: mapserver-6.0.3, postgresql-9.1.6, postgis-2.0.1
I am attempting to integrate a iterative model into a master model to merge all of the outputs of the iterative model into one feature class. However, when I attempt to set the output of the sub-model (iterative) as the input for the merge tool, it will not allow it. I have attached two pictures of each model. Any suggestions?
I am using ArcMap 10.
Thanks,
Eric 
Answer
I have successfully manipulated the sub-model and the outer model to achieve my desired outcome. I ended up getting rid of the Collect Values tool and using the Append tool. I am attaching two images of the models for anyone in the future that might have the same issues. The top image is the iterated sub-model and below it is the composite outer model.
Eric
Is it possible to load external wms services into ArcMap and publish them in ArcGIS Server?
I have tried this and get the error that the layer type is not supported.
This is a follow on from a Question that I asked last year entitled Is there a public OGC Catalog Service (CSW) available?. That Question has an Accept-ed Answer but I would now like to expand upon the question slightly.
Is there a public OGC Catalog Service (CSW) available which serves URLs that can be bound to as Live Data and Maps by ArcGIS for Desktop 10.2?
I would like to learn of a CSW service URL (with its Profile specified) that I could use to Configure the CSW Client of ArcGIS for Desktop 10.2 so that I can see it add a WMS service to ArcMap by using the Add to Map button rather than by using copy/paste of URLs.
A Comment on my Answer to Is there an exhaustive, searchable catalog of all GIS web services (e.g. OGC WMS/WFS, REST, SOAP)? says that it is possible but I do not recall ever being able to do it personally, including during tests I performed yesterday using the latest versions of ArcGIS for Desktop and its CSW Client.
After upgrading to ArcGIS 10.3 for Desktop and @Luke resolving Where to download CSW Client for ArcMap (ArcGIS for Desktop)? at that new version, I was able to conclude my test successfully.
I cannot comment on whether it now works at 10.2.x because I forgot to re-check before upgrading from 10.2.2 and my initial testing was at 10.2.
Answer
URL: http://www.ga.gov.au/geonetwork/srv/en/csw?request=GetCapabilities&service=CSW
Profile: "GeoNetwork CSW 2.0.2 APISO"
Direct link to record highlighted in the screenshots below.
CSW Client config:
Search results:
Added to TOC:
Background
From a known Point, I require to establish the nearest surrounding "visible perimeter" against a table of MultiLineStrings, as shown on the diagram.
I've searched this site with a number of terms (e.g. minimum edge, minimum perimeter, nearest neighbour, clip, containing polygon, visibility, snap, cut nodes, ray-trace, flood fill, inner boundary, routing, concave hull) but cannot find any previous question which seems to match this scenario.
Diagram
Parameters
Methodology Considered
Summary
Whilst this technique will find the nearest intersections, it won't necessarily find all the closest perimeter node points, dependent on the resolution of the radial sweep. Can anyone recommend an alternative method to establish all the perimeter points or supplement the radial sweep technique with some form of buffering, sectoring or offsetting?
Software
My preference is to use SpatiaLite and/or Shapely for the solution but would welcome any suggestions which could be implemented using open source software.
Edit: Working Solution (based on answer by @gene)
from shapely.geometry import Point, LineString, mapping, shape
from shapely.ops import cascaded_union
from shapely import affinity
import fiona
sweep_res = 10 # sweep resolution (degrees)
focal_pt = Point(0, 0) # radial sweep centre point
sweep_radius = 100.0 # sweep radius
# create the radial sweep lines
line = LineString([(focal_pt.x,focal_pt.y), \
(focal_pt.x, focal_pt.y + sweep_radius)])
sweep_lines = [affinity.rotate(line, i, (focal_pt.x, focal_pt.y)) \
for i in range(0, 360, sweep_res)]
radial_sweep = cascaded_union(sweep_lines)
# load the input lines and combine them into one geometry
input_lines = fiona.open("input_lines.shp")
input_shapes = [shape(f['geometry']) for f in input_lines]
all_input_lines = cascaded_union(input_shapes)
perimeter = []
# traverse each radial sweep line and check for intersection with input lines
for radial_line in radial_sweep:
inter = radial_line.intersection(all_input_lines)
if inter.type == "MultiPoint":
# radial line intersects at multiple points
inter_dict = {}
for inter_pt in inter:
inter_dict[focal_pt.distance(inter_pt)] = inter_pt
# save the nearest intersected point to the sweep centre point
perimeter.append(inter_dict[min(inter_dict.keys())])
if inter.type == "Point":
# radial line intersects at one point only
perimeter.append(inter)
if inter.type == "GeometryCollection":
# radial line doesn't intersect, so skip
pass
# combine the nearest perimeter points into one geometry
solution = cascaded_union(perimeter)
# save the perimeter geometry
schema = {'geometry': 'MultiPoint', 'properties': {'test': 'int'}}
with fiona.open('perimeter.shp', 'w', 'ESRI Shapefile', schema) as e:
e.write({'geometry':mapping(solution), 'properties':{'test':1}})
Answer
I have reproduced your example with shapefiles.
You can use Shapely and Fiona to solve your problem.
Point):
from shapely.geometry import Point, LineString
line = LineString([(point.x,point.y),(final_pt.x,final_pt.y)])
from shapely import affinity
# Rotate i degrees CCW from origin at point (step 10°)
radii= [affinity.rotate(line, i, (point.x,point.y)) for i in range(0,360,10)]
from shapely.ops import cascaded_union
mergedradii = cascaded_union(radii)
print mergedradii.type
MultiLineString
import fiona
from shapely.geometry import shape
orlines = fiona.open("orlines.shp")
shapes = [shape(f['geometry']) for f in orlines]
mergedlines = cascaded_union(shapes)
print mergedlines.type
MultiLineString
points = mergedlines.intersection(mergedradii)
print points.type
MultiPoint
from shapely.geometry import mapping
schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
with fiona.open('intersect.shp','w','ESRI Shapefile', schema) as e:
e.write({'geometry':mapping(points), 'properties':{'test':1}})
Result:
points_ok = []
for line in mergeradii:
if line.intersects(mergedlines):
if line.intersection(mergedlines).type == "MultiPoint":
# multiple points: select the point with the minimum distance
a = {}
for pt in line.intersection(merged):
a[point.distance(pt)] = pt
points_ok.append(a[min(a.keys())])
if line.intersection(mergedlines).type == "Point":
# ok, only one intersection
points_ok.append(line.intersection(mergedlines))
solution = cascaded_union(points_ok)
schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
with fiona.open('intersect3.shp','w','ESRI Shapefile', schema) as e:
e.write({'geometry':mapping(solution), 'properties':{'test':1}})
Final result:
I hope that is what you want.
I'm trying to create a polyline shapefile from a text file which provides Name, coordX, coordY. I've created the following script, which I feel will work. However, I can't seem to get it to read the text file properly. I am getting the following error, RuntimeError: Object: CreateObject cannot create geometry from inputs.
import arcpy
arcpy.env.workspace = "C:\\Data"
outFolder = "C:\\Data"
fc = "Paths.shp"
spatRef = arcpy.SpatialReference(26913)
arcpy.CreateFeatureclass_management(outFolder, fc, "POLYLINE", "", "", "", spatRef)
coordinateList = open("C:\\Data\\Paths.txt")
for line in coordinateList.readlines():
print line
pointList = arcpy.Array()
for x, y in coordinateList:
point = arcpy.Point(x,y)
pointList.add(point)
polyline = arcpy.Polyline(pointList)
cursor = arcpy.da.InsertCursor(fc, "SHAPE@")
cursor.insertRow([polyline])
del cursor
Answer
This is the simplest edit to make your code work:
import arcpy
arcpy.env.workspace = "C:\\Data"
outFolder = "C:\\Data"
fc = "Paths.shp"
spatRef = arcpy.SpatialReference(26913)
arcpy.CreateFeatureclass_management(outFolder, fc, "POLYLINE", "", "", "", spatRef)
coordinateList = open("C:\\Data\\Paths.txt")
pointList = arcpy.Array()
for line in coordinateList.readlines():
print line
SplitLine = line.split(',') # break up the string into elements
# 'FredRanch1_1, 529018.125025, 4108038.05548' becomes
# ['FredRanch1_1', '529018.125025', '4108038.05548'] which can be indexed
x = float(SplitLine[1]) # turn the strings into numbers
y = float(SplitLine[2])
point = arcpy.Point(x,y)
pointList.add(point)
polyline = arcpy.Polyline(pointList)
cursor = arcpy.da.InsertCursor(fc, "SHAPE@")
cursor.insertRow(polyline)
del cursor
coordinateList.close() # don't forget to close your file
Using with statements condenses the code and ensures you don't forget to free the cursor and close the file:
import arcpy
arcpy.env.workspace = "C:\\Data"
outFolder = "C:\\Data"
fc = "Paths.shp"
spatRef = arcpy.SpatialReference(26913)
arcpy.CreateFeatureclass_management(outFolder, fc, "POLYLINE", "", "", "", spatRef)
with open("C:\\Data\\Paths.txt",'r') as coordinateList:
with arcpy.da.InsertCursor(fc,'SHAPE@') as cursor:
pointList = arcpy.Array()
for line in coordinateList:
SplitLine = line.split(',') # break up the string into elements
# 'FredRanch1_1, 529018.125025, 4108038.05548' becomes
# [FredRanch1_1, 529018.125025, 4108038.05548] which can be indexed
x = float(SplitLine[1])
y = float(SplitLine[2])
point = arcpy.Point(x,y)
pointList.add(point)
cursor.insertRow(polyline)
With the added complication of starting a new line when the name changes:
import arcpy
arcpy.env.workspace = "C:\\Data"
outFolder = "C:\\Data"
fc = "Paths.shp"
spatRef = arcpy.SpatialReference(26913)
arcpy.CreateFeatureclass_management(outFolder, fc, "POLYLINE", "", "", "", spatRef)
arcpy.AddField_management(fc,'Name','TEXT',field_length=250)
oldName = None
pointList = arcpy.Array()
with open("C:\\Data\\Paths.txt",'r') as coordinateList:
with arcpy.da.InsertCursor(fc,['SHAPE@','Name']) as cursor:
for line in coordinateList:
SplitLine = line.split(',') # break up the string into elements
# 'FredRanch1_1, 529018.125025, 4108038.05548' becomes
# ['FredRanch1_1', '529018.125025', '4108038.05548'] which can be indexed
if oldName == None:
oldName = SplitLine[0] # only should happen on the first iteration
x = float(SplitLine[1]) # turn the strings into numbers
y = float(SplitLine[2])
if oldName.upper() != SplitLine[0].upper():
# the name has changed!
polyline = arcpy.Polyline(pointList)
cursor.insertRow([polyline,oldName]) # insert this line on name change
pointList = arcpy.Array() # reset the array
oldName = SplitLine[0] # make this name the old name
point = arcpy.Point(x,y)
pointList.add(point)
polyline = arcpy.Polyline(pointList)
cursor.insertRow([polyline,oldName]) # clean up the last line
Caveat I have not tested this code, it's from a copy of some of my existing code to do a similar purpose... I may have deleted a necessary line or neglected to change a variable name.
I am trying to update Z values of polygon shape with ArcPy. This script is working without error but increasing the first (0 position vertex) vertex twice. I want to update only once. How can I avoid it?
import arcpy
feature = r"D:\acrpy\VV_test\scripts\map1\ZM_KAD_00.shp" #polygon shape
fields = ["ESRI_OID", "SHAPE@Z"]
whereclause = """{0} = 7""".format(arcpy.AddFieldDelimiters(feature, fields[0]))
z_increase = 2
with arcpy.da.UpdateCursor(feature, fields, whereclause, explode_to_points=True) as cursor:
for row in cursor:
print row[1]
cursor.updateRow([row[0],row[1] + z_increase])
Answer
I wanted to insert a new answer rather than editing my previous one as this comes with a new solution that should handle the problematic of updating the geometries of polygons with more than one inner ring as well.
The OP in the comments pointed out that my other answer didn't work for the case "two donuts in one polygon".
Trying myself I was susprised but today I realized this problem has already come up previously in this other question where a solution is also provided.
I'll highlight the key part that makes the code you are about to see working with all types of polygons (single-to-multipart, zero-to-many inner rings):
arcpy.Polygon() is made up of one or more arcpy.Array()sarcpy.Array() is constructed with at least three arcpy.Point()s (although they require four per definition, as pointed out by @Vince in the comments; yes, in ArcGIS a "Point" means a "Vertex", quite confusing...)arcpy.Polygon() with more than one arcpy.Array() inside is a multipart polygonarcpy.Polygon() with one or more interior rings (alias "donuts"), is NOT a multipart polygon but a singlepart polygon with interior ring(s) (quoting this)a single part arcpy.Polygon() with one or more interior rings is composed by a single arcpy.Array() with a None value where the inner ring(s) begin(s), for example:
[[x1, y1], [x2, y2], [x3, y3], [x4, y4], None, # start of the inner ring coordinatres [x5, y5], [x6, y6], [x7, y7], [x8, y8]]
Last but not least, arcpy.cursors list per each ring (not part, here I mean ring), a starting vertex and an end vertex which are equivalent in terms of coorindates.
As per the previous assessment, one needs to sikp updating (shifting) either the first or last point otherwise the coincident point is shifted twice.
The modified working code of my previous original answer is shown below:
fields = ["OID@", "SHAPE@"]
feature = r"D:\acrpy\VV_test\scripts\map1\ZM_KAD_00.shp" #polygon shape
fields = ["ESRI_OID", "SHAPE@"]
whereclause = """{0} = 7""".format(arcpy.AddFieldDelimiters(feature, fields[0]))
z_increase = 2
desc = arcpy.Describe(feature)
sr = desc.spatialReference
with arcpy.da.UpdateCursor(feature, fields, whereclause) as cursor:
for row in cursor:
print("Feature {}:".format(row[0]))
partnum = 0
# create an array for the updated feature
arr_feat = arcpy.Array()
for part in row[1]:
# create an array for the part
arr_part = arcpy.Array()
print("Part {}:".format(partnum))
i = 0
for pnt in part:
if i == 0:
print("SKIPPING FIRST VERTEX")
arr_part.add(arcpy.Point(pnt.X, pnt.Y, pnt.Z))
i += 1
continue
if pnt:
print("OLD: {}, {}, {}".format(pnt.X, pnt.Y, pnt.Z))
new_z = pnt.Z
new_z += z_increase
i += 1
print("NEW: {}, {}, {}".format(pnt.X, pnt.Y, new_z))
arr_part.add(arcpy.Point(pnt.X, pnt.Y, new_z))
else:
print("Interior Ring:")
null_point = None
arr_part.add(null_point) # MUST DO THIS WHEN INSERTING INTERIOR RING!!!
i = 0
# add part to feature array
arr_feat.add(arr_part)
partnum += 1
polygon = arcpy.Polygon(arr_feat, sr, True)
row[1] = polygon
cursor.updateRow(row)
Note: This question is specifically about installed, desktop software. There is another question specifically about free cloud-based software and services.
What free programs should every GIS user have installed?
I'm not necessarily referring to ESRI extensions or open-source products, but others that increase your productivity and ability to handle GIS tasks.
For example:
Answer
I have a shapefile which contains arround 800 circular polygons with diameter of 60 meters each. Raster is 20 meter resolution so I would have multiple cell values extracted for each polygon. I am using an extract function on multiple raster layers in a for loop, and am putting the results in the dataframe. I know i can get the cellnumbers within the extract function, but I would also like to obtain the feature values (in my case a column in attribute table containing specific polygon code) directly from the shapefile for each pixel value extracted. Thus I would have in one row a pixel cell number, its extracted value and the polygon code.
my code:
dirs=list.dirs(full.names=TRUE)
dirs=dirs[grepl("R20",dirs)]
dirs=as.list(dirs)
DataToWrite=list()
sumirana=data.frame()
for (j in 1:length(dirs)) {
setwd(paste("E:/",as.character(dirs[j]),sep="/"))
files=list.files(full.names=FALSE,pattern="*20m*.tif$")
for (i in 1:length(files)) {
curRaster=raster(paste(getwd(),files[i],sep="/"))
rasterOut=(na.omit(extract(curRaster,azo_plohe, cellnumbers=TRUE,
weights=TRUE, df=TRUE)))
if((length(rasterOut) > 1)) {
DataToWrite=c(DataToWrite,rasterOut[3])
} else {
lista=list(rep("x",15))
tablica_folder=as.data.frame(t(unlist(lista)))
colnames(tablica_folder)=c(1:15)
}
}
if((length(rasterOut) > 1)) {
mylist=as.data.frame(DataToWrite)
naziv_fajla=as.vector(rep((substr((names(DataToWrite)
[1]),50,55)),nrow(mylist)))
tablica_folder=as.data.frame(c(rasterOut[2],naziv_fajla[1],
rasterOut[4],mylist),col.names=c(1:15))
sumirana=rbind(sumirana,tablica_folder)
DataToWrite=list()
} else {
sumirana=rbind(sumirana,tablica_folder)
}
}
I have monthly (and sometimes daily) data on temperatures from about 30 Swedish weather stations in the mid-19th century that I want interpolate to the whole country. I'm using ArcGIS and Stata.
For later years, I have data on additional weather stations (around 100), so I've thought of using the later data as to calibrate a model using the old stations to predict weather at the new stations. Since I know the actual weather at the new stations, I could calibrate to model to achieve the best possible fit. But I'm not sure what a good method for going about achieving a suitable fit (don't want to risk overfitting for example).
Answer
You can likely get a reasonable interpolation using a linear regression (assuming your 30 weather stations are a representative sample) using elevation, latitude and distance from the coast as independent variables with the day as a factor. I've done this using ArcGIS and R previously.
Daily 9am and 3pm temperatures over 10 days in 2003 from weather stations in South Eastern Australia
temp = α + βelev*elev + βlat*lat + βdist*dist) to estimate temperature from elevation, lat and distance. You may need to script this as you'll get different coefficients for each day.From this post, Invalid QgsRasterLayer when attempting to load a Google Static Maps API image into QGIS this block of code can be used to load a Google Static Maps API image into QGIS (using the QGIS Python Console):
>>> reg = QgsMapLayerRegistry.instance()
>>> raster = QgsRasterLayer('http://maps.googleapis.com/maps/api/staticmap?center=51,0&zoom=10&size=2550x3300', 'test')
>>> raster.isValid()
>>> reg.addMapLayer(raster)
This works on Windows (tested versions QGIS 2.14.21 and QGIS 2.18.15). However, this does not appear to work on Mac (tested versions 2.18.9 and 2.14.20) or Linux (version 2.18.13). On Mac or Linux, and invalid raster object is returned.
This is not specific to the Google Static Maps API, an invalid raster object is returned with other http image urls.
Any ideas on how to get this working on Mac/Linux?
When synchronizing an offline Spatialite file with it's corresponding PostGIS database, using the Offline Editing plugin in QGIS, I received an error I hadn't experienced before
The error message is quite clear, I had entered more than 254 characters in a attribute field that in the PostGIS table had this field length.
I looked at the properties of the offline Spatialite table and all the attribute length and precisions are all '0', or undefined. 
In the original PostGIS database table they are defined
Even if the Spatialite attribute that is too long is corrected, the synchronising of the offline Spatialite file now gives a commit error, which means all the changes made offline will now not synchronise.
My question
Is it possible to bring across the attribute field lengths from the PostGIS table to the Spatialite file so as to not allow attributes that are too long to be entered in the first instance?
Answer
With all credit to user30184 who answered in the comments...
SQLite just doesn't care. From sqlite.org/datatype3.html
Note that numeric arguments in parentheses that following the type name (ex: "VARCHAR(255)") are ignored by SQLite - SQLite does not impose any length restrictions (other than the large global SQLITE_MAX_LENGTH limit) on the length of strings, BLOBs or numeric values.
So the answer to the question is No, it is not possible to bring across the attribute field lengths from the PostGIS table to the Spatialite file.
This means that there is no limit to the length of what you can enter in the attribute fields. This is fine if you are working in SQLite and you want to stay working in SQLite, however if you then want to synchronise an SQLite database that was 'offlined' from a PostGIS database you will run into problems if the field lengths are greater than those set up in the PostGIS database.
I want to create a standalone python script using QGIS python API (PYQGIS) for education. First I create a model and export that model to python script.
original python script from the export :
##=name
##raster=raster
##stat=output html
outputs_QGISRASTERLAYERSTATISTICS_1=processing.runalg('qgis:rasterlayerstatistics', raster,stat)
Next I read this Using PyQGIS in standalone scripts and I tried to follow this tutorial to create my script.
Final script :
import sys
from qgis.core import QgsApplication
from PyQt4.QtGui import QApplication
app = QApplication([])
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()
# Prepare processing framework
sys.path.append('/home/username/.qgis2/python')
from processing.core.Processing import Processing
dem='aster.tif'
slope='slope.tif'
aspect='aspect.tif'
outputs_SAGASLOPEASPECTCURVATURE_1=processing.runalg('saga:slopeaspectcurvature', dem,6,1,1,slope,aspect,None,None,None,None,None,None,None,None,None,None)
Processing.initialize()
# Exit applications
QgsApplication.exitQgis()
QApplication.exit()
If I run this script it works fine without error but I don't take some export like aspect and slope (these can't be found?). Any idea why I can't take exports ?
I use Ubuntu and python 2.7 and QGIS 2.18.1.
i try to this script :
import sys
from qgis.core import QgsApplication
sys.path.append('/usr/share/qgis/python/plugins')
import processing
a=processing.alglist()
print a
processing.alghelp("saga:slopeaspectcurvature")
and i take that print any idea why ?
Algorithm not found
but if i use this module in QGIS work fine.
I have a temperature image (interpolated from climatic point data) and NDVI data at the monthly interval.
I need to do correlation analysis between temp and NDVI.
Is it impossible to do this using ArcGIS Desktop, and if so, can you tell me how to do it?
I'm trying to use QgsVectorLayer() function to add delimited text layer in QGIS. So I have a txt file (T_1999_1_2.txt) in path: "D:\PATRICIA\DOCUMENTOS\ESTACOES METEOROLOGICAS\DADOS METEOROLOGICOS\T" and I want to create a shapefile named: T_1999_1_2.shp. The txt file structure is:
par;num;ano;mes;dia;D;btq 3;705;1999;1;2;-999,9;S;40,43885833;-8,43994167
3;619;1999;1;2;25,0;A;41,70972;-8,02699 3;718;1999;1;2;104,0;A;39,78055278;-8,82096667 3;766;1999;1;2;28,0;A;38,76620278;-9,12749444 3;560;1999;1;2;333,0;A;40,71492778;-7,89591667 3;669;1999;1;2;-999,9;S;39,83950000;-7,47866944 3;848;1999;1;2;-999,9;S;38,48486667;-7,47291667 3;555;1999;1;2;8,0;A;41,70655556;-8,80210833
I have tried to run the script but appears the same window error ("Error in file. File cannot be opened or delimiter parameters are not valid") My sript is the following one:
uri = "D:\PATRICIA\DOCUMENTOS\ESTACOES METEOROLOGICAS\DADOS METEOROLOGICOS\T\T_1999_1_2.txt?delimiter=%s&xField=%s&yField=%s" % (";", "x", "y")
vlayer = QgsVectorLayer(uri, "D:\PATRICIA\IG\QGIS\FWI\T\T_1999_1_2.shp", "delimitedtext")
I don't understand, why it doesn't work.
Answer
Change variables indir and outdir according to your need. The code will find each file with extension 'txt' in indir and every subdirectory of indir. If you need another coordinate system than EPSG 4326, please change the EPSG number in line 9. The converted files will be written to directory outdir.
import os
indir = 'G:/LANUV'
outdir = 'G:/LANUV'
for root, dirs, files in os.walk(indir):
for file in files:
if file.endswith('.txt'):
fullname = os.path.join(root, file).replace('\\', '/')
filename = os.path.splitext(os.path.basename(fullname))[0]
uri = 'file:///%s?crs=%s&delimiter=%s&xField=%s&yField=%s&decimal=%s' % (fullname, 'EPSG:4326', ';', 'x', 'y', ',')
layer = QgsVectorLayer(uri, 'my_layer', 'delimitedtext')
QgsVectorFileWriter.writeAsVectorFormat(layer, outdir + '/' + filename + '.shp', 'CP1250', None, 'ESRI Shapefile')
Copy the code and paste it into the Python console.
This code fragment uses code from this thread.
I am working with QGIS Composer 2.0 and I would like to add a second map.
What I would like to do is to have a main map that takes up the majority of the print area. This map, Main Map, will have four layers active from QGIS Desktop.
The second map, Overview Map, is a small map identifying the overview of the area. The overview map will be at a different scale than Main Map, and it will also only have one layer active from QGIS Desktop.
The issue I am having is that when I add a new map, the map areas are updated when I change the layers in QGIS Desktop.
QGIS 2.0 has a new function for Overviews, but the documentation is not overly clear on how to use it.
In ArcMap, you can set up different views, which allows for different data, scale, symbology, etc. ensuring that the different maps in the print layout are independent.
Is there a way to do something like this in QGIS Composer?
Two maps with different layers, from the same QGIS Desktop instance?
Answer
Yes. You can add in a second map to use as an overview map.
You need to lock the features on your your original map first using the Item properties tab. There is a tick box for this in the Main Properties section. You then get QGIS looking how you want the overview map to look and then add it in using the print composer in the same way you added the first map (The "Add New Map" button).
The new map should show up how you want it to look and you can set it as an overview. In the Item Properties for your new map you should find an option for "Overview" (it's about half way down just under the "Show Grid" Option. You need to tell it that the Overview frame you want it to look at is the original map (usually "Map0", but you can set this using trial and error). You can set how you want the overview frame to look.
You can set the scale and extents of this map independently of the first.
You can also add in as many maps as you want.
Using the following two commands yields completely different file sizes:
gdal_translate test.tif -of GTiff -co "COMPRESS=JPEG" -co "JPEG_QUALITY=100" output.tif
and
gdal_translate test.tif -of JPEG -co "QUALITY=100" -co "WORLDFILE=YES" output.jpg.
(I Wrote the commands from my memory, not sure that the syntax is correct)
The jpg file (second command) is 2 times smaller than the compressed tif file. Why?
Answer
By default an RGB image will be written to an RGB color model JPEG image, but this is not actually the most efficient way of writing to JPEG. It is better to convert to the YCbCr color space, and encode that. This is in fact the typical form of standalone JPEGs and what GDAL will produce when writing to a free standing JPEG file. Compressing a 4K x 2.6K image I see:
gdal_translate rgb.tif out.tif -co COMPRESS=JPEG
$ ls -l out.tif
-rw-rw-r-- 1 warmerdam warmerdam 7416934 Nov 4 10:09 out.tif
gdal_translate rgb.tif out_ycbcr.tif -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR
$ ls -l out_ycbcr.tif
-rw-rw-r-- 1 warmerdam warmerdam 2251264 Nov 4 10:09 out_ycbcr.tif
With YCbCr encoding the size of the output is pretty much the same as the JFIF JPEG image.
I am running QGIS 2.18.14 (QGIS code revision: ff83b9a479) on a windows 7 machine, it comes with GRASS 7.2.2 (Code Revision r71496), but it does not come with r.in.lidar. I have noticed QGIS 2.18.16.1 has the same problem and it shares the same GRASS 7.2.2 (Code Revision r71496).
I know there were some issues with it and later solved as I read here, I also certailny used it short ago in 2017, so I wonder which version should I go for before starting the always risky attempt of a new installation.
Answer
When I start editing in ArcMap 10, by default the "Create Features" window is automatically opened and docked. Stopping the edit session hides this window.
When I use the IEditor:StartEditing method this does not display the "Create Features" window however if it is already displayed the IEditor:StopEditing method hides it.
How can I Start Editing programmatically and display the "Create Features" window?
Answer
This should do the trick (C#):
var dockableWindowManager = (IDockableWindowManager)_app; // _app is an IApplication reference
var uid = new UIDClass() { Value = "esriEditor.CreateFeatureDockWin" };
var window = dockableWindowManager.GetDockableWindow(uid);
window.Show(true);
Leaflet and Turf. I am trying find all the layers on a Leaflet map that intersect a point when the user clicks on the map. I am using the turf library to test for this. In the console when I click on the map it returns this error:
turf.min.js:1 Uncaught Error: geojson must be a valid Feature or Geometry Object
at J (turf.min.js:1)
at Object.Lo [as intersect] (turf.min.js:1)
at js2.js:22
at eachLayer (leaflet.js:5)
at e. (js2.js:20)
at e.fire (leaflet.js:5)
at e._fireDOMEvent (leaflet.js:5)
at e._handleDOMEvent (leaflet.js:5)
at HTMLDivElement.r (leaflet.js:5)
Here is the code
map.on('click', function(e) {
var str = "Latitude: " + e.latlng.lat.toFixed(5) + " Longitude: " + e.latlng.lng.toFixed(5) + " Zoom Level: " + map.getZoom();
var lat = e.latlng.lat;
var lng = e.latlng.lng;
var pt = turf.point([lng, lat]);
map.eachLayer(function(layer) {
if (layer instanceof L.GeoJSON) {
var click_intersection = turf.intersect(pt, layer.toGeoJSON());
if (click_intersection) {
console.log(layer.feature.properties)
}
}
});
Reference questions turf.js intersect problem?
Get popup info in multiple layers from one click
update
var c = new L.GeoJSON.AJAX("http://127.0.0.1:8000/childcare_buff_data/",{
style: color(c, "orange", 0.8)})
;
c.addTo(map);
map.on('click',function(e){
lat = e.latlng.lat;
lon = e.latlng.lng;
ProcessClick(lat,lon)
});
var theMarker;
var a;
function ProcessClick(lat,lon){
theMarker = L.marker([lat,lon]).addTo(map);
c.eachLayer(function(layer) {
intersects=turf.intersect(theMarker.toGeoJSON(),layer.toGeoJSON());
if (intersects){
a=layer.feature.properties.buff
console.log(a);
}
})};
I get this error
turf.min.js:1 Uncaught TypeError: Cannot read property 'length' of null
at turf.min.js:1
at turf.min.js:1
at S (turf.min.js:1)
at Pn (turf.min.js:1)
at Object.Lo [as intersect] (turf.min.js:1)
at js2.js:31
at eachLayer (leaflet.js:5)
at ProcessClick (js2.js:30)
at e. (js2.js:22)
at e.fire (leaflet.js:5)
here is the geojson http://www.mediafire.com/file/9fnbz32ib9n1aaj/childcare.geojson/file
Answer
Okay, here is a possible solution, I used, 3 layers, muni, county, state, variables, a,b,c, when I would click on the map and it would grab the polygons feature Name from each layer, Then I posted the three feature Names to a dialog. Here is the process to get the values.I also create a marker to see where I clicked.
map.on('click',function(e){
lat = e.latlng.lat;
lon = e.latlng.lng;
ProcessClick(lat,lon)
});
var theMarker;
var a;
var b;
var c;
function ProcessClick(lat,lon){
if (theMarker != undefined) {
map.removeLayer(theMarker);
};
if (newgeojsonLayer != undefined) {
map.removeLayer(newgeojsonLayer);
};
theMarker = L.marker([lat,lon]).addTo(map); //added to show where you clicked.
muniLayer.eachLayer(function (layer) {
//isInside =turf.inside(theMarker.toGeoJSON(), layer.toGeoJSON());
// Newer method in Turf
isInside =turf.booleanPointInPolygon((theMarker.toGeoJSON(), layer.toGeoJSON());
if (isInside){
a = "City/Town: " + layer.feature.properties.NAME;
console.log("City/Town: " + layer.feature.properties.NAME);
}
});
countyLayer.eachLayer(function (layer) {
isInside =turf.inside(theMarker.toGeoJSON(), layer.toGeoJSON());
if (isInside){
b = "County: " + layer.feature.properties.NAME;
console.log("County: " + layer.feature.properties.NAME);
}
});
stateLayer.eachLayer(function (layer) {
isInside =turf.inside(theMarker.toGeoJSON(), layer.toGeoJSON());
if (isInside){
console.log("State: " + layer.feature.properties.STATE_NAME);
c = "State: " + layer.feature.properties.STATE_NAME;
}
})
}
Basically I used it as a popup for all three layers, I just put the values for a,b,c in the popup dialog.
http://www.gistechsolutions.com/leaflet/DEMO/pointsinpoly/index.html http://www.gistechsolutions.com/leaflet/DEMO/pointsinpoly/indextab.html
For the past few days, I have been working on building a watershed model including drainage basins and find that I am not looking for a global watershed and basin distribution, I am looking to generate those that lead to one specific point.
The idea is to select a point (either by click or creating a point feature or even a polygon feature) and defining the drainage basins as well as watersheds what will inevitable pass by that point.
Using QGIS, I'm wondering if there is a way to create such a tool or if there is a process/algorithm that will give me such a result.
My raw data is a DEM of the area of interest, treated to fill the voids and all.
Answer
Having fiddled some more with the subject, I've come to find an answer to this question and will follow up on it in case someone finds themselves in the same situation.
First off, I have run the GRASS geoalgorithm r.stream.extract, which produces three results: Flow direction (raster), Unique stream ids (point vector) and Unique stream ids (raster). Details on some parameters can be found here.
This previous step was useful for finding the XY points that would work for the next geoalgorithm which is r.water.outlet. This fantastic tool generates the water basin leading to one specific location (exactly what I was looking for, or so I thought) and the results are very impressive and satisfying. I will post an image of this geoalgorithm at the end of this answer.
Basically, the input raster map is the flow direction raster and the Coordinates of outlet point (x, y) are identified either by typing them in or by selecting them with the [...] button beside the text box. Then the GRASS region extent i put in is the extent of my DEM and the cellsize is the cell size of my DEM as well. the output basin looks something like this: (first geoalgorithm, then water basin)
However, in the end I am looking for ALL the water basins that are flowing water into the lake. The progress of this research, for now, is in another post found here
I have this code, most of them getting from this forum. I need to loop over a folder and sub-folders, or just to run a list (putting previously all the raster inside, I don't mind one way or the other).
But I don't find out how to do the loop, I have done it with arcpy in just 30 minutes, but I am fighting with gdal for a couple of days and nothing...
import os, sys
try:
from osgeo import ogr, gdal
from osgeo.gdalconst import *
os.chdir('/home/digd/Desktop/puntos')
except ImportError:
import ogr, gdal
from gdalconst import *
os.chdir('/home/digd/Desktop/puntos')
# open the shapefile and get the layer
driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open('/home/digd/Desktop/puntos/centroids2.shp')
if shp is None:
print 'Could not open puntos.shp'
sys.exit(1)
shpLayer = shp.GetLayer()
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
img = gdal.Open('b5.TIF', GA_ReadOnly)
if img is None:
print 'Could not open landsat.TIF'
sys.exit(1)
# get image size
rows = img.RasterYSize
cols = img.RasterXSize
bands = img.RasterCount
# get georeference info
transform = img.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]
# loop through the features in the shapefile
feat = shpLayer.GetNextFeature()
while feat:
geom = feat.GetGeometryRef()
x = geom.GetX()
y = geom.GetY()
# compute pixel offset
xOffset = int((x - xOrigin) / pixelWidth)
yOffset = int((y - yOrigin) / pixelHeight)
# create a string to print out
s = feat.GetFieldAsString('Name') + ' '
# loop through the bands
for j in range(bands):
band = img.GetRasterBand(j+1) # 1-based index
# read data and add the value to the string
data = band.ReadAsArray(xOffset, yOffset, 1, 1)
value = data[0,0]
#print value
s = s + str(value) + ' '
# print out the data string
print s
# get the next feature
feat.Destroy()
feat = shpLayer.GetNextFeature()
# close the shapefile
shp.Destroy()
Uhmm no, I was trying to do a better explanation of my question.
Anyway, I eventually could do it using numpy and GDAL. I needed to get the pixels values within some polygons, so I convert the polygons to rasters and then I could easily perform different masks and calcs.
I find this process really easy with arcpy, but no idea of how to do this with other tools (well, it's easier and faster and better (at least for me) with the rsgislib, but I have to work in a window enviroment and it's not easy to install this library in windows (at least for me).
Again, I don't know of this is an answer or whatever it's (and truly I don't care), you asked me something and I am just trying to give you the best answer that I can.
I have a recent installation or ArcGIS Desktop on a 64 bit Windows 7 Machine, with 8GB ram, intel core i7 processor, and a solid state drive.
It takes ~2 min to open either ArcMap or ArcCatalog from the start menu (e.g. without loading any data).
As a new installation, I have only the default folder connections.
Is this normal, or is it possible that there is something wrong with my installation?
Answer
After contacting the appropriate IT admin, I recieved the following instructions, resulting in load times closer to my expectation: ArcMap now takes ~4s, ArcCatalog ~3s to load.
run "esriremove.reg" to remove changes made by ArcGIS to the registry, the contents of esriremove.reg are:
Windows Registry Editor Version 5.00
[HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\Session Manager\Environment]
"ARCGIS_LICENSE_FILE"=-
[HKEY_LOCAL_MACHINE\SOFTWARE\ESRI\License]
"LICENSE_SERVER"=-
[HKEY_LOCAL_MACHINE\SOFTWARE\Wow6432Node\ESRI\License]
"LICENSE_SERVER"=-
[HKEY_LOCAL_MACHINE\SOFTWARE\ESRI\License10.0]
"LICENSE_SERVER"=-
[HKEY_LOCAL_MACHINE\SOFTWARE\Wow6432Node\ESRI\License10.0]
"LICENSE_SERVER"=-
re-install the single-user / off-site license (I had started with the single-user / off-site license, so this is not in and of itself the only issue).
I am making a map containing information about public transportation: busses, trams, etc. The map will have a layer containing for example the tram lines. Each line has its own colour, and is represented by a line string geometry.
The problem is, many of the lines have overlapping parts, where more than one tram line cover the same section. To show this to the user, I'd rather want the lines to run parallel to each other instead of being drawn on top of each other. For an example, see how Google Maps show the New York subway lines.
I suspect this is a quite common problem in cartography, but don't know what terminology I should search for.
I am using PostGIS/GeoServer/OpenLayers as my stack, but any open source solution would be acceptable.
I use arcgis 10.1 (arcinfo licence) with arcsde to connect to a SQL Server database. The initial dataset has around 100k polygons (used for land registration, having a resolution that is around 1cm and using an MTM projection). Many geoprocessing operations were made on the dataset to build several related layers (copy / paste / split / merge / join, etc.). Note: I never did any polygon simplification/generalization that would change the shape or amount of vertices in the polygon layers.
At the end of these operations I made a topology (1mm tolerance) between my initial dataset (rank #1) and the final dataset (rank #5) using the "Area Boundary must be covered by boundary of" rule to see if I landed back on my feet. I expected errors only on the polygons that I had modified. However, I got many (several thousand) errors on parts that I had not modified. The symptoms are always the same: they always affect polygon arcs. The width of the sliver/gaps is always between 0,01mm and 10mm. Vertices between my layers are always identical (looking at 1:1 zoom with X,Y coordinates) but the path between vertices is changed for some unknown reason. I have enclosed a few pictures to show the problem (the zoom pictures are at 1:1 zoom). Any idea about what could have happened and how to correct this? The main constraint is that the initial layer must remain unchanged (for legal reasons).
I have tried to change the tolerance of the topology to 1cm - it does solve 50% of my problem, but at the expense of the initial layer being slightly modified (I made an intersection between the original layer and my "initial" layer and they revealed many small changes).
I have also tried the automated "snap" tool using "Edges" and "Vertex" options, but almost nothing was changed (probably because the vertices are already in the right place. Only the arcs between vertices have problems).
Before QgsMapRenderer was deprecated if we wanted to activate OTF we used code like this:
canvas = QgsMapCanvas()
canvas.mapRenderer().setProjectionsEnabled(True)
Currently, QgsMapRenderer is divided into two classes QgsMapSettings and QgsMapRendererJob. I did some research, but with no success. So how we can achieve the same now?
Answer
Not sure if this will be possible as according to one of the QGIS devs @ndawson in response to Feature request #11644:
"Disabling on the fly projection is no longer an option in QGIS 3.0"
I recently wrote a python script that works with raster data and makes some time-consuming operations. My script uses QThread to update a QList and a QProgressBar and everithing works fine if I run the script from PyCharm, the progressbar and the list are updated. Moving the code in a Qgis plugin the gui appears correctly and everything seems to work fine, but if I close and open the plugin again, when I click apply Qgis crashes... To better understand I wrote a code that reproduce the problem. The gui is a QDialog with a QListWidget, a QLineEdit and a QProgressBar. After writing a word in the QLineEdit and clicking on the apply buttom the result should be that the new thread makes a loop using the string and adds every single letter as item in the QListWidget. In the meantime the progressbar increases and I put a sleep time of 1 second to make it visible. I made use also of a QgsMessageBar to alert the user if no word are edited in the QLineEdit. I hope it is enough clear... Here the gui code:
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.gui import QgsMessageBar
import sys, time
class my_worker(QObject):
def __init__(self, word):
QObject.__init__(self)
self.word = word
def run(self):
bar_step = 100/len(self.word)
for i in self.word:
self.emit(SIGNAL("add_item(QString)"), i)
self.emit(SIGNAL("increase_bar(int)"), bar_step)
time.sleep(1)
self.emit(SIGNAL("finished()"))
class my_ui(QDialog):
def __init__(self, parent=None):
"""Constructor."""
super(my_ui, self).__init__(parent)
self.list = QListWidget()
self.messagebar = QgsMessageBar()
input_label = QLabel("Input a word")
self.inputedit = QLineEdit()
self.progress_bar = QProgressBar()
self.progress_bar.setRange(0, 100)
self.progress_bar.setValue(0)
self.button_box = QDialogButtonBox(QDialogButtonBox.Apply | QDialogButtonBox.Close)
layout = QGridLayout()
layout.addWidget(self.messagebar, 0, 0, 1, 4)
layout.addWidget(self.list, 1, 0, 1, 4)
layout.addWidget(input_label, 2, 0)
layout.addWidget(self.inputedit, 2, 1, 1, 3)
layout.addWidget(self.progress_bar, 3, 0, 1, 4)
layout.addWidget(self.button_box, 4, 2, 1, 2)
self.setLayout(layout)
self.connect(self.button_box, SIGNAL("rejected()"), self, SLOT("reject()"))
self.setWindowTitle("Example thread")
def start_worker(self, word):
self.worker = my_worker(word)
self.thread = QThread()
self.worker.moveToThread(self.thread)
self.connect(self.worker, SIGNAL("finished()"), self.workerFinished)
self.connect(self.worker, SIGNAL("add_item(QString)"), self.add_item)
self.connect(self.worker, SIGNAL("increase_bar(int)"), self.increase_bar)
self.thread.started.connect(self.worker.run)
self.thread.start()
def workerFinished(self):
self.progress_bar.setValue(100)
self.worker.deleteLater()
self.thread.quit()
self.thread.wait()
self.thread.deleteLater()
def add_item(self, item):
self.list.addItem(item)
return
def increase_bar(self, bar_step):
self.progress_bar.setValue(self.progress_bar.value() + bar_step)
And here the apply and run methods in the main file of the plugin:
def apply(self):
self.dlg.list.clear()
self.dlg.progress_bar.setValue(0)
if self.dlg.inputedit.text() == "":
self.dlg.messagebar.pushMessage("Missing parameter", 'Please insert a word!',
level=QgsMessageBar.WARNING, duration=2)
return
word = unicode(self.dlg.inputedit.text())
self.dlg.start_worker(word)
def run(self):
"""Run method that performs all the real work"""
self.dlg.list.clear()
self.dlg.inputedit.clear()
self.dlg.progress_bar.setValue(0)
self.dlg.connect(self.dlg.button_box.button(QDialogButtonBox.Apply), SIGNAL("clicked()"), self.apply)
# show the dialog
self.dlg.show()
# Run the dialog event loop
result = self.dlg.exec_()
# See if OK was pressed
if result:
# Do something useful here - delete the line containing pass and
# substitute with your code.
pass
As i said everthing works fine but I get two strange behaviours:
It seems that closing the Dialog it is not deleted and all the changing done during the first opening are preserved.
How can I create a new group in the layer manager using python code?
After I successfully open up an ESRI GDB file/folder with Fiona, what open-source Python package should I use to do a spatial join?
I was told Shapely is not the appropriate package for this.
I am completely new to GIS.
Is there any way in QGIS to apply service area buffers e.g if I have a point on a map I wish to show the service coverage of 1 mile, 2 mile, 3 mile, breaks using the road network?
I have a custom report in ArcGIS that summarizes my fields using basic descriptive statistics. For example, here is a summary field that calculates the maximum value in the group for a particular field:
What I would like to do is calculate the range (|max-min|). However, this function is not included in the summary statistics. I think that I can use the "Fields>Calculated" menu to add custom fields. Here is a screenshot of that:
Does anyone have an idea of what has to go in as the "Formula" for the field calculation, and how do I fill out the element properties for the field I would like to populate with this calculation:
Answer
As recommended by @TimothyMichael I think you will need to follow the help page entitled Calculating fields in reports which advises:
You can use Report Designer to create fields in your report that do not exist in the data source. This is useful if you want to create a dynamic field that is dependent on other attributes within the data source. To create a dynamic field, you can use a C# expression to work with string, date, or numeric field types.
If the issue is that you do not know how to write such a C# expression, then revise your Question to include that detail, and tag it with c# to hopefully attract an Answer from someone who can.
I am new to QGIS. I have been using MapInfo for a couple of years, but I am considering to switch to QGIS. So far I like it, but there are some things that I haven't found out yet.
I want to make multiple graduated chloropleth maps based on a varying set of data but with the same graduated/rule based style.
For example, I have a polygon vector maps with postal code areas and would like to give the polygons a color based on the drivetime to a certain location A.
I can make this style as a graduated style or as a rule based style and save it.
Next, I need to make the exact same map for a different dataset, say drivetime to location B. I want to reuse the same settings (colors en cut-off points for the categories).
In mapinfo this was very easy: Import an excel file with different datasets in each column-> make 1st map and save the thematic map settings -> Make a new map based on a new column but with the saved thematic map settings.
How do I do this in QGIS without -having to remake the graduated style each time (saving a graduated style remembers which column was used and is gone when you change it) or -redefining all rules in the rule-based style?
Or is there a way to make a batch of maps that does work in QGIS 1.8? The tutorial on batching maps does not work anymore in 1.8 (http://www.qgis.nl/2012/07/28/series-kaarten-genereren/?lang=en )
Hope someone can help me.
THanks in advance,
niels
