How can I get a list of all available tables in my own authenticated CartoDB account with a SQL statement?
This returns null:
SELECT table_name as t FROM information_schema.tables WHERE table_schema = 'public';
Alternatively, could it be done with CartoDB.js?
Answer
You cannot do it from unauthenticated queries. You can do it inside the dashboard though with this,
SELECT * FROM CDB_UserTables()
I am trying to generate random points on a route from one point to another. What I have till now is, given a source and a destination, I can get a route using Google's API. But how do I get the coordinates of random points on the route?
I am creating an application where the user picks a value from a combo box, the application then selects that feature and then selects the corresponding features from a second layer, (both selections are in a different colour). What I want to do is as soon as the selection is made, the map will zoom to the envelope of the selected features. I know how to zoom to 2 individual layers (using pEnvlope.union...), but i need to do it on features.
Has anyone got any idea of the best way to do this? (I'm using ArcGIS 9.3.1 and VBA).
Thanks
Answer
You could just use the ArcMap Command, Zoom To Selected Features.
Dim pUID As New UID
Dim pCmdItem As ICommandItem
' Use the GUID of the Save command
' Zoom to all selected features
pUID.Value = "{AB073B49-DE5E-11D1-AA80-00C04FA37860}"
' or you can use the ProgID
' pUID.Value = "esriArcMapUI.ZoomToSelectedCommand"
pUID.SubType = 3
Set pCmdItem = Application.Document.CommandBars.Find(pUID)
pCmdItem.Execute
I have thousand of tif tiles and a vrt file, generated by using gdalbuildvrt. At the moment they are at the same folder level in the system, but if I open the vrt file with Notepad I cannot find any url, address or location. Does this mean that the vrt file and the tiles cannot be separated from each other and have to be together in the same folder? how can I edit that in the vrt file?
Answer
Depending on where you create the VRT, it will either become a relative path, or an absolute path. You can manually set this, by modifying the relativeToVRT="1"to a 0, and then write a complete path in the instead of just the image filename.
See the example below of a full path VRT.
Red
D:/A/Bunch/of/Folders/Image.tif
1
So I have vector data for quad lines, township/range lines and section lines for this area. The data is already projected into NAD83 UTM zone 10N and the lines are no longer straight. They warp pretty severely when going over hilly areas for instance.
What I've been asked to georeference are some 1860'S General land office plats of the area that have township and section lines on them.
If I bring the GLO plats into my MXD it becomes very hard to georeference the very square ancient maps to the modern squiggly lines (caused by modern projection) of the vector data.
I figured I could flatten the vector data out in another environment and georeference the plats to the flattened vector lines then project the newly georeferenced plats into the desired UTM zone 10N projection.
What projection, geographic coordinate system (or tool?) would be optimal for accomplishing this?
Edit: ArcMap 10.1. Student version.
Edit2: By not straight I mean on lets say a printed topo map the section, township and range lines are straight up and down. Those same lines when projected are not straight up and down, they go around elevation changes and are no longer square. I want to georeference the lines from the GLO plats to the lines on the vector data and was hoping I could flatten it out to do this easier. The content of the GLO plats other than the lines is suspect because when they were made the features were only recorded where they crossed section or township/range lines, the rest was essentially guessed at by whoever created the plats.
I've been scouring the internet for the past few hours, and I just can't seem to figure it out!
I am using QGIS 2.4.0-Chugiak.
Here are my simple steps:
1) I create a new ShapeFile (Type: Polygon) and add 3 attributes: 
2) I draw a simple polygon with 6 coordinate points and give it the following attributes: 
3) I now export the Shapefile as a KML and it gives me the following text code:
test1
test1
9
Miami
relativeToGround relativeToGround -73.977273863779757,40.764000085624744 -73.976720692372155,40.764083880580976 -73.976333472386827,40.763925601130516 -73.976480984762205,40.763781287185402 -73.976837473002647,40.763613696404278 -73.977390644410249,40.763674215346221 -73.977273863779757,40.764000085624744
4) Finally I go ahead and re-import what I just exported, using "Add Vector Layer" hoping to get all the attributes back, but I get none of them!
I just want to get all the attributes imported.
What am I doing wrong?
Answer
This seems to be a long running bug in QGIS: How to convert KML to shapefile without losing attributes using QGIS?
I am pretty sure it worked some time, but QGIS 1.8 and 2.2 show the same behaviour you mention.
Add vector Layer is using OGR/GDAL in the background, and that has two drivers for kml reading (KML and LIBKML), but LIBKML is not included in OSGEO4W builds.
On the other side, ogr2ogr from gisinternals/sdk is able to read the kml file, and convert it into a valid shapefile with all attributes. They have added the LIBKML driver. The same applies for the ubuntugis unstable QGIS 2.4.0 version on Ubuntu 14.04.
There are some bug report for it: http://hub.qgis.org/issues/8273 and http://trac.osgeo.org/osgeo4w/ticket/291
In ArcMap 10.2, I have a layer that has two attributes: depth and magnitude.
Depth should be represented by color: deeper points being lighter and shallower being darker.
Magnitude represented by symbol's size: points with a higher magnitude have a bigger symbol and points with a lower magnitude have a smaller one.
How do i do that?
So, many people has told me to use "multiple attributes" in the symbology dialog. I wish i could upload a screenshot but my firewall seems to have the image service blocked. Let's try with words.
Basically, whenever i select something on "symbol size" or "color ramp", the same selection is duplicated on the other button.
Answer
The answer is here! I am also having the same problem and there is an official answer.
https://support.esri.com/en/technical-article/000010616
In short, in the 'value fields' choose only ONE variable, in the 'variation by-symbol size' choose another variable
Full procedure:
In Table Of Contents, right-click the layer and click Properties.
In the Layer Properties dialog box, click the Symbology tab. On the left side (the Show: box), click Multiple Attributes.
In the Value Fields section, select the field containing the attribute to base the color symbology on. Leave the other two fields blank.
Click the Add All Values button at the bottom and uncheck .
Click the Symbol Size button. The Draw quantities using symbol size to show relative values dialog box is displayed.
Select the Value field and set it to the field that has quantities. Leave the Normalization as none. Click OK.
I have a PostGIS table of polygons where some intersect with one another. This is what I'm trying to do:
select the_geom from the_table where ST_Intersects(the_geom, (select the_geom from the_table where source_id = '123'))Any ideas?
Answer
Since you said you get a group of intersecting polygons for each polygon you're interested in, you may want to create what is referred to as a "polygon overlay".
This isn't exactly what Adam's solution is doing. To see the difference, take a look at this picture of an ABC intersection:
I believe Adam's solution will create an "AB" polygon that covers both the area of "AB!C" and "ABC", as well as an "AC" polygon that covers "AC!B" and "ABC", and a "BC" polygon that is "BC!A" and "ABC". So the "AB", "AC", and "BC" output polygons would all overlap the "ABC" area.
A polygon overlay produces non-overlapping polygons, so AB!C would be one polygon and ABC would be one polygon.
Creating a polygon overlay in PostGIS is actually pretty straightforward.
There are basically three steps.
Step 1 is extract the linework [Note that I'm using the exterior ring of the polygon, it does get a little more complicated if you want to correctly handle holes]:
SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines
Step 2 is to "node" the linework (produce a node at every intersection). Some libraries like JTS have "Noder" classes you can use to do this, but in PostGIS the ST_Union function does it for you:
SELECT ST_Union(the_geom) AS the_geom FROM (...your lines...) AS noded_lines
Step 3 is to create all the possible non-overlapping polygons that can come from all those lines, done by the ST_Polygonize function:
SELECT ST_Polygonize(the_geom) AS the_geom FROM (...your noded lines...)
You could save the output of each of those steps into a temp table, or you can combine them all into a single statement:
CREATE TABLE my_poly_overlay AS
SELECT geom FROM ST_Dump((
SELECT ST_Polygonize(the_geom) AS the_geom FROM (
SELECT ST_Union(the_geom) AS the_geom FROM (
SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines
) AS noded_lines
)
)
I'm using ST_Dump because the output of ST_Polygonize is a geometry collection, and it is (usually) more convenient to have a table where each row is one of the polygons that makes up the polygon overlay.
Sometimes the .gpx file that you obtain from your GPS-Device does not deliver the speed. The file looks like this:
[...]
518.38
[...]
The speed is missing. How is it possible to obtain the speed?
Answer
Download GPSBabel - for Ubuntu just use:
sudo apt-get install gpsbabel
With this command you will get a .gpx file of all track points with the speed (meter/second):
gpsbabel -t -i gpx -f input.gpx -x track,speed -o gpx -F output.gpx
EDIT: please see comment by drnextgis - the command requres the addition of ',gpxver=1.0' and should therefore be as below:
gpsbabel -t -i gpx -f input.gpx -x track,speed -o gpx,gpxver=1.0 -F output.gpx
The .gpx will look like this then:
[...]
624.630000
2.338139
[...]
It is also possible to obtain other output formats like .csv:
gpsbabel -t -i gpx -f input.gpx -x track,speed -o unicsv -F output.csv
The result will look like this then:
No,Latitude,Longitude,Altitude,Speed,Date,Time
1,59.414697,3.693597,521.3,0.00,2016/07/17,10:30:14
2,59.414704,3.693641,521.3,3.28,2016/07/17,10:30:15
3,59.414789,3.694040,524.6,2.34,2016/07/17,10:30:28
[...]
I am trying to accomplish PyQGIS simple map rendering via the QGIS 2.0 Python console, but end up with only a blank image. I've followed instructions of this link. My objective is to export a simple image of the current extent of layers of interest to a png. Where can I be going wrong?
#identify crs
crs = QgsCoordinateReferenceSystem(3297, QgsCoordinateReferenceSystem.EpsgCrsId)
# create image
img = QImage(QSize(800,600), QImage.Format_ARGB32_Premultiplied)
# set image's background color
color = QColor(255,255,255)
img.fill(color.rgb())
# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)
render = QgsMapRenderer()
#set crs
render.setDestinationCrs(crs)
# set layer set
layers = qgis.utils.iface.legendInterface.layers()
lst = []
for layer in layers:
lst.append(layer.())
render.setLayerSet(lst)
# set to CURRENT extent**map of layer's full extent is successfully rendered when QgsRectangle(render.fullExtent() is used)
rect = QgsRectangle(render.extent())
rect.scale(1.1)
render.setExtent(rect)
# set output size
render.setOutputSize(img.size(), img.logicalDpiX())
# do the rendering
render.render(p)
p.end()
# save image
img.save("render.png","png")
I know that ArcGIS is available for Android, but does anyone know of any other GIS apps that are available for Android tablets?
What I have to do for migration of my workspace from one GeoServer to another (from 15.2 on server#1 to 16.2 on server#2)? What I do: -replace directory tomcat/webapps/geoserver/data/my_workspace_dir -restarted tomcat. What I have had? this: -have had workspace in GUI geoserver -have had styles in GUI geoserver (but all are broken, when I opened them, they throw this: org.apache.wicket.WicketRuntimeException: Can't instantiate page using constructor 'public org.geoserver.wms.web.data.StyleEditPage(org.apache.wicket.request.mapper.parameter.PageParameters)' ... - no one source (DB) or layers from my workspace
What I have to do?
Answer
You cannot move a single worskspace around by simple file copies, the unit that can be moved is the data directory.
You may want to try the backup/restore plugin, but mind, it's an unsupported community module, not sure it it will work or not: https://docs.geoserver.org/latest/en/user/community/backuprestore/usagegui.html
I am using QGIS 2.18.13.
I have a polygon layer with building footprints, and a point layer with a given number of points. I would like to draw lines (randomly) between two points and see how many buildings the line crosses. How I can do that with QGIS?
Answer
One approach is through Virtual Layers / SpatiaLite syntax.
An example with 12 building and three points.
For illustration purpose add connecting lines between points with the Virtual Layer SQL:
select p1.id as id1, p2.id as id2, MakeLine(p1.geometry, p2.geometry) as geometryline
from point p1
cross join point p2
where p1.id <> p2.id
Use the Add / Edit Virtual Layer button:
In the QGIS Create a virtual layer dialog it looks like:
It will give you:
Add another virtual layer doing the complete query. Remember to rename the layer so you won't overwrite the first virtual layer:
with
crossinglines(geometry) as (
select distinct MakeLine(p1.geometry, p2.geometry) as geometry
from point p1
cross join point p2
)
select distinct b.geometry
from building b
inner join crossinglines cl on st_crosses(cl.geometry, b.geometry)
In the layer control context menu Show feature count for the last virtual layer - showing the number of buildings with a line crossing.
The lines of the first SQL are doubled. This is handled in the second SQL.
Notice that the virtual layers are dynamic to their base layers. Adding a point will recalculate the virtual layers when they are refreshed with the QGIS refresh button.
UPDATE:
Now for adding the building id to the building being crossed add a new virtual layer with the following SQL.
select v2.geometry, b.osm_id
from virtual_layer2 v2
inner join building b on st_equals(b.geometry, v2.geometry)
Where virtual_layer2 is you joined building layer that crosses the lines. Change the column name osm_id to your building id.
In what coordinate system and unit measures the coordinates returned by the Google Maps API are? I ask this because I have a page that returns the geocode lat / long for me to find in my postgis basis. But the query does not return anything, I need to know which projection should I work with my data.
Answer
Google Maps, like most web maps, works in WGS84, Web Mercator (Auxiliary Sphere) - EPSG 3857. http://spatialreference.org/ref/sr-org/epsg3857-wgs84-web-mercator-auxiliary-sphere/
Planar units are meters.
Also, see this similar question.
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3857"]]
I have a buffer layer (green polygon) which I want to split to two polygons whenever it crosses a barrier (blue line). I have been trying to use "splitGeometry" method, but I just can't get it to work. My code so far is this :
while ldbuffprovider.nextFeature(feat):
while barprovider.nextFeature(feat2):
if feat.geometry().intersects(feat2.geometry()):
intersection = feat.geometry().intersection(feat2.geometry())
result, newGeometries, topoTestPoints=feat.geometry().splitGeometry(intersection.asPolyline(),True)
Which returns 1 for result( error) and an empty list for newGeometries. Any help is greatly appreciated.
Answer
You can use the reshapeGeometry function of the QgsGeometry object for this, which cuts a polygon along its intersection with a line.
The following will intersect the buffer polygons with the lines, and add the split polygon features to a memory layer (QGIS 2.0 syntax):
# Get the dataProvider objects for the layers called 'line' and 'buffer'
linepr = QgsMapLayerRegistry.instance().mapLayersByName('line')[0].dataProvider()
bufferpr = QgsMapLayerRegistry.instance().mapLayersByName('buffer')[0].dataProvider()
# Create a memory layer to store the result
resultl = QgsVectorLayer("Polygon", "result", "memory")
resultpr = resultl.dataProvider()
QgsMapLayerRegistry.instance().addMapLayer(resultl)
for feature in bufferpr.getFeatures():
# Save the original geometry
geometry = QgsGeometry.fromPolygon(feature.geometry().asPolygon())
for line in linepr.getFeatures():
# Intersect the polygon with the line. If they intersect, the feature will contain one half of the split
t = feature.geometry().reshapeGeometry(line.geometry().asPolyline())
if (t==0):
# Create a new feature to hold the other half of the split
diff = QgsFeature()
# Calculate the difference between the original geometry and the first half of the split
diff.setGeometry( geometry.difference(feature.geometry()))
# Add the two halves of the split to the memory layer
resultpr.addFeatures([feature])
resultpr.addFeatures([diff])
I have many raster files each one containing information that is different from the others (e.g. one contains information on rain, the other on altitude, and another on land use). The cells in each raster file may or may not line up with the cells of the other raster files.
My goal is the following: I need to create a text file where each row represents a cell and, along the columns for that row, contains the information in that cell on rain, altitude, land use, etc.
What workflow would you suggests to "merge" these raster files together so that their cells line up and contain all the information available across all of the rasters?
I am using ArcGIS 10.1
I have created a fishnet grid (3km x 3km) over a vector polygon.
I now want to clip the fishnet over the vector polygon retaining all grid squares within the boundary of the polygon.
Therefore I want all squares that would typically be clipped into partial squares to remain as whole squares.
Does anybody know how to do this?
Answer
You are mixing a lot of terminology which is making your question rather confusing. You talk about "clipping" but then state "within the boundary" and finally you say you want all whole squares, that's 3 different things!
I'm assuming that you want all squares that are intersecting the polygon, which would include all whole squares completely within the polygon and whole squares touching the boundary of the polygon. Therefore, what you want is to use the Select By Location tool with the relationship of Intersects.
This will select those squares which you then export to a new dataset.
My team is building a web GIS portal using GeoNode. We're trying to add a "go to XY" function or widget into MapLoom or GeoExplorer, so users can enter coordinates and pan/zoom to that location on map.
Does anyone know particular script, tool, or API to make this work?
I have to get an extent of a given raster so that I can save it and use the same to clip another raster with that extracted extent. As of now I am stuck in the extent part.
import arcpy
from arcpy import env
desc= arcpy.env.extent("D:\GIS @ UTD\Sem 1\GISC 6317\Lab\Lab 10\temp\Tahoe\Tahoe\Emer\erelev.grid")
print desc
The following gives an output 0 0 0 0 NaN NaN NaN NaN.
Any suggestions.
Answer
Dont use env.extent you need to get raster extent.
import arcpy
elevRaster = arcpy.sa.Raster('C:/data/elevation')
myExtent = elevRaster.extent
print myExtent
i hope it helps you...
I am trying to import a CSV file with Lat Long values I geocoded in Google Earth into QGIS, and I'm getting some very strange results.
A well-formatted decimal-degree Lat/Long pair (like 30.402716, -88.867475) gets translated into a tiny tiny XY pair like (-0.0008,0.0003). The points show up basically in the right position relative to each other, but these decimal values are way off relative to my other layers, why is this decimal conversion happening and why is it so wrong?
I'm using the "Create Layer from Delimited Text" dialogue, and I set X to Long, and Y to Lat, and selected the CRS to Google Maps Global Mercator EPSG:900913
Answer
If they're truly lat/lon coordinates, set the CRS to EPSG:4326 and see if that helps. If you want EPSG:900913, you'll need to ("re")project them.
I'm having trouble using the Dimension Values (Optional) argument in arcpy.MakeNetCDFRasterLayer_md.
Whatever I seem to do the output raster is always the first layer which is the default value.
Below is my attempt to export all layers to raster in the time dimension.
Code:
def extractAllNetCDF():
variable = "RRt_10m"
x_dimension = "lon"
y_dimension = "lat"
band_dimension = ""
dimension = "time"
valueSelectionMethod = "BY_VALUE"
outLoc = "E:/New Folder/"
inNetCDF = "E:/netCDFFiles/RRt.nc"
nc_FP = arcpy.NetCDFFileProperties(inNetCDF)
nc_Dim = nc_FP.getDimensions()
for dimension in nc_Dim:
top = nc_FP.getDimensionSize(dimension)
for i in range(0, top):
if dimension == "time":
dimension_values = nc_FP.getDimensionValue(dimension, i)
nowFile = str(dimension_values)
arcpy.MakeNetCDFRasterLayer_md(inNetCDF, variable, x_dimension, y_dimension, nowFile, band_dimension, dimension_values, valueSelectionMethod)
arcpy.CopyRaster_management(nowFile, outLoc + nowFile + ".img", "", "", "", "NONE", "NONE", "")
print dimension_values, i
The print method at the end will show the dates as they should be and the index i is also moving alon so there is no reason to think that there are other problems with the code other than the Dimension Value being incorrect and reverting to the default.
Does anyone have any idea how to get the subsequent layers to export?
Is there any code online that has a specific example of this using the Dimension Value argument other than empty quotes?
Answer
Okay, so this was a little tricky as there seemed to be a few different ways of using what is descibed as a Value Table on the ESRI Help Page:
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//004300000006000000.htm
There are also a few uses of brackets that could be tricky if you are not really sure of the syntax so hopefully if you can just use the example below it will work.
So what you really want is answers so here it is with a fix
def extractAllNetCDF():
variable = "RRt_10m"
x_dimension = "lon"
y_dimension = "lat"
band_dimension = ""
dimension = "time"
valueSelectionMethod = "BY_VALUE"
outLoc = "E:/New Folder/"
inNetCDF = "E:/netCDFFiles/RRt.nc"
nc_FP = arcpy.NetCDFFileProperties(inNetCDF)
nc_Dim = nc_FP.getDimensions()
for dimension in nc_Dim:
top = nc_FP.getDimensionSize(dimension)
for i in range(0, top):
if dimension == "time":
dimension_values = nc_FP.getDimensionValue(dimension, i)
nowFile = str(dimension_values)
#THIS IS THE NEW CODE HERE
dv1 = ["time", dimension_value]
dimension_values = [dv1]
#END NEW CODE
arcpy.MakeNetCDFRasterLayer_md(inNetCDF, variable, x_dimension, y_dimension, nowFile, band_dimension, dimension_values, valueSelectionMethod)
arcpy.CopyRaster_management(nowFile, outLoc + nowFile + ".img", "", "", "", "NONE", "NONE", "")
print dimension_values, i
So that's it essentially. There is no need to create an instance of 'Value Table' type e.g.
vtab = arcpy.ValueTable(2)
as is seemed to be implied by the fact the argument was labelled 'Value Table'.
There is no need to use all of the brackets that they show in the examples whether curly, round or otherwise. Follow the above and it should work.
I'm working with this downloadable shapefile that has X Y coordinates. They are akin to '3672187.92698000, 534175.72095400'.
I would like to convert them to longitude latitude so they are more like '-90.097017, 29.963176'.
I've seen this question tackled using ArcMap however I don't have that software. I was able to download and install QGIS but I am unfortunately perplexed by its complicated interface. Would like to do the conversion with it, if possible.
Answer
The dataset you mention is a shapefile, a format invented by ESRI, but understood by most GIS software, including QGIS.
After extracting the zip, you can add it with Add vector layer and point to the .shp file. The CRS information is stored in the .prj file, and the layer CRS will automatically set right by QGIS. In your case, NAD_1983_StatePlane_Louisiana_South_FIPS_1702_Feet with US feet as units.
With the openlayers plugin, you can add a Openstreetmap or Google background layer. For doing that, you have to set the project CRS to EPSG:3857.
If you want coordinates in lat/lon degrees, just rightclick on the shapefile layer, and Save as ... to a new file under a different name, selecting EPSG:4326 as CRS for that, and check to add that layer to the canvas. Saving may take some time.
For the next step, you better zoom in to see just a couple of points. Open the attribute table, and click on the pencil symbol at the bottom to enter the edit mode, and then the field calculator icon bottom right. Create a new field named degx, type real, precision 6, and select $x from geometry. After saving (which takes some time), do the same for degy and $y. Leave edit mode, then the attribute table.
The new columns in the attribute table give you lat and lon in degrees.
I've been scouring the internet for a fast, open-source cost-distance function that I can embed inside my code. The functionality should basically work like the ArcGIS implementation where I can specify a source raster and cost raster and output a cost distance raster. That said, the implementation does not need to be as sophisticated as ArcGIS in that I don't need to necessarily input "rasters" with embedded geo metadata - a signature that takes simple numeric arrays will suffice since I can verify that the data overlaps correctly at the time I call the function.
My specific intent is to calculate the cost-distance to a single point in the center of a 1000 x 1000 cost raster based on a user click, so fast execution is very important.
A C# implementation would be ideal but I'll look at anything that's out there.
Does anyone know of an open-source library that supports this? Thanks for your help!
Answer
GRASS GIS has a C implementation in r.cost (source, documentation) which uses a min-heap. Alternatively, you could use a graph package like QuickGraph and Floyd-Warshall to compute the cost.
Recent changes in GRASS 6.4 have made r.cost significantly faster, so perhaps performance may be good enough: on my laptop, it takes about 3s for a 1M cell region, or 5s with knight's move enabled. GRASS is a C application, not a drop-in solution for a C# codebase. If you're OK adding to your stack, you could use PyWPS to make calls to GRASS, and then use the result elsewhere in your application.
I have a set of aerial photos of vegetation which needs to be classified. Instead of using manual interpretation, I want to try using Object-Based classification based on the Orfeo Toolbox in QGIS Processing.
However, I cannot find any tutorial that explains how to do it in QGIS. In the Orfeo guide, the steps are
But I cannot find any specific step-by-step, which geo algorithm to use. I tried to follow LSMC (Large-Scale Mean-shift Classification) workflow, and also run Segmentation (Meanshift). The LSMC result was quite good, but I have no idea what to do next to label the object (or make a training object) and classify it. 
I tried the TrainImageClassifier(SVM), TrainORGLayerClassifier, ORGLayerClassifier, but there were no outputs came out. I tried to find some tutorials, but could not find any. I am totally confused about the workflow.
Answer
I found a tutorial here
But it is not that helpful, as when I am in the preparation of reference data (join attributes by location), it results in shapefile, and there is no XML file. Meanwhile, the next step requires an XML file. Still stuck.
In the end, I classified them based on ruleset like eCognition, but have to write the script down in field calculator
The XML file in the next step is for output ... doc. says:
"Output XML file: XML filename where the statistics are saved for future reuse."
SORRY, this manual is of different tool... really seems, that something is missing in that tutorial. I will try to use scikit-learn. As I have now layer with segments and its features (in 4th step of segmentation I have used as an input stack of layers, which I want to use as classification features). I will report here.
Right now I am doing address resolution with a ESRI file called Street_Addresses_US.loc which works. I have been given a shape file that contains streets that was compiled by the state. Is it theoretically possible that I could use this streets shapefile for address resolution? Or do shapefiles never contain the necessary data for address resolution?
By address resolution I mean put in an address get a lat / long. ArcEngine 10 C# VS2010.
ESRI support suggested converting a shapefile into a geodatabase file (.gdb) or (.mdb) so I went looking for that, and I'm going to post it here for the community
public void ShapeFileToAccess()
{
//output
IWorkspaceName pWorkspaceName = new WorkspaceNameClass() as IWorkspaceName;
pWorkspaceName.WorkspaceFactoryProgID = "esriCore.AccessWorkspaceFactory";
pWorkspaceName.PathName = "c:\\Data\\test.mdb"; //this is the output file that is created
IFeatureClassName pFeatureClassName = new FeatureClassNameClass();
IDatasetName pDataSetName = pFeatureClassName as IDatasetName;
pDataSetName.WorkspaceName = pWorkspaceName;
pDataSetName.Name = "test";
//input
IWorkspaceName pInShpWorkspaceName = new WorkspaceNameClass() as IWorkspaceName;
pInShpWorkspaceName.PathName = "c:\\Data\\Shapefiles";
pInShpWorkspaceName.WorkspaceFactoryProgID = "esriCore.ShapefileWorkspaceFactory";
//define the dataset
IFeatureClassName pFCName = new FeatureClassNameClass() as IFeatureClassName;
IDatasetName pShpDatasetName = pFCName as IDatasetName;
pShpDatasetName.Name = "Roads.shp"; //this is your input file
pShpDatasetName.WorkspaceName = pInShpWorkspaceName;
//convert
IFeatureDataConverter pShpToFc;
pShpToFc = new FeatureDataConverterClass();
pShpToFc.ConvertFeatureClass(pFCName, null, null, pFeatureClassName, null, null, "", 1000, 0);
}
Answer
It is possible that the attributes necessary are available in the shapefile to do a reverse address lookup, but it's not required. I've seen shapefiles that didn't have any data other than the shape, so there's no guarantee that enough data will be available for you to perform a lookup with the shapefile alone.
Is it possible to get the main bearing of buildings by floor plan?
I want to have a building (by ID) and it's main bearing (e.g. E/90° or W/270°). The idea is to change the polygon-feature class into a line-feature class. Then I can calculate the lengths of each side. I can also calculate the bearing of each side. Buildings with a more complex floor plan I planned to cumulate the side's lengths on the bearing.
The problem is that one building with four parallel sides get four different bearings because of the starting point for the bearing-calculation (see image). Furthermore the calculation doesn’t always start at the same end of a line.
Answer
ArcGIS version 10 has "Minimum Bounding Geometry" in Data Management, Features toolset. I see an option RECTANGLE_BY_WIDTH that produces a new field MBG_Orientation—The orientation of the longer side of the resulting rectangle. If this doesn't give you the correct orientation, there is also "Update COGO attributes", for each line. You can then use Summary Statistics or arcpy.Statistics_analysis to get the angle of the Max Length, for each polygon ID.
We are using QField 0.7.5 that seems a really great project to bring QGis on a mobile device.
I prepared a QGis 2016 project : orthophoto + 1 polygon layer (casdatre) + 1 point layer to input data.
I transferred the project and the data to a Galaxy Tab A (Android 5.0.2). I open the project with QField and everything displays fine, I can move and query the data, really great.
But when I create a point, I can't save the new point (save don’t save the new point).
Do you have any idea what’s wrong?
Answer
There are potentially two things that could trigger this behavior.
The project might be located on an external SD card. Android prevents apps from writing to external SD cards except for an app specific folder. If you save the project on an external SD card, put your data into the folder /Android/data/ch.opengis.qfield/files See also: Documentation, Issue report
The layer in which you digitize is not in the same CRS as the project itself. This is mentioned in the documentation. It is also reported as a bug in the current version that newly digitized geometries are not reprojected. Update: This issue is fixed in QField >= 0.8
I was just wondering if there was a good way to convert GeoJSONdata to WKT data using OpenLayers.
The method I've tried is to create a Feature.Vector with the GeoJson and then tried using WKT.write(feature) to get the WKT but it's not working. I'm obviously missing something.
var geojson_format = new OpenLayers.Format.GeoJSON();
var testFeature = geojson_format.read(val);
var wkt = new OpenLayers.Format.WKT(wkt_options);
wkt.write(testFeature);
Could anybody point me in the right direction?
My GeoJSON is
{
"type": "FeatureCollection",
"features":
[
{"type":"Feature", "geometry": {"type":"MultiPolygon","coordinates":[[[[-0.12072212703174,51.51899157882951],[-0.119842362474856,51.51708218165009],[-0.117203068804201,51.51314296254295],[-0.117803883623537,51.51300942372214],[-0.118125748705324,51.51283582266979],[-0.118533444475588,51.51254203476644],[-0.118855309557375,51.51192774484595],[-0.118898224901614,51.51126002901546],[-0.117846798967776,51.50983108425668],[-0.119220089983401,51.50931024277883],[-0.120335888933596,51.508549002215815],[-0.121151280474123,51.5077610380613],[-0.122095418047365,51.50653232100095],[-0.126429867815432,51.50729359525858],[-0.127266717028078,51.507440506370465],[-0.127266717028078,51.507801104364084],[-0.127266717028078,51.508094922841394],[-0.127331090044436,51.50848222612105],[-0.127567124437746,51.50939037262451],[-0.128146481584963,51.509764310040985],[-0.128511262010989,51.51021837277772],[-0.128597092699465,51.51120661132621],[-0.128597092699465,51.51191439062526],[-0.129004788469729,51.51260880491084],[-0.129584145616946,51.51374388239237],[-0.130120587419924,51.51494569831],[-0.130614113878664,51.51653471734371],[-0.126794648241457,51.51701541806259],[-0.125507187914309,51.51718900318654],[-0.121001076769289,51.519178508517115],[-0.12072212703174,51.51899157882951]]]]}}
]
}
Answer
Not sure what you mean by "it's not working", I'm getting output that looks correct:
var val = ''; // your GeoJSON goes here
var wkt_options = {};
var geojson_format = new OpenLayers.Format.GeoJSON();
var testFeature = geojson_format.read(val);
var wkt = new OpenLayers.Format.WKT(wkt_options);
var out = wkt.write(testFeature);
alert(out);
See it in action: http://jsfiddle.net/FUeJj/
I'm trying to calculate EVI in Arcmap using raster calculator. I put in this calculation
2.5 * ((Float("Band4") - "Band3")) / (Float("Band4") + 6 * "Band3" - 7.5 *"Band1" + 1))
and the output is syntax error. Do you know what I'm doing wrong?
I make a new column in my attribute table and it has a default value (for example 0). I want to (probably using the field calculator) have the result that row 1 has the value 1, row 2=2 and so on, just like a numbered index.
I'd be thankful for a code example for the field calculator (including the use of $rownum (I think I need this for switching to the line)).
Answer
Just put $rownum (QGIS 2) or @row_number (QGIS 3+) as the expression. Simple as that. :)
An up-to-date list of all the field calculator functions can be read in the official docs.
I am attempting to allow for the pop-up to be toggled on and off by pressing a checkbox. Here are some snippets of the code:
var clickHandler
function popupOn(evt){
var query = new esri.tasks.Query();
query.geometry = pointToExtent(map,evt.mapPoint,10);
var deferred = featureLayer.selectFeatures(query,esri.layers.FeatureLayer.SELECTION_NEW);
map.infoWindow.setFeatures([deferred]);
map.infoWindow.show(evt.mapPoint);
}
function clickConnect(connect){
if(connect){
//perform the identify task on click
clickHandler = dojo.connect(map, "onClick", popupOn);
}
else{
//disconnect the click handler
dojo.disconnect(clickHandler);
clickHandler = null;
}
}
We want to display several markers on a static map and want to calculate the optimal zoom-level like Google Maps does. We already calculated the bounding rectangle and the centerpoint of the map but now we have a hard time to calculate the correct zoomlevel to display the whole bounding rectangle. Can someone please point us into the right direction?
I am trying to build my first widget, for use in the Configurable Map Viewer (CMV), using code from ESRI's ArcGIS API for JavaScript feature layer with selection sample. I've been successful in configuring a template and css file so the select and clear buttons show up in the viewer as desired, but I am having trouble modifying the ESRI developer sample for use in the CMV. I've adapted enough of the sample to be able to draw a box on the map with the draw tool but nothing is selected. I'm not following enough of what's going on in some of the other widgets I've looked at to have any confidence that I've successfully adapted the code for PostCreate, initSelectToolbar or sumGasProduction (obviously I haven't). I'm wondering if anyone has tackled the ESRI sample and has any suggestion as to what PostCreate, initSelectToolbar and sumGasProduction should look like or have any suggested references for someone getting started?
Code:
define([
'dojo/_base/declare',
'dijit/_WidgetBase',
'dijit/_TemplatedMixin',
'dijit/_WidgetsInTemplateMixin',
'esri/toolbars/draw',
'esri/InfoTemplate',
'esri/layers/FeatureLayer',
'esri/symbols/SimpleFillSymbol',
'esri/symbols/SimpleLineSymbol',
'esri/Color',
'esri/tasks/query',
'dijit/form/Button',
'dojo/_base/lang',
'dojo/on',
'dojo/text!./Select/templates/Select.html',
'xstyle/css!./Select/css/Select.css'
], function (declare, _WidgetBase, _TemplatedMixin, _WidgetsInTemplateMixin, Draw, InfoTemplate, FeatureLayer, SimpleFillSymbol, SimpleLineSymbol, Color, Query, Button, lang, on, SelectTemplate, css) {
return declare([_WidgetBase, _TemplatedMixin, _WidgetsInTemplateMixin], {
widgetsInTemplate: true,
templateString: SelectTemplate,
selectTools: null,
postCreate: function () {
this.selectTools = new Draw(this.map);
this.selectTools.on('load', lang.hitch(this, 'initSelectToolbar'));
this.fieldsSelectionSymbol =
new SimpleFillSymbol(SimpleFillSymbol.STYLE_SOLID,
new SimpleLineSymbol(SimpleLineSymbol.STYLE_SOLID,
new Color([255, 0, 0]), 2),
new Color([255, 255, 0, 0.5]));
var content = "Status: ${STATUS}" +
"
Cumulative Gas: ${CUMM_GAS} MCF" +
"
Total Acres: ${APPROXACRE}" +
"
Avg. Field Depth: ${AVG_DEPTH} meters";
var infoTemplate = new InfoTemplate("${FIELD_NAME}", content);
this.featureLayer = new FeatureLayer("http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Petroleum/KSPetro/MapServer/1",
{
mode: FeatureLayer.MODE_ONDEMAND,
infoTemplate: infoTemplate,
outFields: ["*"]
});
this.featureLayer.setDefinitionExpression("PROD_GAS='Yes'");
this.featureLayer.setSelectionSymbol(this.fieldsSelectionSymbol);
this.featureLayer.on("selection-complete", 'sumGasProduction');
this.featureLayer.on("selection-clear", function () {
document.getElementById('messages').innerHTML = "No Selected Fields";
});
this.map.addLayer(FeatureLayer);
},
select: function () {
this.selectTools.activate(Draw.EXTENT);
},
clearSelection: function () {
this.featureLayer.clearSelection();
},
initSelectToolbar: function (event) {
//selectionToolbar = new Draw(event.map);
this.selectQuery = new Query();
on(this.selectTools, "DrawEnd", function (geometry) {
this.selectTools.deactivate();
this.selectQuery.geometry = geometry;
this.featureLayer.selectFeatures(this.selectQuery,
FeatureLayer.SELECTION_NEW);
});
},
sumGasProduction: function (event) {
this.productionSum = 0;
//summarize the cumulative gas production to display
arrayUtil.forEach(event.features, function (feature) {
this.productionSum += feature.attributes.CUMM_GAS;
});
document.getElementById('messages').innerHTML = "Selected Fields Production: " + productionSum + " mcf. ";
}
});
});
Answer
The following configuration below seems to be working fine. I'd appreciate any issues being pointed out.
The configuration of the widget:
select: {
include: true,
id: 'select',
type: 'titlePane',
canFloat: true,
path: 'gis/dijit/Select',
title: 'Select',
open: false,
position: 10,
options: {
map: true,
mapClickMode: true,
featureLayerURL: 'http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Petroleum/KSPetro/MapServer/1',
content: "Status: ${STATUS}" +
"
Cumulative Gas: ${CUMM_GAS} MCF" +
"
Total Acres: ${APPROXACRE}" +
"
Avg. Field Depth: ${AVG_DEPTH} meters",
title: '${FIELD_NAME}',
defExpress: "PROD_GAS='Yes'",
clearSelectMsg: 'No Selected Fields',
fieldToSum: 'CUMM_GAS',
selectSumMsgPrefix: 'Selected Fields Production:',
selectSumUnits: 'mcf'
}
},
The widget:
define([
'dojo/_base/declare',
'dijit/_WidgetBase',
'dijit/_TemplatedMixin',
'dijit/_WidgetsInTemplateMixin',
'esri/toolbars/draw',
'esri/InfoTemplate',
'esri/layers/FeatureLayer',
'esri/symbols/SimpleFillSymbol',
'esri/symbols/SimpleLineSymbol',
'esri/Color',
'esri/tasks/query',
'dijit/form/Button',
'dojo/_base/lang',
'dojo/on',
'dojo/_base/array',
'dojo/dom',
'dojo/text!./Select/templates/Select.html',
'xstyle/css!./Select/css/Select.css'
], function (declare, _WidgetBase, _TemplatedMixin, _WidgetsInTemplateMixin, Draw, InfoTemplate, FeatureLayer, SimpleFillSymbol, SimpleLineSymbol, Color, Query, Button, lang, on, arrayUtil, dom, SelectTemplate, css) {
return declare([_WidgetBase, _TemplatedMixin, _WidgetsInTemplateMixin], {
widgetsInTemplate: true,
templateString: SelectTemplate,
selectTools: null,
postCreate: function () {
this.selectTools = new Draw(this.map);
this.selectQuery = new Query();
on(this.selectTools, "DrawEnd", lang.hitch(this, function (geometry) {
this.selectTools.deactivate();
this.selectQuery.geometry = geometry;
this.featureLayer.selectFeatures(this.selectQuery,
FeatureLayer.SELECTION_NEW);
}));
this.fieldsSelectionSymbol =
new SimpleFillSymbol(SimpleFillSymbol.STYLE_SOLID,
new SimpleLineSymbol(SimpleLineSymbol.STYLE_SOLID,
new Color([255, 0, 0]), 2),
new Color([255, 255, 0, 0.5]));
var infoTemplate = new InfoTemplate(this.title, this.content);
this.featureLayer = new FeatureLayer(this.featureLayerURL,
{
mode: FeatureLayer.MODE_ONDEMAND,
infoTemplate: infoTemplate,
outFields: ["*"]
});
this.featureLayer.setDefinitionExpression(this.defExpress);
this.featureLayer.setSelectionSymbol(this.fieldsSelectionSymbol);
this.featureLayer.on("selection-complete", lang.hitch(this, 'sumGasProduction'));
var clearMsg = this.clearSelectMsg;
this.featureLayer.on("selection-clear", function () {
dom.byId('messages').innerHTML = "" + clearMsg + "";
});
this.map.addLayer(this.featureLayer);
},
select: function () {
this.selectTools.activate(Draw.EXTENT);
},
clearSelection: function () {
this.featureLayer.clearSelection();
},
sumGasProduction: function (event) {
var productionSum = 0;
var sumField = this.fieldToSum;
arrayUtil.forEach(event.features, function (feature) {
productionSum += feature.attributes[sumField];
});
dom.byId('messages').innerHTML = "" + this.selectSumMsgPrefix + " " + productionSum + " " + this.selectSumUnits + ". ";
}
});
});
I used to use spyder as python IDE and was learning QGIS developing recently. But I cannot import qgis.core in spyder although the library can be imported in QGIS Desktop Console successfully.
Here is my setting environment:
My OS is Windows 7;
Spyder was integrated in Anaconda2-4.1.1-Windows-x86_64 and install together, in which default interpreter is Python 2.7.12(64bit)
QGIS 2.18 was installed by Standalone installer 64bit but python is 2.7.5 64bit. Then I modified the environment variables by appending
“C:\Program Files\QGIS 2.18\apps\qgis\python” in PYTHONPATH and
“C:\Program Files\QGIS 2.18\apps\qgis\bin;
C:\Program Files\QGIS 2.18\apps\qgis\python;
C:\Program Files\QGIS 2.18\apps\qgis\python\qgis\core;” in Path.
However, I still cannot import it in Spyder’s IPython, error says”DLL load failed:” 
Same error even I change the interpreter to QGIS’s 2.7.5 with Desktop open.
In term of OSGeo4W Shell, that is another problem.
Still failed to import it when finding no module at first opening, but it works only once if run Germán Carrillo’s BAT file. 
Which means that I have to run this batch file every time I start OSGeo4W Shell.
So it is weird my addition path in PATH and PYTHONPATH have no effect at all.
My question is similar to Merging points (bus stops) in QGIS manually or programmatically?, but it's not answered sufficiently
I have OSM-transit data with bus stops. Some stops have two or more bays for different directions. In this example, the stop 'Steinerstraße' has two points: 
I want to merge all points belonging to one stop to one point. I assume that they can be identified by an identical name in the 'name'-column of the attribute table.
Best case would be to create a new point in the middle of the points, however, I'd be satisfied with a simple deletion of the 2nd feature as well.
Can anyone give me a hint how to do it?
I've never used either python nor SQL in a GIS-environment and would prefer a GUI-based solution. However, I'm open for everything (some xp in python and SQL in other contexts).
Answer
Not an easy task, but this way it works for me:
Create Hublines. Select your two layers, and the name attribute for both as ID.After executing, all busstops with the same name are connected with lines. Unfortunately, the points are also connected to themselves, leading to lines of zero length. So:
Vector -> Geoprocessing -> Dissolve with the name field as Dissolve field, to a new fileVector -> Geoprocessing -> Buffer to create a buffer of 5 (meters) in a new layerVector -> Geometry Tools -> Polygon centroids to create single points from the buffer polygonsAs you can see, it works for single and multiple bus stops: 
I have coordinate system with x and y. Also I need to show there x,y point but I have coordinates latitude, longitude.
How can I transform latitude, longitude into x,y to show on x,y scale?
Currently I have a layer with polylines; in another layer, I allowed the user to create 2 points on the polyline.
I want to know if there is a function for the polyline feature that I can use to get the clipped portion of the polyline based on the 2 user-added points. Thanks.
Answer
The testmethod only works with arcmap, but the other methods should work with arcengine.
private void TestGetsubCurve()
{
IPoint pnt1 = null;
IPoint pnt2 = null;
IPolyline polyLine = null;
var gc = ArcMap.Document.FocusMap as IGraphicsContainer;
gc.Reset();
IElement elem;
while((elem = gc.Next())!= null)
{
if (elem.Geometry is IPoint && pnt1 == null)
pnt1 = elem.Geometry as IPoint;
else if (elem.Geometry is IPoint && pnt2 == null)
pnt2 = elem.Geometry as IPoint;
else if (elem.Geometry is IPolyline)
polyLine = elem.Geometry as IPolyline;
}
if (pnt1 == null || pnt2 == null || polyLine == null)
{
Debug.Print("missing some geometry");
return;
}
var subCurve = GetSubCurve(polyLine, pnt1, pnt2);
elem = new LineElementClass();
((ILineElement)elem).Symbol = ((IDocumentDefaultSymbols)ArcMap.Document).LineSymbol;
elem.Geometry = subCurve;
gc.AddElement(elem, 0);
((IActiveView)ArcMap.Document.FocusMap).Refresh();
}
private IPolyline GetSubCurve(IPolyline inpolyLine, IPoint pnt1, IPoint pnt2)
{
double d1 = GetDistAlong(inpolyLine,pnt1);
double d2 = GetDistAlong(inpolyLine, pnt2);
var c = inpolyLine as ICurve;
ICurve outCurve;
c.GetSubcurve(d1, d2, false, out outCurve);
if (c == null || c.IsEmpty)
throw new Exception("unable to get subcurve");
var outPolyline = outCurve as IPolyline;
if (outPolyline == null)
{
// this didn't happen in testing, but one never knows ...
outPolyline = new PolylineClass() as IPolyline;
var sc = outPolyline as ISegmentCollection;
sc.AddSegment((ISegment)outCurve);
((IGeometry)sc).SpatialReference = outCurve.SpatialReference;
}
return outPolyline;
}
private double GetDistAlong(IPolyline polyLine, IPoint pnt)
{
var outPnt = new PointClass() as IPoint;
double distAlong = double.NaN;
double distFrom = double.NaN;
bool bRight = false;
polyLine.QueryPointAndDistance(esriSegmentExtension.esriNoExtension, pnt, false, outPnt,
ref distAlong, ref distFrom, ref bRight);
return distAlong;
}
I have an ArcGIS colormap (*.clr) and I want to open it in QGIS to apply a predefined style.
The *.clr file is really simple, only has four columns representing [value] [red] [green] [blue], as this example:
19 161 161 161
21 152 181 129
22 114 168 144
23 124 142 173
Is there a way to read these columns and convert them in a QGIS style file (*.qml)?
Answer
Nice answer from @whyzar! You could load .clr files into QGIS and then save it as a .qml file. As described in this post, the standard text format is:
Value R G B Alpha Label
So in your case, you could create a text file with:
19,161,161,161,255,19
21,152,181,129,255,21
22,114,168,144,255,22
23,124,142,173,255,23
And load it from the menu:
You can then edit and save the style as a .qml.
From @Stefan's comment, I've included some quick code which reformats the input text file in a way that QGIS can read it (hopefully!):
main_path = 'path/to/directory/'
with open(main_path + 'infile.txt') as infile, open(main_path + 'outfile.txt', 'w') as outfile:
for line in infile:
value = line.split(' ')[0]
line = line.replace(' ', ',').rstrip('\n') + ',255,' + value + '\n'
outfile.write(line)
I want to add only a few specific columns from the join table but there is no drop down list. When selecting manually I'm able to select only 1 column when trying more than 1 I'm not getting the desired results.
I have tried writing like
1] column1, column2
2] 'column1', 'column2'
3] "column1", "column2"
4] (column1), (column2)
5] (column1, column2)
But none of them works. Please help how I can add only few selected columns
PS: the drop down option is available in join attributes by field tool when I try to use it without the processing modeler
How do I get mouse coordinates after a SelectFeature, which returns the feature by itself?
I am trying to project a Lat/Long value I get from a table to a x/y value that I can use to plot a graphic.
I have seen one way to do this may be to use something like this:
feature.Geometry = mercator.FromGeographic(new MapPoint(lat, lon));
feature.Symbol = LayoutRoot.Resources["TRUCK"] as ESRI.ArcGIS.Client.Symbols.Symbol;
graphicsLayer.Graphics.Add(feature);
However my graphics aren't showing up overlapping my tiled basemap service. I tried the mercator.FromGeographic and Mercator.ToGeographic with unsuccessful results.
My Lat/Long value: 37.1223133333333,-93.27893
Mercator.FromGeographic converted value: X=-10382795.4374339, Y=4469266.40579955, WKID=102100
Mercator.ToGeographic converted value: X=-0.000837860806483038, Y=0.000334317600035687, WKID=4326
Neither of these coordinates put my graphics in the correct location. It should be closer to something like this: X=1424565.77, Y=501970.98. I am trying to learn more about coordinate systems, conversions between systems, etc. Any thoughts on what I can try next?
I am still in the beginner phase of learning programming so thanks for your help. Would you have this as a standalone class in your project?
Here is the code I have so far:
C#
using System;
using System.Collections.Generic;
using System.Linq;
using System.Net;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Documents;
using System.Windows.Input;
using System.Windows.Media;
using System.Windows.Media.Animation;
using System.Windows.Shapes;
using ESRI.ArcGIS.Client;
using ESRI.ArcGIS.Client.Geometry;
using ESRI.ArcGIS.Client.Symbols;
using System.Runtime.Serialization;
using ESRI.ArcGIS.Client.Tasks;
using System.Windows.Threading;
using ESRI.ArcGIS.Client.Bing;
namespace TruckGPS_Test
{
public partial class MainPage : UserControl
{
//private static ESRI.ArcGIS.Client.Projection.WebMercator mercator = new ESRI.ArcGIS.Client.Projection.WebMercator();
public MainPage()
{
InitializeComponent();
//AddMarkerGraphics();
QueryTask queryTaskTruckGPS = new QueryTask("http://gisaprd/ArcGIS/rest/services/TruckGPS/MapServer/8");
queryTaskTruckGPS.ExecuteCompleted += QueryTask_ExecuteCompletedTruckGPS;
ESRI.ArcGIS.Client.Tasks.Query queryTruckGPS = new ESRI.ArcGIS.Client.Tasks.Query();
queryTruckGPS.OutFields.Add("*");
queryTruckGPS.ReturnGeometry = false;
queryTruckGPS.Where = "1=1";
queryTaskTruckGPS.ExecuteAsync(queryTruckGPS, "initial");
//This timer updates the outage feature count after a 30 second wait
System.Windows.Threading.DispatcherTimer myDispatcherTimer = new System.Windows.Threading.DispatcherTimer();
myDispatcherTimer.Interval = new TimeSpan(0, 0, 0, 10, 0);
myDispatcherTimer.Tick += new EventHandler(Each_Tick);
myDispatcherTimer.Start();
}
public void Each_Tick(object o, EventArgs sender)
{
QueryTask queryTaskTruckGPS = new QueryTask("http://gisaprd/ArcGIS/rest/services/TruckGPS/MapServer/8");
queryTaskTruckGPS.ExecuteCompleted += QueryTask_ExecuteCompletedTruckGPS;
ESRI.ArcGIS.Client.Tasks.Query queryTruckGPS = new ESRI.ArcGIS.Client.Tasks.Query();
queryTruckGPS.OutFields.Add("*");
queryTruckGPS.ReturnGeometry = false;
queryTruckGPS.Where = "1=1";
queryTaskTruckGPS.ExecuteAsync(queryTruckGPS, "initial");
}
private void QueryTask_ExecuteCompletedTruckGPS(object sender, ESRI.ArcGIS.Client.Tasks.QueryEventArgs args)
{
//If Counter isn't working check the web service, make sure the table is added to the map service
FeatureSet featureSet = args.FeatureSet;
GraphicsLayer graphicsLayer = MyMap.Layers["TruckGPS"] as GraphicsLayer;
graphicsLayer.ClearGraphics();
if (featureSet != null && featureSet.Features.Count > 0)
{
////Check the system time and subtract the time amount wanted from the last gps systime read for movement or idle truck
//System.DateTime dTime = DateTime.Now;
//// tSpan is 0 days, 0 hours, 0 minutes, and 0 seconds.
//System.TimeSpan tSpan = new System.TimeSpan(0, 1, 0, 0);
//System.DateTime result = dTime - tSpan;
foreach (ESRI.ArcGIS.Client.Graphic feature in featureSet.Features)
{
if ((feature.Attributes["LONGITUDE"] != null) && (feature.Attributes["LATITUDE"] != null))
//&& (feature.Attributes["S"])
{
var Machine_Name = feature.Attributes["MACHINE_NAME"].ToString();
var User_Name = feature.Attributes["USER_NAME"].ToString();
double lat = Convert.ToDouble(feature.Attributes["LONGITUDE"]);
double lon = Convert.ToDouble(feature.Attributes["LATITUDE"]);
feature.Geometry = new MapPoint(lat, lon);
//if (User_Name == "CU656")
//{
//feature.Geometry = mercator.ToGeographic(new MapPoint(lat, lon));
feature.Symbol = LayoutRoot.Resources["TRUCK"] as ESRI.ArcGIS.Client.Symbols.Symbol;
graphicsLayer.Graphics.Add(feature);
//}
//else
//{
// feature.Symbol = LayoutRoot.Resources["RedMarkerSymbol"] as ESRI.ArcGIS.Client.Symbols.Symbol;
// graphicsLayer.Graphics.Add(feature);
//}
}
}
{ (MyMap.Layers["TruckGPS"] as GraphicsLayer).Refresh(); };
}
else
{
graphicsLayer.ClearGraphics();
//MessageBox.Show("No features returned from query");
}
}
private void MyMap_MouseMove(object sender, System.Windows.Input.MouseEventArgs args)
{
if (MyMap.Extent != null)
{
System.Windows.Point screenPoint = args.GetPosition(MyMap);
ScreenCoordsTextBlock.Text = string.Format("Screen Coords: X = {0}, Y = {1}",
screenPoint.X, screenPoint.Y);
ESRI.ArcGIS.Client.Geometry.MapPoint mapPoint = MyMap.ScreenToMap(screenPoint);
if (MyMap.WrapAroundIsActive)
mapPoint = ESRI.ArcGIS.Client.Geometry.Geometry.NormalizeCentralMeridian(mapPoint) as ESRI.ArcGIS.Client.Geometry.MapPoint;
MapCoordsTextBlock.Text = string.Format("Map Coords: X = {0}, Y = {1}",
Math.Round(mapPoint.X, 4), Math.Round(mapPoint.Y, 4));
}
}
}
}
At what point would I use the geometry service or call the class you stubbed out?
I'm looking for a complete, freely available file showing UTM Grid zones. Preferably not a simplified version, but one that shows the full grid zones and all 'exceptions' (e.g. in Northern Europe) - as per http://whatutmzoneamiin.blogspot.com/
Does anyone know where I can download such a file?
Answer
Will this do: http://www.baruch.cuny.edu/geoportal/data/esri/world/utmzone.zip I haven't taken a close look at it, but it looked like it has the exceptions on there.
I downloaded raster data sets together with the legend qml files from the SoilGrids page (ftp://ftp.soilgrids.org/). Assignment of the qml legend (TAXOUSDA_1km/TAXN_WRB_1km) didn't work in QGIS 3, as the legend was discretized in floating-point numbers and not integers (to represent the xx unique identifiers of the raster). However, in QGIS 2.18. loading and proper assignment of the legend worked. Thus, I was thinking that this might be a technical issue of the new QGIS version, which could be addressed/solved in the future. Or is there another way to solve the problem in the legend/qml properties?
I am creating a model in gis ModelBuilder. I want to check the type of a particular field and
Is there a way to check the field or any already available python script which could help me do this?
I want my script to perform actions on a particular layer. Online advice I've found so far suggests using the iface method "setActiveLayer()". Makes sense based on the name alone.
Scenario: In the Layers window I have a layer named "Province" highlighted (active). Then I run the following code:
>>> vl = QgsMapLayerRegistry.instance().mapLayer( 'na_roads_Prov' )
>>> iface.setActiveLayer(vl)
Console Output = False
As far as I understand executing this code in the console should move the highlighted "active layer" from "Province" to 'na_roads_Prov'. Yes? Instead the "Province" layer remains highlighted as active, and the console output says "False". Any ideas? Am I not referring to the layer correctly?
Answer
Use mapLayersByName method to get the layer by name
vl = QgsMapLayerRegistry.instance().mapLayersByName('na_roads_Prov')[0]
iface.setActiveLayer(vl)
I'm working on a custom web-based application that generates spatial data and that includes a web-based map viewer. We have a custom servlet that handles WMS requests (not a complete WMS implementation, just enough for our needs) and uses MapServer/SWIG to generate image data. The web front end uses OpenLayers. So far, so good: everything works and gives us nice-looking interactive maps.
We also built a feature that creates a KML file with a
The problem comes when the user starts to tilt the camera and pan around in Google Earth. Google Earth starts mis-registering our data, at first by just a few meters at a time but eventually by tens of kilometers, if you tilt, zoom in and out, and pan around enough.
I did some searches and I found several sites that indicated that Google Earth's support for tilting WMS imagery was poor, but little in the way of specifics. Has anyone else tried this and experienced similar issues?
More importantly, how should I approach this problem? Given that our users are very interested in doing exactly this kind of exploration (plot data with topography, tilt the camera down, and fly around in our subject area) and that our users are already pretty committed to using Google Earth as the client, how should we make the data available so that they can view it?
Answer
I am guessing that you are trying to render a single WMS image to cover the entire view. Remember that when the image is tilted it is very difficult to render a rectangular image of the correct resolution for the entire area. Locations close to the camera need to be rendered differently to features on the horizon. Therefore, direct use of WMS is completely unsuitable for use in Google Earth and should be avoided.
The best approach is to create a network link that divides the area into a set of tiles, each with its own unique link to a fixed extent of the WMS as a Super Overlay. This allows Google Earth to decide which tiles to load and when. It will also make it much simpler to implement caching at the server end. Create larger tiles for different zoom levels using the Region tag. When the view is tilted Google Earth can use lower resolution tiles on the horizon, and higher resolution for closer areas. For an idea of a good structure try converting your data manually using the KML Super Overlay GDAL driver.
For large tiles the WMS response needs to be in WGS84 due to projection issues. For smaller tiles it is sometimes possible to use a local projection system by individually calculating the WGS84 coordinates of the tile corners. There are obvious projection issues, although the results are often better than interpolating.
I have a SRTM DEM and I want to create a shaded relief from it. I created the shaded relief in GRASS and the result is very nice, but a little rough because the area is near flat and the DEM is 90m in resolution.
What I want is to make the DEM smoother in order to generate a smooth shaded relief. Is there an algorithm or interpolation method to do that?
Here is the shaded relief to get an idea, I want to flatten these small bumps:
Answer
How about John Stevenson's r.denoise, from the GRASS AddOns wiki:
r.denoise denoises (smooths/despeckles) topographic data, particular DEMs derived from radar data (including SRTM), using Xianfang Sun's denoising algorithm. It is designed to preserve sharp edges and to denoise with minimal changes to the original data.
I read further from this website (that I also give credit for the above animation) that a more generic method would be to use an Esri ASCII Grid file. The location of mdenoise (downloaded from Sun's website) needs to be in your PATH variable (e.g., Windows users: drop MDenoise.exe in the bin folder with your OSGeo4w or FWTools install). Then, for example, you can use the following shell command to process the ASCII grid file:
# gdal_translate -of AAIGrid my_dem.tif my_dem.asc # convert to .asc
mdenoise -i my_dem.asc -n 5 -t 0.99 -o my_dem_DN.asc # denoise
# gdal_translate -of GTiff my_dem_DN.asc my_dem_DN.tif # convert back to .tif
Denoise is under GNU license, see here
I have geoserver installation that is accessable on both, port 80 and port 8080.
Now, I try to install another instance of geoserver on a different server.
I can access it on port 8080, but not on port 80. I could not find out where to flip the switch to make it run on port 80 as well.
How do I do that?
Answer
Ciao Ferenjito, I would not expose both port 8080 and 80 and above all I would not make Tomcat work on port 80.
My best suggestion is to have something like apache httpd in front of tomcat and have the two talk to each via AJP leaving port 80 in the hands of httpd. Basic/startup docs here.
Hope that helps.
I'm trying to convert a map I found here into lat/lon compatible with another data set I created from Google Maps lat/lons.
From the directory with that .shp file, I run:
library(rgdal)
SG = readOGR('../data', 'MP14_PLNG_AREA_NO_SEA_PL')
A subset of the data I'm trying to co-plot is:
dim = c(5L, 2L)
# polygon coordinate matrices
poly_mats = list(structure(c(103.843, 103.843, 103.854, 103.854, 103.843,
1.318, 1.324, 1.324, 1.318, 1.318), .Dim = dim),
structure(c(103.843, 103.843, 103.854, 103.854, 103.843,
1.324, 1.329, 1.329, 1.324, 1.324), .Dim = dim),
structure(c(103.854, 103.854, 103.865, 103.865, 103.854,
1.318, 1.324, 1.324, 1.318, 1.318), .Dim = dim),
structure(c(103.854, 103.854, 103.865, 103.865, 103.854,
1.324, 1.329, 1.329, 1.324, 1.324), .Dim = dim))
# nest properly to create SP object
polySP = SpatialPolygons(lapply(seq_along(poly_mats), function(ii) {
Polygons(list(Polygon(poly_mats[[ii]])), ID = ii)
# assign lat/lon CRS
}), proj4string = CRS('+init=epsg:3857'))
Both of these objects look fine when plotted separately:
png('separate.png')
par(mfrow = c(1, 2))
plot(SG)
plot(polySP)
dev.off()
Not that it should matter, but for certain, the objects are currently in different CRS:
identical(proj4string(SG), proj4string(polySP))
# [1] FALSE
So we should be able to transform SG like so:
SG = spTransform(SG, proj4string(polySP))
But this fails:
png('together_no.png')
plot(SG)
plot(polySP, col = 'red', add = TRUE)
dev.off()
It seems spTransform has failed, as the coordinates in the output are non-sensical as lat/lons:
coordinates(SG)[1:5, ]
# [,1] [,2]
# 0 11559649 153646.0
# 1 11569258 147405.4
# 2 11559462 150848.3
# 3 11544079 146374.9
# 4 11549925 150925.5
What's going wrong/how can this be rectified?
Answer
You've created your polygons with epsg:3857:
proj4string = CRS('+init=epsg:3857'))
which is Google's Web Mercator, not lat-long. You mean epsg:4326, WGS84 lat-long, most likely.
Some more detailed reading:
Client: OpenLayers / Server: Geoserver
When I update a WFS Element,
JAVASCRIPT:
element.feature.attributes.ENAME = "Mystring";
saveStrategy.save();
HTTP-POST:
ENAME
MyString
some features of Type String do not get fully updated, only the first letter ends up in the WFS layer.
...
M
//only M ends up on the server, not Mystring as expected
the respective line in the logfile (catalina.out) says:
name = ENAME
value = Mystringproperty[0]:
strangely, the attached property[0] is on all attributes, not only on the one which is reduced to its first letter, but also on some strings that do not get shortened as well as some integers.
I would expect the type of the respective feature to be something like String[1], but as far as I can see the Type System, there are no length attributes to Strings.
Who knows what I can do?
UPDATE
When I change the value of any feature's attribute (ENAME) with OpenJump and upload the changed data, my application shows the feature correctly, which means, the geoserver itself CAN store a string longer than 1 in the respective attribute.
However, everytime I make an update via OpenLayers, regardless on which feature, any other feature's ENAME attribute gets shortened (like from 'MyString' to 'M'), even if it was a totally different feature that was updated on OpenLayers.
The log (catalina.out) says NOTHING about this perfidious update operation.
UPDATE
Content of the Firebug Post TAB, Content of the respective Log Section:
2015-01-20 17:01:00,671 INFO [geoserver.wfs] - Request: transaction service = WFS version = 1.1.0 baseUrl = http://myserver.de:80/geoserver/ group[0] = wfs:update=net.opengis.wfs.impl.UpdateElementTypeImpl@53a7e8d6 (filter: [ Graeber.203 ], handle: null, inputFormat:
STILL not solved
WFS Layer:
Javascript definition:
layer_selectable = new OpenLayers.Layer.Vector("WFS", { strategies : [new OpenLayers.Strategy.BBOX(), saveStrategy], styleMap: new OpenLayers.StyleMap({ 'default': style, 'select': style_selected }), protocol : new OpenLayers.Protocol.WFS({ url : geoserver + config.wfs_layer, version : "1.1.0", featureType : "Graeber", featureNS : config.featureNS, srsName : epsg, }), renderers : renderer //tileOptions: {crossOriginKeyword: 'anonymous'} });
Answer
You could edit your WFS-T Service directly with QGIS for example to see if it's an error of your javascript-application or of geoserver itself.
http://docs.qgis.org/2.2/en/docs/user_manual/working_with_ogc/ogc_client_support.html section "WFS and WFS-T Client" ( or in german: "QGIS als OGC Datenclient")
In QGIS you could also check which limit is set for the character-length of the column that gets cut off (in Layer-Properties):
Have you considered importing the shapefile to postgis and then publish the layer directly from an postgis-store?
If you like to you can post the javascript-code on http://jsfiddle.net/ or http://www.codeshare.io/
The Content of the Firebug Post TAB and Content of the respective Log Section you posted shows different versions(1.0.0 vs. 1.1.1):
2015-01-20 17:01:00,671 INFO [geoserver.wfs] -
Request: transaction
service = WFS
version = 1.1.0
but I guess these were two different tries, weren't it?
