texts: Using postgis to rasterize a set of points grouping by time

texts:
I have used postgis for vector operations for years, but have recently begun a project mixing vector and raster data and I am bit confused as to what is the best way to do what I want to do. I have the trajectory of a boat as lon,lat points that I want to rasterize, counting the number of points in each grid cell. The trick is that I want to groups the points according to time, counting a set of points only once if they all occur in the same grid cell within the same time window (e.g., within an hour of each other). Efficiency is important because in reality I have many many such trajectories each with a large number of points.

I can think of several ways of attacking this problem, a couple of which I detail below, but in all cases there is a step that I am not sure how to do. More generally, I am wondering if there is a standard approach to these types of problems.

Approach 1:

Group the points of a trajectory by time using ST_Union grouping over a truncated version of the timestamp for each point
Use ST_AsRaster on each of these groups of points to produce one raster layer for each group with 1's wherever the points intersect the reference raster layer.
Sum together these individual raster layers to get the final result.

My problems with this approach are (1) that I am not sure how to sum all the rasters together and (2) it seems inefficient to me to create all these raster layers. Regarding the first problem, I believe I can use ST_SummaryStatsAgg to do this (though I would need to update my version of postgis to have this function), but an example would help. Regarding the second problem, I will have many many time periods and trajectories and calculating all of these rasters and then computing a full suite of summary statistics over these rasters seems daunting.

Approach 2:

This approach essentially passes through vector calculations only returning to raster maps at the end.

Convert the raster layer to a vector layer using ST_PixelAsPolygons (adding a spatial index to these polygons)
Intersect these polygons with the points to associate a grid cell with each point.
Group the points by time period and grid cell so that we have a table with one row per time period and intersected grid cell combination.
Aggregate across time periods so that we count the number of "positive" (boat present) time periods per grid cell.
Convert these counts per grid cell back into a raster layer.

The problem here is that I am not sure how to do the last step. I can image using ST_AsRaster on the centroid of each grid cell to generate a raster with one non-zero pixel, then using ST_MapAlgebraExpr to multiply this 1-pixel raster by the counts, an then using ST_SummaryStatsAgg to put all these 1-pixel rasters back together, but there must be a better way.


texts: How to extract speed from gpx with "track points badly ordered"

texts:
I'm trying to convert gpx (without speed information) to csv (with speed information). I have gpx file that contains multiple <trk> tags, each one contains one trkseg with track points. All track points seem to be correctly ordered (from quick scrolling trough it). (Data are recorded by human, so there are segments probably because human paused recording.)

I have used following command from this answer:

gpsbabel -t -i gpx -f input.gpx -x track,speed -o unicsv -F output.csv

I got following error:

trackfilter-init: Track points badly ordered (timestamp 2008-10-09T16:16:06Z > 2008-10-09T16:15:55Z)!

Is there a way to automatically fix this? I don't necessarily need to process file as multiple segments, I need to dump everything to one csv.


texts:
Just installing Qgis 1.8 Lisboa, but I can't find "Quick Print" Extension.
Any idea where to catch it again?


texts: CartoDB - How do you do a line fill of polygons in CartoDB?

texts: I would like to show overlaping layers with hash marks, is this possible by editing the css in some way?


texts: Programmatically Selecting CRS in QGIS using Python?

texts:
When you run the following Python code in the Python Console of QGIS 2.4,

uri = "file:///C:/data.csv?type=csv&xField=lng&yField=lat&spatialIndex=no&subsetIndex=no&watchFile=no"
vlayer = QgsVectorLayer(uri, 'Points', "delimitedtext")

you are presented with the Coordinate Reference System Selector popup window as shown below.

Is there a way to programmatically select a CRS instead of having a user manually select it?

Tried the following but the CRS Selector dialog window still appears.

uri = "file:///C:/data.csv?type=csv&xField=lng&yField=lat&spatialIndex=no&subsetIndex=no&watchFile=no"
vlayer = QgsVectorLayer(uri, 'Points', "delimitedtext")
vlayer.setCrs( QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId) )


texts: Convert downloaded LiDAR info to DEM in GRASS

texts:
I have downloaded LiDAR information as a .asc file and a .tif file, and even though I have tried multiple methods to convert them to a dem I am having no luck. Grass GIS says there is an invalid northpoint on the ascii data. Could someone here point me in the right direction?

Here are the two datasets: 1 and 2.


texts: How to close QGIS polygon without a mouse?

texts:
One of my students is using a Windows laptop that does not have physical trackpad buttons; left and right clicks are completed by one and two-finger taps, respectively. The one/two-finger taps work for everything else in QGIS except closing polygons—when a two-finger tap is made (i.e., a right click), the polygon simply disappears. Is there a key command that can be used instead?

This is QGIS 3.10.2 on Windows 10.


texts:
I would like to export a plan from MapInfo into DWG format. This exported DWG will be used by a laser cutter to engrave the vector outlines into wood.

I'm aware of the Universal Translator tool included within MapInfo but this tool only converts a single layer to an individual DWG file. I would like to complete the same process but convert a complete Plan view into  a single DWG file.

Is there a way in MapInfo to export a full plan as DWG, if not what other combinations of software could I use to create the same output?


texts: Add geometry attributes tool in ArcMap successfully runs but no fields are added to attribute table

texts:
I am attempting to add a line bearing attribute column to my shapefile in ArcMap. I'm using the Add Geometry Attributes tool to accomplish this. Running this tool is successful, however when I inspect the attribute table, no new fields are added. If I repeat the process, the tool runs successfully again, but this time with a warning:

  WARNING 001097: BEARING already exists -- overwriting.

Again, there is no new attribute table field for line bearing.

This is not an issue to my specific dataset. Any polyline shapefile I try this on gets the same result.

If you have experienced this problem, how did you solve it?


texts: How to include the metadata in a shapefile download from geonetwork?

texts:
My web colleagues have set up geoserver and a geonetwork data portal website for which I am creating content.

I have used FTP to upload a shapefile to geoserver including its xml metadata file.
Created a store for this dataset on geoserver.
Gone into geonetwork, imported metadata to create an entry in geonetwork.
To enable the dataset to be downloaded I have added an online resource as shown below:

This creates a download link on the main metadata page, which if you click on download, downloads the shapefile as a zip file.

If one opens the zip file you see something like this:

There is no XML file which seems a bit crazy to me as now you have distributed a dataset without any metadata which is obviously a bad thing!

I can't work out how one ensures that the XML (which I had uploaded using FTP along with all the other files that make up a shapefile) is included. I've tried to search the geonetwork online help and failed to find a solution. Does anyone know how to do this?


texts: Why different transformation behaviour ArcGIS Pro 2D/3D layers?

texts:
I am working with ArcGIS Pro and have noticed a strange difference between the co-ordinate transformation of 2D and 3D layers.

Here are my transformation settings:

I have the same layer represented in 2D and 3D but you can see they do not align:

3D layer: apparently no transformation

2D layer: correct transformation

Has anyone else experienced this issue before?


texts: GeoWebCache 1.11.0 application/json not supported for GetFeatureInfo

texts:
I have been trying to find out whether or not I can serve a response to a GetFeatureInfo request from my GeoWebCache in a JSON format, but I have been coming up empty. My capabilities list doesn't include it, but I don't see why GeoWebCache wouldn't support it when Geoserver does.

My work around consists of altering the URL of the GetFeatureInfo call to point to my Geoserver instead of my GeoWebCache, but it feels a bit hacky if it's possible to get the JSON format through GWC.

I'm using Openlayers 4 and the built-in WMS GetFeatureInfo URL functionality.


texts:
Does anyone have the proper codes for line simplification and smoothing algorithms which is good to implement instead of ArcGIS algorithms?

For example any of these:

Visvalingam–Whyatt
Reumann–Witkam
Opheim simplification
Lang simplification
Zhao-Saalfeld


texts: Publishing WFS and Feature Access through ArcPy script?

texts:
I'm trying to publish a large number of MXD files to my ArcGIS Server with an ArcPy script. Here's the overview of how I'm doing this

mapDoc = arcpy.mapping.MapDocument(...)
con = r'C:\...'
service = currentMXD[:-4]
sddraft = wrkspc + '\\sdfiles\\' + service + '.sddraft'
sd = wrkspc + '\\sdfiles\\' + service + '.sd'
arcpy.mapping.CreateMapSDDraft(mapDoc, sddraft, service, 'ARCGIS_SERVER', con, True, 'etc')
arcpy.mapping.AnalyzeForSD(sddraft)
arcpy.StageService_server(sddraft, sd)
arcpy.UploadServiceDefinition_server(sd, con)

Pretty much following the example code of how to publish using an MXD file. Now,
I want to add WFS and Feature Access capabilities when publishing these MXDs, which seems like it should be simple because you just need to check the boxes off when publishing manually, but I can't find how to do so. Does anyone know of a way to do so?


texts: What's the difference between gdalwarp and gdal_merge for mosaicing

texts:
I wanted to merge around 20 dem tiffs each around 100MB, and I used gdalwarp (I know about virtual file option):

gdalwarp $(list_of_tiffs) merged.tiff

However this command took very long time, and after some 15min it was processing 5th image from set. I terminate it.

I then found that gdal_merge.py script can be used also for mosaicing, and tried it:

gdal_merge.py $(list_of_tiffs)

which completed in less then 3min.

As I expected that both commands would produce same result, I wonder what's the difference between the two, as to why gdalwarp takes so much time if result is same?


texts: Setting multiple colors transparent in raster map in ArcMap?

texts:
I have an aerial scene in jpeg2000 Format in which I would like the black and white areas to be transparent:

In QGIS I can set transparent pixels in the layer properties by simply adding colors picked from the map to the transparency percentage I desire:

In ArcGIS Desktop (10.6) at first sight in the layer properties I can only define one color (here: either black or white) that is displayed transparent:

Is there an Option in ArcGIS Desktop to display more than one color transparent?

I do not have the Spatial Analyst extension licensed.


texts: Create attributes from unique values of another features classes category attribute with FME

texts:
I have a table (say landuse polygons for example) structured as follows:

ID | landuse_category
---------------------
1  | woodland
2  | industry
3  | farmland
4  | woodland
5  | lake
6  | farmland
7  | industry
8  | lake

Now I would like to create attributes from the distinct values (farmland,industry,lake,woodland) of landuse_category to a table districts for example like

ID | farmland | industry | lake | woodland

I succeeded in obtaining the distinct values with a DuplicateFilter with key attribute landuse_category on the Unique port and then create a list of these unique values with a ListBuilder like

Since input data may vary, this is a generic task and thus the attribute values themselves should not be part of the process.

How to create attributes on another table using this list?


texts: PyQGIS - Print Layout Legend. setLegendFilterByMapEnabled()

texts:
I'm creating a legend in my print composer and only want "visible", or active, layers shown in the legend. This means I want only layers with the box checked in my layers panel.

Here is how I create my legend.

legend = QgsLayoutItemLegend(layout)
legend.setTitle("Legend")
layout.addLayoutItem(legend)

Looking at the documentation for QgsLayoutItemLegend there is a method: setLegendFilterByMapEnabled() whose description is:

This leads me to believe if I do something like:

legend.setLegendFilterByMapEnabled(True)

My legend will only include "active" layers. But my legend still includes all layers, active and inactive.

What am I missing?


texts: ST_DISTANCE isn't giving the expected result

texts:
I don't get why ST_DISTANCE is giving me the same closest node when clearly it is not the one nearest to the starting point. Also, using 2 different starting points, I still have the same result.

My data is this one, I run in my pgadmin

SELECT id, ST_AsText (ST_Transform (the_geom, 4326)) FROM edges_noded order by id;

and I run this commands to get the closest node

SELECT id,
       ST_Distance(the_geom,ST_SetSRID(ST_MakePoint(124.263770, 8.265921),3857)) AS DIST,
       ST_AsText (ST_Transform (the_geom, 4326))
FROM edges_noded
LIMIT 1;

SELECT id,
       ST_Distance(the_geom,ST_SetSRID(ST_MakePoint(124.264505, 8.264620),3857)) AS DIST,
       ST_AsText (ST_Transform (the_geom, 4326))
FROM edges_noded
LIMIT 1;

these two commands give me the same nearest node.

Does ST_DISTANCE work for this kind of situations?

attached also is a visual of the situation.


texts: How to get the Features of FeatureCursor?

texts:
I am unsure how to correctly use IFeatureCursor.NextFeature method.

I've selected a few features:

I need to assigne the selected features into few pFeature variables (pFeature1, pFeature2, pFeature3)

The code below seems to be the solution, but it doesnt.

Dim pMxDoc As IMxDocument
Dim pFeatureLayer As IFeatureSelection
Dim pFeatureSelection As ISelectionSet
Dim pFeatureCursor1 As IFeatureCursor

Dim pFeature1 As IFeature
Dim pFeature2 As IFeature
Dim pFeature3 As IFeature

Set pMxDoc = ThisDocument
Set pFeatureLayer = pMxDoc.FocusMap.Layer(0)
Set pFeatureSelection = pFeatureLayer.SelectionSet

pFeatureSelection.Search Nothing, True, pFeatureCursor1

Set pFeature1 = pFeatureCursor1.NextFeature
Set pFeature2 = pFeatureCursor1.NextFeature
Set pFeature3 = pFeatureCursor1.NextFeature

MsgBox pFeature1.Value(2) & pFeature2.Value(2) & pFeature3.Value(2)

To check the result I used a MsgBox containing the value of "NAME_MST" field.

I can't explain it, but every PFeature contains the last Feature of pFeatureCursor.

To solve this problem I have to write a preposterous code as shown below (by using a few pFeatureCursors).

Dim pMxDoc As IMxDocument
Dim pFeatureLayer As IFeatureSelection
Dim pFeatureSelection As ISelectionSet
Dim pFeatureCursor1 As IFeatureCursor
Dim pFeatureCursor2 As IFeatureCursor
Dim pFeatureCursor3 As IFeatureCursor

Dim pFeature1 As IFeature
Dim pFeature2 As IFeature
Dim pFeature3 As IFeature

Set pMxDoc = ThisDocument
Set pFeatureLayer = pMxDoc.FocusMap.Layer(0)
Set pFeatureSelection = pFeatureLayer.SelectionSet

pFeatureSelection.Search Nothing, True, pFeatureCursor1
pFeatureSelection.Search Nothing, True, pFeatureCursor2
pFeatureSelection.Search Nothing, True, pFeatureCursor3

Set pFeature1 = pFeatureCursor1.NextFeature

Set pFeature2 = pFeatureCursor2.NextFeature
Set pFeature2 = pFeatureCursor2.NextFeature

Set pFeature3 = pFeatureCursor3.NextFeature
Set pFeature3 = pFeatureCursor3.NextFeature
Set pFeature3 = pFeatureCursor3.NextFeature

MsgBox pFeature1.Value(2) & pFeature2.Value(2) & pFeature3.Value(2)

The result is correct but I don't like the way I've reached it. What do I do wrong in first code?


texts: Why are my DMS points plotted way off from their true locations?

texts:
I'm having a bit of a nightmare plotting some GPS coordinates in QGIS. My coordinates came from a PDF so it took a little wrangling to extract the from the PDF and convert them to a CSV file. I've done this now, and have ended up with a CSV containing the coordinates in DMS format. When I attempt to plot in QGIS (add delimited layer, using DMS option) they are plotted way off from their true locations. Is it possible that the coordinates were captured in a different projection to that which they are initially plotted in QGIS (WGS84)?

Or perhaps the coordinates were just incorrect in the first place. For example, If I copy and paste Northing and Eastings into Google Maps the resulting point is about 1000 miles away from where I would expect it to be.

Anyone have any experience of importing GPS coordinates into QGIS with similar results?


texts: Displaying raster images vertically in QGIS 3.x (Geological cross section)

texts:
I'm aware this has been asked before to an extent, but the answer here doesn't really cover raster images.

I want to 'hang' a raster image in 3D alongside a surface. The raster is a geological cross section, and the surface is topography. You can do this in ArcScene rather hackily, I was wondering if there is the possibility of doing something similar in QGIS.

Just to be clear the cross sections are hand drawn and there is no actual data to sample or interpolate in 3D to make a section. I'm not slicing a 3D model, as the area around the cross section is almost 100% uncertain.

If you can extrude between shapefiles/vectors I could digitise the cross section lines and extrude them but I've no idea how to do this in QGIS and produce something that looks publication ready and better than simply displaying the sections in 3D space.

How can I do this?


texts: Recognizing when road is uphill and when downhill using QGIS?

texts:
How do I automatically recognize and cut roads of a great network in places when the roads are changing the path uphill and downhill using QGIS?

I have a geotiff dem.

Schematic truck's transport


texts:
Using QGIS with GRASS I have found a route through r.drain but it has many sharp curve. the maximum radius of curve allowed to me is 176m.
Is there any to reduce the curve sharpness to 176m or at least to find the radius of a curve?


texts: Print Composer: How to use Print Composer Atlas to swap between columns of linked data

texts:
I have a map with an openstreetmap layer and a point-layer, in which linked data (rows of a csv-file) is classified in 5 classes (the values are between 1 and 100).
The linked data (csv-file) looks like this:
HEADER: name indicator1 indicator2 indicator3 ..... indicator50
Row1:   A    100        75         50.7       ..... 65
Row2:   B    80         65         24         ..... 43
.                         .                           .
.                         .                           .
.                         .                           .
R1000:  Z    20         25         94         ..... 73

The layer contains coordinates for points with the names (A,B,...,Z) (i joined the corresponding fields in layer>properties>joins) and is linked to a selected column (indicator1, indicator2 etc), so based on which column it is linked to, the points are colored following the class borders (e.g. the point with the name A would be displayed green if linked to indicator1 or yellow if linked to indicator3).
The column of the linked csv-file can be chosen manually in layer>properties>symbology>column, what I want to do now is to use the print composer atlas-function to do this automatically.
So basically what I want is an automized output for each indicator, so I want the content of the map (the colors of the points) to swap based on the page of the atlas. I only found out that I can apply the atlas function on the rows of the indicator that I choose before I create a new layout for the map. So the map will just jump from point A-Z based on the pages but its all linked to one specific indicator.
Is there any way to change the value of a layer (which I mentioned can be done manually in layer>properties>symbology>column) with the atlas function or is there another way to do this?


texts: Reading response of GetFeature request with OpenLayers

texts:
I am using OpenLayers to send a GetFeature request to GeoServer. With the help of Firebug I can see that the response contains some of the attribute data I need. I am interested in showing some attributes once the user clicked on a feature.

How can I make sure that GeoServer sends all attributes and how can I read them in OpenLayers?

Here is my code, for now I only succeeded at displaying the feature selected on my map using an OpenLayers example:

Feature = new OpenLayers.Layer.WMS(
        "Ramsar sites",
        "http://ramsardev.grid.unep.ch:8080/geoserver/wms",
        {
            layers:'ramsar_sdi:features',
            transparent:'true',
            visibility: true,
            isBaseLayer: false
        });
map.addLayer(Feature);

select = new OpenLayers.Layer.Vector("Selection", {styleMap: new OpenLayers.Style(OpenLayers.Feature.Vector.style["select"])});
map.addLayer(select);

control = new OpenLayers.Control.GetFeature({
        protocol: OpenLayers.Protocol.WFS.fromWMSLayer(Feature),
        box: false,
        hover: false,
        multipleKey: "shiftKey",
        toggleKey: "ctrlKey"
    });

control.events.register("featureselected", this, function(e) {
        select.addFeatures([e.feature]);

    });

control.events.register("featureunselected", this, function(e) {
        select.removeFeatures([e.feature]);
    });

map.addControl(control);
control.activate();

The problem is that I can't see any of my attributes anywhere in the event object e when the event featureselected is triggered.


texts: I have a dwg file containing some polylines that is georeferenced in UTM coordinates, but unfortunately has been converted to mm instead of metres, so has 3 extra digits on each x and y point. Is there a way to move the decimal place back by 3 and convert to GIS format? I'm using FME.


texts: How to override the memory limits in GeoServer cross layer filtering?

texts:
According to this document retrieving of minimum features and coordinates can be override by using QUERY_LAYER_MAX_FEATURES and
GEOMETRY_COLLECT_MAX_COORDINATES. This can be done using JVM system variables, servlet context variables, or environment variables.

My doubt is where to set these values. Can any one explain briefly?
What is meant by JVM system variables?


texts: Passing selection through relates in ArcGIS ModelBuilder?

texts:
I have a layer of facilities and a feature that has three related layers (Point, Area, Line). I'm trying to automate extracting data from the Line layer that's related to a selected set of facilities via relationship classes. I'm using Arcmap 10.3.

When I do this manually, I select all features in Facilities, open its attribute table then select the Facilities / Point relationship from a set of saved relationships. This opens the Point attribute table, which shows selected records related to the selected items in Facilities. Next, I select the Point / Area relationship, which opens the Area attribute table, also showing selected records. Finally, I select the Area / Line relationship. When that attribute table opens, it shows a selection of Line features that are related back to the selected Facilities features. From here, I can export the selected data from the Line layer.

I've created a model that runs successfully and creates all three relationships. The records in Facilities are selected but when I open the Line attribute table, none of them are.

What am I doing wrong?


texts: Editing a map created with PyQgis3 in Qgis Print Layout

texts:
I create a new print layout using PyQgis and then I open it in the layout designer. Everything seems ok till I try to refresh the page or update the map extent, then map disappears.

project = QgsProject.instance()
lmg = QgsProject.instance().layoutManager()
l = QgsPrintLayout(project)
l.initializeDefaults()
l.setUnits(QgsUnitTypes.LayoutMillimeters)
lmg.clear()
lmg.addLayout(l)
page = l.pageCollection().pages()[0]

lm = 20         # left margin
tm = 32         # upper margin

w, h = 900, 600 #width and height
page.setPageSize(QgsLayoutSize(1189, 841))
theMap = QgsLayoutItemMap(l)
theMap.updateBoundingRect()
theMap.setRect(QRectF(lm,tm, w, h))
theMap.setPos(lm,tm)
theMap.setFrameEnabled(True)
theMap.updateBoundingRect()
theMap.setLayers(project.mapThemeCollection().masterVisibleLayers())
theMap.setExtent(QgsRectangle(363600.0, 376800.0, 364500.0 ,377400.0))
theMap.attemptSetSceneRect(QRectF(lm,tm, w, h))
l.addItem(theMap)
l.updateBounds()

iface.openLayoutDesigner(l)


texts: Extracting raster values at points locations in R, but resulting in a lot of NA value

texts:
I use point data with 28 points, to extract the values from 3 raster layers, but got a lot of NA value.
The code is like:

#read the point and change its coordination system to fit the raster layers
file<-paste("/Volumes/Distribution/point0.csv")
file
Point1<-read.table(file, header=TRUE, sep=",")
sppoints<-SpatialPoints(Point1, proj4string=CRS("+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
slope1<-brick("/Volumes/rasters/slope.tif")
slope1
tp<-spTransform(sppoints, crs(slope1))
tp

#Read the raster layers
files<-list.files("/Volumes/rasters", pattern="tif$", full.names=TRUE)
files
predictors<-stack(files)
predictors

#Extract values
presvals<-extract(predictors, tp)
set.seed(0)
backgr<-randomPoints(predictors, 500)
absvals<-extract(predictors, backgr)
pb<-c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
sdmdata<-data.frame(cbind(pb, rbind(presvals, absvals)))
head(sdmdata)
tail(sdmdata)

Then I got head(sdmdata) results like:

  pb aspect elevation slope
1  1     NA        NA    NA
2  1     NA        NA    NA
3  1     NA        NA    NA
4  1     NA        NA    NA
5  1     NA        NA    NA
6  1     NA        NA    NA

and tail(sdmdata) results like:

    pb         aspect     elevation          slope
523  0 213.0868417062 3.00000000463  4.90379238129
524  0 294.1608609572 4.00000000000 10.47768211365
525  0 341.7759552833 7.00000000000 24.43413352966
526  0 246.7905575610 7.00000000000 29.83376884460
527  0 247.4420388641 3.99999999843 28.22425270081
528  0  23.7686386468 3.99999999749  4.22268724442

I plotted the points and raster layers in QGIS, the points fell exactly in the extent of these raster layers. I don't know why I got so many NA values.


texts: What is a good formula for calculating if a LatLong is within a radius (in miles) of a second LatLong?

texts: If I have two LatLongs, what is the best way to determine if the first one is within a 10 mile radius of the second?


texts:
I need to create a buffer around a polygon.

I know how to do this, but my problem is I need the buffer to be 250 feet to the south, east, and west, but 500 feet to the north.

Is there any way to specify the distance of a buffer in a certain direction from a feature?


texts: How to Update Multiple Fields of one Layer Using WFS Service

texts:
Trying to update two fields of a shapefile using a WFS update service published from GeoServer 2.4.3 with OpenLayers 2.13.1. I have a JavaScript with all the OpenLayers fun and a HTML document with all the interface fun. Essentially the user clicks on an object and updates a form and then saves it.

Assuming its possible? It's probably really simple but I can't crack it. The fields do not update and in the demo requests the service times out.

Here is the code I am trying;

function saveHandler() {
newDim = document.getElementById("dimension").value;
var num = parseFloat(newDim);
var newDimRound = num.toFixed(3);
newMaterial = document.getElementById("MaterialDD").value;

var saveDetails = OpenLayers.Request.POST({
        url: '<url of wfs service>',
        username: '<username>',
        password: '<password>',
        data: '<wfs:Transaction service="WFS" version="1.1.0"\n' +
                        'xmlns:topp="http://www.openplans.org/topp"\n' +
                        'xmlns:ogc="http://www.opengis.net/ogc"\n' +
                        'xmlns:wfs="http://www.opengis.net/wfs">\n' +
                        '<wfs:Update typeName="Store:Layer">\n' +
                        '<wfs:Property>\n' +
                        '<wfs:Name>Field 1</wfs:Name>\n' +
                        '<wfs:Value>' + newDimRound + '</wfs:Value>\n' +
                        '<wfs:Name>Field 2</wfs:Name>\n' +
                        '<wfs:Value>' + newMaterial + '</wfs:Value>\n' +
                        '</wfs:Property>\n' +
                        '<ogc:PropertyIsEqualTo>\n' +
                        '<ogc:PropertyName>Spatial_ID</ogc:PropertyName>\n' +
                        '<ogc:Literal>' + getSpatial_ID + '</ogc:Literal>\n' +
                        '</ogc:PropertyIsEqualTo>\n' +
                        '</wfs:Update>\n' +
                        '</wfs:Transaction>\n',
        callback: handler
    }); // end of Post

    function handler(updateSQL) {
        alert("Details Saved!")
    } // end handler
};

I have tried putting the fields in their own property sections and different orders to no avail.


texts: Example for Coordinate Transformation Suite (abridged CTS) Android

texts:
I want to transform my EPSG:25833 coordinates to google WGS84 at my android app.
After long time searching i found the Coordinate Transformation Suite (abridged CTS) library.

Has anybody a working example

if i try the example from orbisgis i get the error Caused by: java.lang.ClassNotFoundException: Didn't find class "org.slf4j.LoggerFactory" on path:

   CRSFactory cRSFactory = new CRSFactory();

    // Add the appropriate registry to the CRSFactory's registry manager. Here the EPSG registry is used.
    RegistryManager registryManager = cRSFactory.getRegistryManager();
    registryManager.addRegistry(new EPSGRegistry());

    // CTS will read the EPSG registry seeking the 4326 code, when it finds it,
    // it will create a CoordinateReferenceSystem using the parameters found in the registry.
    CoordinateReferenceSystem crs = cRSFactory.getCRS("EPSG:25833");
  //  CoordinateReferenceSystem crs2 = cRSFactory.getCRS("EPSG:4326");

    double[] coord = new double[2];
    coord[0] = 54.321;
    coord[1] = 9.876;

    CoordinateReferenceSystem  crs1 = cRSFactory.getCRS("EPSG:4326");
    CoordinateReferenceSystem  crs2 = cRSFactory.getCRS("EPSG:25832");
    GeodeticCRS sourceGCRS = (GeodeticCRS) crs;
    GeodeticCRS targetGCRS = (GeodeticCRS) crs2;
    List<CoordinateOperation> coortzdOps = CoordinateOperationFactory.createCoordinateOperations(sourceGCRS,targetGCRS);
    if (coortzdOps.size() != 0) {
        for (CoordinateOperation op : coortzdOps) {
            double[] dd  = op.transform(coord);
            for (int i = 0; i < dd.length; i++) {

            }
        }
    }

Error log:

  java.lang.NoClassDefFoundError: Failed resolution of: Lorg/slf4j/LoggerFactory;
                            at org.cts.registry.RegistryManager.(RegistryManager.java:40)
                            at org.cts.CRSFactory.(CRSFactory.java:68)
                            at com..transform.trans(transform.java:169)
                            at com..MapGoogle.onMapReady(MapGoogle.java:75)
                            at com.google.android.gms.maps.SupportMapFragment$zza$1.zza(Unknown Source)
                            at com.google.android.gms.maps.internal.zzp$zza.onTransact(Unknown Source)
                            at android.os.Binder.transact(Binder.java:392)
                            at xz.a(:com.google.android.gms.DynamiteModulesB:82)
                            at maps.ad.u$5.run(Unknown Source)
                            at android.os.Handler.handleCallback(Handler.java:815)
                            at android.os.Handler.dispatchMessage(Handler.java:104)
                            at android.os.Looper.loop(Looper.java:207)
                            at android.app.ActivityThread.main(ActivityThread.java:5683)
                            at java.lang.reflect.Method.invoke(Native Method)
                            at com.android.internal.os.ZygoteInit$MethodAndArgsCaller.run(ZygoteInit.java:888)
                            at com.android.internal.os.ZygoteInit.main(ZygoteInit.java:749)
                         Caused by: java.lang.ClassNotFoundException: Didn't find class "org.slf4j.LoggerFactory" on path: D


texts:
Currently doing a final project for school on ARCMAP 10.

I have a set of unsuitable and suitable buffers over a big area which is Quebec City, Canada. The final potential sites give me areas between 900 sq meters and over 33,000,000 sq meters. We want to keep areas between 500,000 and 1,000,000. The problem is that some polygons become too thin in certain areas. So, although a polygon is big enough in area, in some spots it is too narrow for all the infrastructure we need.

Any help would be much appreciated.
Thanks

(Original at http://postimage.org/image/ftwsumkt1.)


texts: Removing filepath from ArcGIS geoprocessing result?

texts:
This is my python script code in ArcGIS 10.3;

projectResult = arcpy.Project_management(fc, outFolder + "\\" + fc, template)
#Name of shapefiles re-projected
arcpy.AddMessage(projectResult)

Which gives me this in the geoprocessing dialog result box;

C:\GIS_Result\abc.shp
C:\GIS_Result\1234.shp
C:\GIS_Result\zyx.shp
C:\GIS_Result\987.shp

Can I remove the filepath to give me just the filename?;

abc.shp
1234.shp
zyx.shp
987.shp


texts:
let me explain my SQL server issue

We have two tablesManholes and Basis, which are related to each other by GlobalID/Guid fields.

In the Manhole table i set coded value domain for owner field (1.Company 2.Unknown) and in Basis table another domain for category field (1.CRM_ID 2. None...).

My issue is to restrict our GIS editors to do following action: when they'll choose CRM_ID from basis category field, editors couldn't be able to choose Unknown value from Manholes owner field.


texts: Get ComboBox selected Layer's path as a string to use in OGR/GDAL method

texts:
I'm making a plugin that converts vector .shp to raster .tif. This is the UI of the plugin. The ComboBox is the custom mMapLayerComboBox widget.

After selecting a loaded shapefile in QGIS, what is the type of the layer object that I'll use in the python code as a shapefile? This is the code:

Also, OSGEO/OGR/GDAL needs a string as an argument (the pathfile of the shapefile to be exact) as input. So ogr.Open(layer) gives an error.
This is the error:

But running ogr.Open(layer) in Jupyter, where layer = 'datasets/input/Agriculture.shp' gives no error and runs the code properly.

So,
2. How do I get the pathfile of the shapefile from the ComboBox? (if not ComboBox, what other widgets can I use)?


texts: Carto Builder map not filtering dynamically

texts:
I have a Carto Builder map with several layers and widgets in it. I assumed that when filtering the widgets or the time series that the points on the map would also update to show only the filtered data. For my animated layer, the dynamic filtering is working but the static point layers are not. I have the dynamic option selected under Behavior within the Widgets section and each layer is sourced from the same dataset. What can I do in order for the map to update when the data is filtered through the widgets?

https://acabrera0819.carto.com/builder/843b66b8-0490-433c-a560-3f2b2fa3b386/embed


texts: Calculating multiple text fields based on other fields in ArcGIS Pro?

texts:
I have an attribute table with two types of text fields (34 of each).  One type attributes an "r" or a "d" for each feature.  These fields represent the voting results in either presidential or gubernatorial elections of each county for each election year from 1948 to 2014.  These fields are called prez48, gov50, prez52, gov54, etc.  The second set of fields have "_lng" appended to them (prez48_lng, gov50_lng, etc.)  These fields have no data in them yet.

I would like to automate the process of calculating this second set of fields to turn "r" (from the first set) into "Republican" and "d" into "Democrat".  Example: if prez48 = "r" then prez48_lng = "Republican".

Can this be done without coding?


texts: Using GeoServer REST Configuration to add a datastore

texts:
I have used the example from the Geoserver REST Docs to add a datastore from at PostGIS database.
I have used the following xml doc to add the new datastore:

<dataStore>
  <name>my_datastore</name>
  <connectionParameters>
    <host>localhost</host>
    <port>5432</port>
    <database>postgis</database>
    <user>username</user>
    <passwd>secret</passwd>
    <schema>public</schema>
    <dbtype>postgis</dbtype>
  </connectionParameters>
</dataStore>

The datastore was created with the connection parameters as I expected. But the datastore.xml created looks different from the datastore that I have created with the GeoServer Web Admin.
The connection parameters I am missing is:

Connection timeout
validate connections
max connections
Loose bbox
Expose primary keys
fetch size
Max open prepared statements
preparedStatements
Estimated extends
min connections

So my question is: What are the above mentioned connection parameters corresponding xml elements?

EDIT: I am using GeoServer 2.1-RC1.


texts: How to improve the performance WPS vec:Clip with srsName in wfs:Query?

texts:
I need to export a shape in .zip file with a no native CRS. With the WPS vec:Clip call, the response time is 2.5 minutes.

<?xml version="1.0" encoding="UTF-8"?><wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>vec:Clip</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>features</ows:Identifier>
      <wps:Reference mimeType="text/xml" xlink:href="http://geoserver/wfs" method="POST">
        <wps:Body>
          <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2">
            <wfs:Query srsName="EPSG:3857" typeName="rv:c0506111_ccs2007splus"/>
          </wfs:GetFeature>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>clip</ows:Identifier>
      <wps:Data>
        <wps:ComplexData mimeType="application/json"><![CDATA[{"type":"MultiPolygon","coordinates":[[[[1342724.8071,5836286.1877],[1342718.4154,5836243.9089],[1342716.2899,5836177.9661],[1342712.3319,5836069.7865],[1342697.1044,5835958.9979],[1342715.8319,5835928.4304],[1342748.7523,5835867.138]]]]}]]></wps:ComplexData>
      </wps:Data>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="application/zip">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

Is there a way to get the same result with greater speed?


texts: Performing multiple layer intersection in QGIS?

texts:
I am trying to find a way with QGIS to get the intersection area of more than 2 layers.

All the tools in QGIS take just the 2 input layers.

For example in the image below, I have 4 layers and I want to retrieve just the black area.

Furthermore, I need to retrieve the count attribute and add it to the intersected area attribute table (in this case 4)

I know there is a tool that does the trick in ArcGIS, but is there a QGIS solution?


texts: java vividsolutions jts wgs-84 distance to meters

texts:
my (java) project uses com.vividsolutions.jts.geom, which provides a distance function. After some googling I finally found that the return value from this function is in central angle degrees, when the CRS is WGS-84 (as it is in my project).

What is an easy to implement and accurate on short distances (no need for long distances) conversion from central angle degrees to meters?

I tried converting by dividing by 360 (from degrees), multiplying by 2 * PI (to radians) and multiplying by 6371000 meters, but my answer is off by quite a lot (176m (QGis ellipsoidal ruler tool) vs 204m (mine)).

I found that GeoTools has some conversion functions included, but I would like to not have to depend on another library, since this code is meant to run on the server and an Android client and is meant to be easy to use by other devs.

EDIT: I also emailed the JTS mailing list, and established that wrapping the distance functions for every possible Geometry.distance(geometry) would be a little tedious.

I decided to use the 111 111 m per degree hack described here: Algorithm for offsetting a latitude/longitude by some amount of meters

Thanks for the help, though!


texts: Why I obtain different values calculating distance between two point in the space using this "Great-circle distance" formula and Google Earth tool?

texts:
I am not so into database and GIS and I have the following doubt about this function used on a MySql database to calculate the distance between 2 points.

I am started from this tutorial where there is this function used to calculate the distance between 2 point field:

CREATE FUNCTION earth_circle_distance(point1 point, point2 point) RETURNS double
    DETERMINISTIC
begin
  declare lon1, lon2 double;
  declare lat1, lat2 double;
  declare td double;
  declare d_lat double;
  declare d_lon double;
  declare a, c, R double;

  set lon1 = X(GeomFromText(AsText(point1)));
  set lon2 = X(GeomFromText(AsText(point2)));
  set lat1 = Y(GeomFromText(AsText(point1)));
  set lat2 = Y(GeomFromText(AsText(point2)));

  set d_lat = radians(lat2 - lat1);
  set d_lon = radians(lon2 - lon1);

  set lat1 = radians(lat1);
  set lat2 = radians(lat2);

  set R = 6372.8; -- in kilometers

  set a = sin(d_lat / 2.0) * sin(d_lat / 2.0) + sin(d_lon / 2.0) * sin(d_lon / 2.0) * cos(lat1) * cos(lat2);
  set c = 2 * asin(sqrt(a));

  return R * c;
end

I think that it should return a value in km as unit measurment.

Then I have used it to calculate the distance between 2 point in the space (taken from Google Maps for example Rome and Milan cities, here in Italy).

I am obtaining reasonable values. My doubt is that, then I have check the distance using the ruler tool in Gooogle Earth and I am obtaining different values.

For example considering 2 points that are near (2 km of distance) I obtain something like 300 m of difference between the difference calculated by Google Earth and my previous function.

Considering 2 points that are very distant (for Example Rome and Milan) I obtain something like 100 km of difference between the value calculated by Google Earth and the one calculated by the previous function.

The previous function always give the bigger value.

So I have to use it to calculate on small distance (maximum 30 km) so the error should be little and could be ok for my pourpose.

But why I obtain these difference between my funcion and the distance calculated using Google Earth?


texts:
I'm developing a GP Tool with .Net ArcObjects SDK.

I would like to debug my GP tool but I haven't figured how.

I have tried to enable background processing in Geoprocessing Options and attach process to my Visual Studio but the GP tools crash with the ERROR 000820: The parameters need repair.

If I disable the background processing I can't debug anymore.

How do I debug GP tools developed in .Net ArcObjects?


texts: Calculating attribute with formula in QGIS

texts:
I want to build a QGIS formula where I calculate one attribute from 10 others, and I only need to sum them up.

My Problem is, that I don't get the actual values of the attributes. I thought about using the function: current_value() which works fine if I build a value relation but is not present in the default value expression builder.

Here my example:

The "RMR-Wert" should be calculated of the upper values (the values in the parentheses are saved to the attribute table and database).

Is there another option than depending on the function: current_value()?

The Entrys only look like strings in the formula, in the attribute Table they appear as int. This is because I use "Wertabbildung" (Value references (?)):

So the Columns I want to calculate with are all int:

For the Formula I use the Drag and Drop option and it looks like this:

Until now none of the proposed approaches worked for me.

It could be, because the layer i want to alter is a view from a postgres database. Also in the QGIS attribute table the values are displayed as strings:

while in the database they are stored as integers:

Also i can't use a virtual column, because then my QGIS crashes, i guess it's also because its a database view, but not sure about this.

If i want to change my values to int values i get an error message:

evaluation error: not possible to convert "..." to integer (hope thats a propper translation...)

That could be, because in the attribute table the values appear as strings...

My plan now is to move on from QGIS to Postgresql and try to tackle the problem there, using Generated Columns.


texts:
I'm using leaflet with OSM for two different web mapping projects and have the same problem with both. On both, I set bounds and zoom limits and create the map, and they work nicely on my laptop (Chrome 56.0.2924.87 (64-bit) on Windows 8.1). However, there are problems on mobile: the scaling appears to be wrong, as if it were zoomed to 150%, so that text in the map and images both in the map and elsewhere are blurry, and leaflet controls positioned in bottom-right or bottom-left appear slightly off the edge of the screen. However, this doesn't change if I manually change the scaling in the browser on mobile, and isn't affected by having the line:

<meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=no" />

in my code or

body{zoom:100% !important;}

in my css. This happens exactly the same on Chrome on Android on a couple of different phones and on Safari on iOS.

You can see one of the two sites here (still very beta in lots of ways): https://cymraeg.ling.cam.ac.uk/?map=6

Any ideas why this might be happening?


texts: Google Earth Engine Reprojection to Interrupted Goode Homolosine

texts:
There was a great answer to Google Earth Engine Reprojection to non-EPSG defined CRS about the support of Interrupted Goode  Homolosine projection in GEE.
It seems that since release 22.3 the Homolosine projection is now available in GeoTools. However when I try in Google Earth Engine I get the error Projection: The CRS of a map projection could not be parsed. No problems with Mollweide.
The code I am running is here.
// Get a sample image (Landsat 8) which uses a UTM projection.
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_044034_20170614');

// Get the bounding coordinates of the Landsat image.
var bounds = image.geometry().bounds().getInfo()['coordinates'];

print('Original image (UTM) URL:',
  image.getThumbURL({
    bands: ['B4', 'B3', 'B2'],
    min: 0,
    max: 0.3,
    dimensions: '400'
  })
);

// Project the image to Mollweide.
var wkt = ' \
  PROJCS["World_Mollweide", \
    GEOGCS["GCS_WGS_1984", \
      DATUM["WGS_1984", \
        SPHEROID["WGS_1984",6378137,298.257223563]], \
      PRIMEM["Greenwich",0], \
      UNIT["Degree",0.017453292519943295]], \
    PROJECTION["Mollweide"], \
    PARAMETER["False_Easting",0], \
    PARAMETER["False_Northing",0], \
    PARAMETER["Central_Meridian",0], \
    UNIT["Meter",1], \
    AUTHORITY["EPSG","54009"]]';

var igh = ' \
  PROJCS["Homolosine", \
    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["Interrupted_Goode_Homolosine"], \
    UNIT["Meter",1]]';

var proj_mollweide = ee.Projection(wkt);
var image_mollweide = image.reproject({
  crs: proj_mollweide,
  scale: 300
});

var proj_igh = ee.Projection(igh);
var image_igh = image.reproject({
  crs: proj_igh,
  scale: 300
});

Is the GeoTool version used in GEE less recent than the one integrating this projection?


texts: arcpy.mapping.ListLayoutElements() takes longer in tool than in ArcMap?

texts:
Has anyone else noticed that accessing elements from arcpy.mapping.ListLayoutElements() through a tool script takes longer than within ArcMap? Whether using a specific name, wildcard or nothing, the call takes on the order of 90 seconds or longer to return:

arcpy.AddMessage("  Start: " + str(datetime.now()))
mapDoc = arcpy.mapping.MapDocument(mapDocumentPath)
elements = arcpy.mapping.ListLayoutElements(mapDoc, "GRAPHIC_ELEMENT")
arcpy.AddMessage("  Finished: " + str(datetime.now()))

Start: 2011-05-24 12:29:28.479000
Finished: 2011-05-24 12:31:10.510000

When the call is made within ArcMap the return is instantaneous.

mapDoc = arcpy.mapping.MapDocument("CURRENT")
elements = arcpy.mapping.ListLayoutElements(mapDoc, "GRAPHIC_ELEMENT")

In both instances, I'm referencing the same map document.

Using ArcGIS 10 SP2.

It may not be obvious but the results are instantaneous when the script is run against "CURRENT". Specify any map document other than "CURRENT" and the script will slow down.
I executed a loop as Jason suggested and found little difference between the first call to ListLayoutElements() and the last. The first selection took 28 seconds while the other nine averaged 24 seconds.
The map document has nine data frames and nine corresponding legend elements that were converted to graphics. When these "legends" are present in the document then ListLayoutElements() is slow; when they're deleted then ListLayoutElements() returns immediately.


texts: Limiting Downtime on Map Services Shared Directly to ArcGIS Portal While Doing Schema Changes

texts:
I have a number of feature layers shared directly to Portal for ArcGIS.

If I need to update the schema of these feature classes, can I do so while the feature service is still in use, or must I somehow stop it?


texts: FME Surface Modeller "minimum number of points"

texts:
IN FME workbench 2016.1 I want to create a DEM from Multipatch-geometry contained in a geodatabase, for doing that I thought using the SurfaceModeller could work, where in SurfaceModeller I group the points by the UUID generated attribute:

But the errors are

The specified cell size of '2' units in the y axis is invalid. Cell size must be less than raster extents

Surface Model cannot be constructed. Either all points are linear or initial set of points did not consist of 3 non NAN elevations

For the first error I think it's due to objects with 3 points in less than a diameter of 2 meters
Second error is due to objects with less than 3 points
Questions:
1) how could I exclude/drop/ignore geometries with less than 3 points ?
2)how could I still add a cell value to objects less than 2 m ?


texts: How to Enabling the option to convert layer transparency to color transparency for ArcGIS Server?

texts: i have a trouble for publishing map service with transparent color, but when i make analyze map document, there error "warnings",told me (Enabling the option to convert layer transparency to color transparency for ArcGIS Server).


texts: Convert Linestrings from OSMnx into Multipolygon

texts:
I'm retrieving some geometries from OpenStreetMap with OSMnx library with the following code:
G = ox.graph_from_place('Casco Viejo, Bilbao, Spain', network_type='walk',
                        retain_all=True, buffer_dist = 50, which_result = 2,
                        infrastructure = 'relation["highway" = "pedestrian"]')

which yields the following graph composed by shapely linestrings:

Then I convert the graph into a geopandas geodataframe:
ped = ox.graph_to_gdfs(G, nodes = False)

I've tried this to convert Linestrings to Points and then Points to Multipolygon
Is there a way to convert this linestrings into a shapely Multipolygon:
from shapely import geometry, ops

# combine them into a multi-linestring
multi_line = geometry.MultiLineString(list(ped['geometry']))

merged_line = ops.linemerge(multi_line)

from shapely.geometry import Point, MultiPoint

points = []
for i in range(0, len(merged_line)):
    points.append((Point(list(merged_line[i].coords[1]))))

coords = [p.coords[:][0] for p in points]

poly = Polygon(coords)

This yields a weird wrong geometry:
shape(poly)

If I try:
MultiPolygon(points)

It gives this Error Message: TypeError: 'Point' object is not subscriptable
Is there a way to transform Linestrings into Multipolygon and this into GeoDataFrame?


texts: Cracks between ortho images when floating on 3D surface - ArcScene

texts:
I am trying to create a seamless 3D visualisation in ArcScene 10.1. I have a TIN surface created from 1m contour lines. I have 4 tiles of ortho imagery I would like to drape over the surface. So far I have been using the Base Heights tab in the properties dialog box on the images to "float on a custom surface" so that the imagery takes on the 3D shape of the TIN surface.

However, the result is that most of the image looks great draped over the surface except where the images come together, I get a very jagged seam, as if data is missing. Here is a screenshot:

I have tried increasing the layer offset to see if I could raise the images above any noise or elevation irregularities but it doesn't make a difference, I still get the seam gaps.

When displayed flat (not floating on custom surface) the images fit together perfectly with no gaps.

How can I display this imagery in 3D using ArcScene or other software so that it appears seamless and without gaps?


texts: Failing Unions due to Non-Noded Intersections

texts:
I'm trying to perform lots of unions in GeoPandas. I've got around fifty sets of between 7-12 GeoDataFrames of the same place at different times that I'm trying to union into one layer per set. Right now, about 20% of the unions are failing due to non-noded intersection exceptions. I'm not sure why snapping the layers doesn't solve this issue. The union looks like this:

# Snap the add layer to the base
flat = flatten(base)
add.geometry  = [snap(g, flat, 0.0001) for g in add.geometry]
base.geometry = [snap(g, flat, 0.0001) for g in base.geometry]

# Ensure the geometries are still
add  = ensureValid(add)
base = ensureValid(base)

# Union the add layer to the base
base= gpd.overlay(add, base, "union")

Here are the helper functions:

def flatten(gdf):

    # Try performing a unary union
    try:
        flat = unary_union(gdf)

    except:

        # Try performing a unary union with buffers
        try:
            flat = unary_union(gdf.buffer(0.00001).buffer(-0.00001))

        # Break the MultiPolygons down into LineStrings and do the unary union on those
        except:
            lines = []
            for f in mapping(ensureValid(gdf).geometry)["features"]:
                 coords = json.loads(json.dumps(f["geometry"]["coordinates"]))
                 cpairs = []
                 for a in range(len(coords)):
                     for b in range(len(coords[a][0])):
                         cpairs.append(tuple(coords[a][0][b]))
                 lines.append(LineString(cpairs))
            flat = unary_union(lines)

    # Return the result
    return flat

and

def ensureValid(gdf, geom_type=MultiPolygon):

    # See if there are invalid geometries
    if len(gdf[~gdf.is_valid]):

        # Use a zero buffer to fix problematic intersections
        gdf.geometry = gdf.geometry.apply(lambda g: g.buffer(0.0) if not g.is_valid else g)

        # Remove remaining invalid geometries because I have no idea how to fix them
        gdf = gdf.loc[gdf.is_valid & ~gdf.is_empty]

    # Ensure that every geometry has the correct type
    gdf.geometry = gdf.geometry.apply(lambda g: geom_type([g]) if type(g) != geom_type else g)

    return gdf

Relevant imports:

import geopandas as gpd
import json, shapely
from shapely.geometry import mapping, LineString
from shapely.geometry.multipolygon import MultiPolygon
from shapely.ops import snap, unary_union

Traceback:

    base = gpd.overlay(add, base, "union")

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 371, in overlay
    result = _overlay_union(df1, df2)

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 299, in _overlay_union
    dfsym = _overlay_symmetric_diff(df1, df2)

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 273, in _overlay_symmetric_diff
    dfdiff1 = _overlay_difference(df1, df2)

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 260, in _overlay_difference
    [geom] + list(df2.geometry.iloc[neighbours]))

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 259, in <lambda>
    new = reduce(lambda x, y: x.difference(y).buffer(0),

  File "C:\ProgramData\Anaconda3\lib\site-packages\shapely\geometry\base.py", line 616, in difference
    return geom_factory(self.impl['difference'](self, other))

  File "C:\ProgramData\Anaconda3\lib\site-packages\shapely\topology.py", line 70, in __call__
    self._check_topology(err, this, other)

  File "C:\ProgramData\Anaconda3\lib\site-packages\shapely\topology.py", line 39, in _check_topology
    raise err

TopologicalError: This operation could not be performed. Reason: unknown

What is wrong here?


texts:
I'm using Mapbox to create webmapping, I want to get GeoJSON Features. I tried map.queryRenderedFeatures but it returns array empty.
there are my code:

map.on('load', function(e) {
map.addSource("places", {
    "type": "geojson",
    "data": "data/polygon.geojson"
});

map.addLayer({
    "id": "Surfaces",
    "type": "fill",
    "source": "places",
    "sources-layer": "layer",
    "layout": {},
    "paint": {
        "fill-color": "#088",
        "fill-opacity": 0.8
    }
});

var feature = map.queryRenderedFeatures({ layers: ['Surfaces'] });
console.log(feature);
})


texts: How to choose a fixed position to labels in QGIS 2.14.2?

texts: How can I manually position these numbers for a better position in a polygon layer?


texts: Image-to-image registration to link the 2 images (RGB AND NIR)

texts: I have two images (RGB and NIR) in png format from field and i need to link them together somehow to be able to see the overlapping areas, and cut out the wrong ones. How to do this with QGIS 2.4.0?


texts:
I am trying to use " Create Enterprise Geodatabase" option in arcmap 10.1:

I am getting two following errors:

  You must copy the latest ST_GEOMETRY and dependent libraries to the PostgreSQL software location. Refer to the ArcGIS help topics for more details.

and

  Connected RDBMS instance is not setup for Esri spatial type configuration.
  Failed to execute

I have put C:\Program Files (x86)\ArcGIS\Desktop10.1\DatabaseSupport\PostgreSQL\Linux64\st_geometry.so in /var/lib/pgsql/9.2 directory.


texts: Project points to a vector layer and get z elevation at QGIS

texts:
I have one point file and one dtm vector file (with elevation contours as well as point elevation information).

How can I project points to the vector file and get z elevation of them?


texts:
I would like to know the architecture of QGIS: relationships between QGIS and other open sources, and what elements are implemented in the QGIS.

Usually, There is a diagram showing the architecture of a software such as arcobjects. Even it is simple like Geotools's architecture, it is useful to understand internal structure of a software.


texts:
I recently started having a problem with the google street view addin in ArcMap. When I press on the "Street view" icon this shuts ArcMap down entirely! it used to work fine and I have no idea what is causing this problem.

Anyone face this issue before?

I also get the following error :

  no gui components found in this add-in. add-in version does not match

when attempting to install the add-in from the file


texts: How to get replaceDataSource to work with personal geodatabase?

texts:
I am having trouble using the arcpy.mapping's replaceDataSource with a personal geodatabase. The script runs, but does nothing if I run the script in the python window of ArcMap. If I run it in PyScripter, PyScripter stops working and gives the message: EOFError: [Errno 10054] An existing connection was forcibly closed by the remote host.

It ought to be a simple fix; the feature classes are the exact same, just from a personal geodatabase that has been copied to a different location.

import arcpy

fcSource = r'c:\users\athom\Documents\GIS\GIS_computer_06_20_2013\GIS_LOCAL_MASTER\DATA\Vector_Data\Toponomy\TH_toponomy_250k.mdb'
fcName = 'sl_anno250k'
mapDoc = r'c:\users\athom\Documents\GIS\GIS_computer_06_20_2013\GIS_LOCAL_MASTER\MAPS\Map_Projects\TEST\AAAJ2011-004-02 TH Traditional Territory.mxd'
MxD = arcpy.mapping.MapDocument (mapDoc)
for lyr in arcpy.mapping.ListBrokenDataSources(MxD):
    if lyr.supports ("DATASOURCE"):
        if lyr.dataSource == r"N:\DATA\Vector_Data\Toponomy\TH_toponomy_250k.mdb\sl_anno250k":
            lyr.replaceDataSource(fcSource, 'ACCESS_WORKSPACE', fcName)
            arcpy.RefreshTOC()


texts: Adding another spatial DataBase in PostGIS?

texts:
Is there a way to add another Spatial DataBase in PostGIS in the same SERVER?

Becuse it gives me this ERROR, when i´m trying to do ...

I disconnected the DB postgis_23_Sample and create my new DB using advice from @Vesanto.


texts: Esri js api 4.6 fires animation event on fast drag

texts:
We are having a problem with fast drag event on the new esri js 4.6 api. When a user drags the map too fast, the map extent's changes rapidly.
I've detected that on a fast drag, there's an animation event, which is not documented.
The goTo() function also fires an animation, but it accepts duration, speed and easing as parameters.
A hacky solution is submitted as an answer, I'm still looking for a straight api approach.


texts:
I am not sure if someone already asked this question, but I am having trouble to understanding MapInfo and I only recently started working with ArcMap.

I have done a map (different shp files) in a coordinate system MGI Balkan 6 (the coordinate system I need), but when the client opens this map (shp.files) in MAPInfo the coordinate system they see is not in the "fitting" coordinate system.

How can I solve this problem in ArcMap so when I send it to them they could see it appropriately?


texts: How to export a raster into ascii without nodata values?

texts:
I need to export a grd raster file into ascii .xyz without nodata values. This is not a problem in itself...I´m using Saga-GIS export grid to xyz module.(ArcGIS 10.01 produces a very strange table of xyz values with lots of fields??)...but the person who wants to use the data would like it without the no-data values. The project area is a state border - so there are always going to be no-data values within the bounding box due to the unregular nature of the border.

The data will be used in WINDPRO (which I do not have so I can´t test it).

is there a way to export a raster without getting no-data values?
Does anyone know windpro, and what could the problem be with no-data values?

thanks for any suggestions,

rob


texts: QField and Inputapp with PostgreSQL available on the internet

texts:
I have been using both tools in field work, making use of forms to facilitate filling in information using spatial databases in geopackage format. It has worked very well, but we have to take care to download it to the desktop and incorporate the collected data into the PostgreSQL database on the intranet.
I would like to know if there is any way to do the same field operation, but changing the geopackage for PostgreSQL exposed on the internet, but with some kind of protection to the data and the database server, protecting them from possible hacker invasions.
The technical support team says there is no way to do this without exposing the systems to intranet network security.
Is there any solution to do this?


texts: Using QGIS or GRASS API in custom application

texts:
I've already read in this link that is not possibile to use QGIS or other GIS application API (GRASS, ArcGIS) without rebuild all libraries into my custom code. First I would like to know if that corresponds to the truth or if there is a way to integrate these APIs without the whole application. What I would like to do is to convert a DXF file into a GeoJSON and use it for my application. A possibile workaround is install the server version of QGIS somewhere (I'm considering to add it to my application's docker container) and create a Python application that can leverage QGIS API previously installed.

Is there any other library for this type of conversion?


texts:
I'm hoping to accurately align hundreds of aerial images that are slightly off-kilter.  The images cover the same area on different dates (several acres).  They contain identifiable features around the edge (roads) but vegetation in the middle.

Is there any way to do this?

I know ArcGIS has the automatic georeference toolbar - which I've used with limited success.

Anything else out there?


texts:
Looking to connect to MySQL DB from ArcCatalog.

As of now I have ArcCatalog 10.5 running on machine 1.
I have a MySQL 8 instance running on machine 2 (server).

I am able to connect to machine 2 from machine 1 via HeidiSQL.  I can modify tables, edit users, etc.

I have installed the MySQL OCDB Drivers, done the Microsoft OLE DB Provider for ODBC Drivers in Arc Catalog on machine 1.  However, I am unable to connect to the MySQL instance in any fashion.
I am greeted with this error: Test connection failed because of an error in initializing provider. [Microsoft][ODBC Driver Manager] Invalid connection string attribute.

This error leads me to believe that I might have typed something wrong.  However, I'd like to not use the connection string at all; ideally, it would show me the MySQL instance in the Data Source Name drop down box.

When I load ArcCatalog on machine 2 I see my DBs in the drop down box.


texts: Limit the number of rows allowed in a table in PostgreSQL

texts:
I would like to set a maximum number of rows on a PostgreSQL table. I want users to be able to delete and recreate rows if need be; but there should never be more than x rows in my table. Is there a way to add such a constraint to a table - or a workaround to achieve the same goal?

The purpose of this is to prevent people from adding too many new polygons in a table. I'm looking for either a general SQL solution or a PostGIS-related solution (based on geometry for example).


texts: Stop QGIS opening multiple taskbar items?

texts:
I have QGIS pinned to my taskbar, right clicking shows that it's named "QGIS Desktop 3.8.2". When I run it and open projects some of them open a new taskbar item called "QGIS 3". This also happens when I open an attribute table or layout view.

When I open their properties they are the same except for the "Start in:" field:

They also show different Recent projects when I right click on them.

When I installed the latest version of QGIS I made sure to do a clean install and remove all previous versions. When I go to Help > About they are all running 3.8.2.

How do I stop duplicate taskbar items from opening?

EDIT: I am using Windows 10 but this issue has been around since before I updated from Windows 7


texts: Why is my aspect not displaying the full data?

texts: I am trying to get the slope aspect for the altitude shown above. Can someone please advise me on why this may not be working properly?


texts: Transferring Data From One Table to Another Using Multiple Keys in ArcPy?

texts:
I have data in one field and I want to transfer it to a specific field in another table based on specific keys Such as ID, road_type and veg_type.
My original table looks like this

ID  Road_Type    veg_type      DENSITY
1   1            A             5.26
2   2            B             2.12
3   2            A             .564
3   1            C             1.23
4   3            B             .5
5   2            B             1.5

The Transfer table will have each possible combination of Road_Type and veg_type as a field and I want the density value to be placed under the field based on the ID, Road and veg keys.

ID  R1_VA  R2_VA  R3_VA  R1_VB  R2_VB  R3_VB R1_VC R2_VC R3_VC
1   5.26   0      0      0      0      0     0     0     0
2   0      0      0      0      2.12   0     0     0     0
3   0      .564   0      0      0      0     1.23  0     0
4   0      0      0      0      0      .5    0     0     0
5   0      0      0      0      1.5    0     0     0     0

My code just gives me zeros

Here's what it looks like:

for fc in arcpy.ListFeatureClasses('*M'):

table1 = fc + 'table'
flds = ['ID', 'Road_Type', 'veg_type', 'DENSITY']

search_fc = {f[0:2]:f[3] for f in arcpy.da.SearchCursor(table1, flds)}

uc_flds = ['ID', 'R1_VA', 'R2_VA', 'R3_VA', 'R1_VB', 'R2_VB', 'R3_VB',
         'R1_VC', 'R2_VC', 'R3_VC']

with arcpy.da.UpdateCursor(fc, uc_flds) as Cursor:
    for row in Cursor:
        ID = row[0]

        if search_fc.get(row[0]) == ID and search_fc[1] == '1' and search_fc[2] == 'A':
            row[1] = search_fc[row[0]]

        elif search_fc.get(row[0]) == ID and search_fc[1] == '2' and search_fc[2] == 'A':
            row[2] = search_fc[row[0]]

        #continue with elif statments for each field

        Cursor.updateRow(row)

        print str(row) +str(fc)

Is there anyway to fix it or is there a more efficient solution?


texts: QGIS Gazetteer Plugin connection to internal database

texts: A few posts on here mention that the Astun QGIS-Gazetteer-Plugin could be adapted to reference a PostGIS database. Has anyone done this adaptation and could share it? I want to use a gazetteer stored in a PostGIS database, but don't want to have to create a webserver.


texts: Error generating pdf from QGIS web client under Windows?

texts:
If I push the print button on the toolbar, I get the menu which works fine. However, when I try to print the map, the pdf doesn't generate.
That's why I look into firebug and I get following comment:

"NetworkError: 500 Internal Server Error - ...localhost/cgi-bin/qgis_mapserv.fcgi?map=C:/OSGeo4W/apache/htdocs/qgis-web-client-master/projects/demo.qgs&SERVICE=WMS&VERSION=1.3&REQUEST=GetPrint&FORMAT=pdf&EXCEPTIONS=application/vnd.ogc.se_inimage&TRANSPARENT=true&SRS=EPSG:3857&DPI=300&TEMPLATE=test&map0:extent=1511664.8,6839757.4,1525015.2,6846582.6&map0:rotation=0&map0:scale=75000&map0:grid_interval_x=5000&map0:grid_interval_y=5000&LAYERS=Hausnummer%2Cinfrastruktur%2CGeb%C3%A4ude&OPACITIES=255%2C255%2C255"

Apache access log give this response:

"POST /cgi-bin/qgis_mapserv.fcgi?map=C:/OSGeo4W/apache/htdocs/qgis-web-client-master/projects/demo.qgs& HTTP/1.1" 200 798
"GET /cgi-bin/qgis_mapserv.fcgi?map=C:/OSGeo4W/apache/htdocs/qgis-web-client-master/projects/demo.qgs&SERVICE=WMS&VERSION=1.3&REQUEST=GetPrint&FORMAT=pdf&EXCEPTIONS=application/vnd.ogc.se_inimage&TRANSPARENT=true&SRS=EPSG:3857&DPI=300&TEMPLATE=test&map0:extent=1509449.3,6838361.9,1527249.8,6847462.1&map0:rotation=0&map0:scale=100000&map0:grid_interval_x=5000&map0:grid_interval_y=5000&LAYERS=Hotel%20-%20Ferienhaus%2CCamping&OPACITIES=255%2C255 HTTP/1.1" 500 538

If I manually change the FORMAT=pdf to FORMAT=png it runs correctly.
Does somebody now what I have to change so I can export a pdf?
I'm using Windows.


texts: Is there any navigation dijit or map option , used for panning the map?

texts: I want to add an add-on to my map that will allow the user pan the map by clicking on it. I am thinking something like the zoom toolbar where the user csn click on the +- buttons to zoom in/out. Any ideas?


texts: get WKT for OSM node, Way and relation for a specific tag

texts:
How can I obtain semantically correct WKT geometries for a OSM way or relation or node?

Well node is simple and way as well, but I am struggling with the relations which can contain many ways:

For example the relation https://www.openstreetmap.org/relation/2764083 contains two ways https://www.openstreetmap.org/way/26128247 and https://www.openstreetmap.org/way/205746186#map=18/48.21971/16.35130 which form the border.

However, how can I ensure that when I create a line string the whole (single) border is returned in a semantic correct way? A simple GROUP BY then collect_list for the ID of the relation on any of the nodes matching can only return the separate relations as LINESTRINGs. How could these be merged in the right way or better a proper single geometry be returned in the first place?

I looked at the SQL from https://github.com/mojodna/osm2orc for inspiration but still struggle with the relations.


texts: WMS for Active Fire Data from NASA site - Can't connect

texts:
I am trying to connect to the WMS service of the Active Fire Data that NASA provides. This is the site: https://firms.modaps.eosdis.nasa.gov/web-services/#firms-wms where I get the links from. My procedure is: copy one of the example links

  http://firms.modaps.eosdis.nasa.gov/wms-t/viirs/?TIME=2016-02-01&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=VIIRS_Hotspots&width=512&height=512&BBOX=-20037508.3427892480,-20037508.3427892480,20037508.3427892480,20037508.3427892480&SRS=EPSG:3857

and paste it into the WMS connection in QGIS. It just takes too long, and never actually load anything.I don't know if I just have to leave it the whole night loading, or there is something going on (something wrong).


texts: R - Find "n" closest points to each point in SpatialPointsDataFrame

texts:
I am working with a dataset that is similar to this one:

library(raster)
library(sf)

# Define CRSs
crs.latlon <- "+init=epsg:4326"
crs.sirgas <- "+init=epsg:5880"

# define edges in somewhere within brazil
(spbb <- matrix(c(-53,-44,-25,-20), nrow=2, byrow=F))

# create spatial points in lat/lon
spbb <- SpatialPoints(spbb, CRS(crs.latlon))

# convert to sirgas 2000
(spbb <- spTransform(spbb, CRS(crs.sirgas)))

# make a raster extent object rounded to km:
ext <- extent(5200000, 6046000, 7234000, 7756000)

# create raster
r <- raster(ext, ncol=252, nrow=156, crs=crs.sirgas) # ~3 km resolution
r[] <- runif(ncell(r)) * 10

# create station points
pts <- sampleRandom(r, 500, na.rm=TRUE, sp=TRUE, cells=T) #%>%
#plot(st_geometry(sf::st_as_sfc(pts)), col=sf.colors(4), axes = TRUE, cex=2, pch=19)

#check it
plot(r)
plot(pts, add=T, col="red", pch=19)

For each point in pts, I would like to find the n (i.e. 5) closest surrounding ones. The resulting data frame would be something roughly similar to:

id      nearest     distance
pt1       pt3         14607
pt1       pt204       8540
pt1       pt301       6306
pt1       pt455       4956
pt1       pt337       2145

And so on for each point in the dataset.

What is the fastest way to achieve that?


texts: Using Copernicus Global Land Cover Layers: CGLS-LC100 collection 3 in Google Earth Engine

texts:
I have some excperience of using GEE, but I would like to get to know the dataset named Copernicus Global Land Cover Layers: CGLS-LC100 collection 3. I would like to get more familiar about its possibilities and examples.
Is it possible to find any tutorials or training materials to follow based on this dataset or maybe someone has done something using this dataset to show more possibilities of using it? I am familiar with GEE developers page example, but it doesn't say much.
For example, what do I get, when I don't select 'discrete_classification' band? Is it possible to calculate the land cover change during the monitoring period? How? Is it possible to add a legend to the map for better visualisation?
var dataset = ee.Image("COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019")
.select('discrete_classification');

Map.setCenter(-88.6, 26.4, 1);

Map.addLayer(dataset, {}, "Land Cover");


texts: Delete all non-end point vertices in a line ArcGIS Pro

texts: I have some culvert polylines that were collected and I would like them to be straight without having to manually delete the vertices between end points. Is there a tool to do this for the whole shapefile without Python.


texts: Is there a way to snap 2 lines together from 2 shapefiles?

texts:
I am currently having 2 shapefiles of a roadnetwork. One of them has the best network (lines intersecting correctly etc.) and the other has some attribute info I want. I have no similar attributes to join the two shapefiles and the road networks are slightly different from each other. I therefore cannot use the select by location method, because they are either are not intersecting completely and using a search distance will select other road parts as well (to illustrate my idea: I added a new variable in the blue road network and want to add a 'Yes' for every road segment that more or less overlaps the red road network). Is there a solution for this problem? Maybe it is possible to make the two road segment overlap 100% in some way? I work in ArcMap, but qgis is also possible.

I have added a picture to illustrate the slight difference between both networks:


texts: Configuring OTB for Processing on Linux Mint

texts: I'm using QGIS 2.18.3 in Linux Mint 18. I have downloaded Orfeo Toolbox and want to add it to the processing toolbox in QGIS. I have done it under Windows indicating the bin and lib path, but on Mint doesn't work. I read https://wiki.orfeo-toolbox.org/index.php/QGIS_access_to_OTB_applications but it doesn't seem logic that I have to build it from source just to do that. Anyone knows a easy/better way to do it?


texts: Compute Angle of line in conventional bearing using ArcGIS API for JavaScript?

texts:
I am trying to compute angle between two lines (line segment) using ESRI (ArcGIS, JS) in conventional bearing (Eg: N29'E)

I am using Add Toolbar to draw a line segment where point A is a fixed point (already drawn) and when the user is hovering over the map looking for point B, it should display projected angle.

Has anybody tried anything like this before or any idea?


texts:
I work with several MXD'S and i try to update a group layer named "allHa" containing several layers with this code:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.workspace = r"D:\PROJECTS\road_20\gis"
sourceLayer = arcpy.mapping.Layer(r"D:\PROJECTS\road_20\gis\layers\16_12_15\allHa.lyr")
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"D:\PROJECTS\road_20\gis\\" + mxdname)
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    updateLayer = arcpy.mapping.ListLayers(mxd, "allHa", df)[0]
    dfList = arcpy.mapping.ListDataFrames(mxd, "*")
    for df in dfList:
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if lyr.name == u"allHa":
                 arcpy.mapping.UpdateLayer(df, lyr, sourceLayer, True)
                 print 'UpdateLayer'
    mxd.save()
del mxd

but i get an error:

>>>
antiquities__55-30__403.mxd

Traceback (most recent call last):
  File "C:\Program Files\CCleaner\yaron\shonot\software\---gis---  \tools\YARON_SCRIPTS\UpdateLayer 1 df.py", line 16, in <module>
    arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, True)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\utils.py",    line 182, in fn_
    return fn(*args, **kw)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\mapping.py",     line 1897, in UpdateLayer
    tl._update(rl, symbology_only)
    ValueError: LayerObject: Unexpected error
>>>


texts:
I am looking for a GeoJSON file with all the capitals in this format:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Point",
        "coordinates": [
          37.628173828125,
          55.7642131648377
        ]
      }
    }
  ]
}


texts: Removing erroneous symbology bands from export?

texts:
Every time I export my mxd to PDF (or jpeg, EMF, ai, gif, etc) or print to PDF, I am getting a band of symbology across the page (orange lines across picture).

I found this article on how to fix things, but the suggestions are either not feasible or not helping to resolve the problem (https://support.esri.com/en/technical-article/000002772). I can't really move the layer around in the TOC or change it to a solid fill because the map is too complex and it would cover other things up. Opening and exporting the MXD in ArcPro does remove the lines, but it also removes/changes several other cartographic elements to the map (I am using Data driven pages in the mxd with the data frame clipped to the feature of interest. Data Frame clipping is not incorporated into ArcPro 2.1). I am downloading the ArcPro update, but am not proficient at Pro at all yet so am hoping for another solution.

Has anyone come across this problem and how was it resolved?


texts: Writing log message in QGIS's script algorithms GUI tab for PyQGIS 3?

texts:
I want to write some log in the log tab of processing algorithm User Interface like this post but for QGIS 3 (PyQGIS 3).

It seems that "progress.setText('Youpi')" does not work anymore.

What's the good way ?


texts: Get GeoJson feature by ID and zoom on entity

texts:
I have got a leaflet map with Geojson layers :

I use jquery Autocomplete to search an entity by name or id, I have no problem with that.
I need your help because I'm not able to refer to a specific point in my layer  (without a click on the map). Ideally, I would zoom in on it and open a popup.

Here is my code for loading Geojson markers :

var markers = L.markerClusterGroup({showCoverageOnHover: false,disableClusteringAtZoom: 9});
        var points_rand = L.geoJson(points, {
            onEachFeature: function (feature, layer) //functionality on click on feature
                {
                    //ARRAY FOR AUTOCOMPLETE
                    musees_names_array.push(feature.properties.nom_musee, feature.properties.id);
                    //ARRAY KEY VALUE
                    musees_key_value[feature.properties.nom_musee] = feature.properties.id;

                    {
                var title_ = "
"+feature.properties.nom_musee+"
";
                var subtitle = "";
                var ouverture = "";
                if(feature.properties.site_web != null ){
                    subtitle +="Site web : "
                    +""+feature.properties.site_web+"
";
                }
                if(feature.properties.periode_ouverture != null ){
                    subtitle +="Horaires : "
                    +"
"+feature.properties.periode_ouverture+"

";
                }
                    }
                layer.bindPopup(
                    title_ + subtitle + ouverture
                    )
                layer.setIcon(IconPerso);
                }
            });

After a research with Autocomplete I have a function with the id selected but I don't know how zoom on a specific marker :

$(function() {
$("#search_musee").autocomplete(
    {
        source:musees_names_array,
        minLength: 3,
        select: function( event, ui ) {
                var id_musee = musees_key_value[ui.item.value];// ID selected
                //HERE IS MY PROBLEM
                //I would like something like :
                // map.fitBounds(SPECIFIC_MARKER.getBounds());
        }
    })

});

I did some research like trigger event leaflet but without success,


texts: reverse clip (erase) for spatial point data frame (R)

texts:
I am trying to recreate the suggestion from this link but instead of using it on a sPolygondf, to do the same for a sPointdf. All is well, but I get an unexpected output - when I do the difference between my one point feature and the other I get a result of X (59104); when I do the same for the attribute information in order to join these and create a sPointdf for later steps I get a result of Y (59488).

Can anyone explain why this is so?

Code:

    #------load data
cities_nl <- readOGR('/home/iordamo/Documents/GIS_Workload/DutchCities/allcities_dis.shp')
imgs <- readOGR('/home/iordamo/Documents/GIS_Workload/DutchAgriculturalParcels/mapillary-public.geojson')

#------selection
imgs <- spTransform(imgs, '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs')
imgsdf <- data.frame(imgs)

imgs_in <- imgs[cities_nl,]
imgs_indf <- data.frame(imgs_in)

imgs_out <- gDifference(imgs, imgs_in)
imgs_outdf <- setdiff(imgsdf, imgs_indf)

    > nrow(imgs_outdf)
[1] 59488

> imgs_out
class       : SpatialPoints
features    : 59104

So in the case of the spatial operations I am using the gDifference function from the rgeos package and in the other the setdiff from the dplyr package. To my mind these should result in the same number of features in the end.

link to datasets


texts: Raster and vector layer in the same CRS don't match

texts:
While adding multiple raster and vector layers (and setting them all in CRS WGS84) in QGIS 2.2, they don't seem to overlap. Interestingly enough, only raster files overlap between each other, the same for vector layers. My 'on the fly' projection is included.

As a "hint", raster layers coordinates are 6 or 7 digit numbers, while it's 2-digit numbers for vector layers (area in central Italy). Also to mention is that I'm completely new to GIS so every answer is good enough. Thanks


texts:
I have created scatter plot in GEE but for some reason it doesn't  display the R2 and the way the linear regression displayd is also not good:

as you can see. the R2 isn't displayed at all and the linear regression has NDVI inside it.

This is the code I have used:

// Generate chart from sample
var chart = ui.Chart.feature.byFeature(sample, 'SAR', 'NDVI')
    .setChartType('ScatterChart')
    .setOptions({title: "NDVI vs SAR(db)", pointSize: 3,
    legend: {maxLines: 5, position: 'top'},
    series: {
          0: {targetAxisIndex: 0},
          1: {targetAxisIndex: 1}
        },
        vAxes: {
          // Adds titles to each axis.
          0: {title: 'NDVI'},
          1: {title: 'NDVI '}
        },
                hAxes: {
          // Adds titles to each axis.
          0: {title: 'db'},
          1: {title: 'db '}
        },
    trendlines: {
        0: {
          type: 'linear',
          color: 'lightblue',
          lineWidth: 3,
          opacity: 0.7,
          showR2: true,
          visibleInLegend: true
        },

      },
});

As you can see, the showR2 is true, what is the mistake?


texts: Finding distance between points in two different layers along route in ArcMap?

texts:
I need to find the distance between points in separate layers along a stream network. I have source sites for macro-invertebrates and sites where they are lacking. I need to find the distance between the source sites and the sites where they are not through the stream network.

I tried the linear referencing toolbox but the locate features along route tool can only use one layer of points and has not way to differentiate between sources and not source sites.

Any ideas?

I am using ArcMap 10.3?


texts: Qgis2web mouse cursor changing symbology

texts:
Hovering mouse cursor over polygon on QGIS created webpage changing symbology.

I´m have exported webpage from QGIS, thru qgis2web plugin. On settings I have checked "Highlight on hover". When I move cursor over polygon with complicated symbology (2x simple line, simple fill)(fig. 1), polygon is highlighted (fig.2). When I move out of polygon, symbology is changed and simplified, no line symbology.(fig.3). Same problem occur with polylines.

When I don´t use "Highlight on hover" option, symbology doesn´t change.


texts: Looping using Iterate Features in ModelBuilder of ArcGIS Pro

texts:
I have a model that iterates through a dataset (iterate features), each feature has a defined set of steps that gets processed against another features.
ie - take each feature from the dataset, run a euclidian distance, clip out the euclidian distance based on another dataset and so on.
The model contains many "branches" that does the same process with the exception that the clipping feature is different for each of the processes.
The model runs fine "if the clipping" datasets exists.  If the dataset doesn't exist, the model will fail - so I put a "if data exists" in there and set the true portion to be a precondition of the first step.  Again the model works fine, but...if the clipping feature doesn't exist, the model will still run fine but I get warnings at the end of the process.
"The process did not execute because the precondition is false." and then "All the inputs are not current"
Again, the model runs fine for clipping features that are present but I get the warnings.  I'm assuming because the if dataset exists returns false, but false does nothing except continue on with model.
I'd like to be able to get rid of the warning message that appears after the model completes.  If I run the model, that's fine because I know what the warnings are, but if someone else runs the model, they won't know what the warnings are for.
Any way of handling that?


texts: DSAS Can't Recognize Baseline and Shorelines

texts: I have two shorelines and baseline which were created in accordance with the necessary rules such as attributes, coordinate system etc. All data in a personal geodatabase for DSAS analysis. The problem is that neither the shorelines nor the baseline appear in Set Default Parameters. Where am I making the mistake?


texts: Tracking field personnal using android app and showing in ArcGISonline

texts:
Planning to equip field personnel with GPS enabled phones (android) and collect the location and photograph of assets.

Can we use ArcGIS for android in the phone, use this to upload GPS and photos , and the use ArcGIS online to show these locations ?

Do we need ArcGIS for server or simple we can do it using ArcGISOnline and android SDK ?

Adarsh


texts: Z-value using pointFromAngleAndDistance method in arcpy

texts:
I want to access x, y and z values of a new point from a given angle and distance to a known point using arcpy. I have tried to do it like shown below:

import arcpy
sr=arcpy.SpatialReference(4326)
pt1=arcpy.PointGeometry(arcpy.Point(0,0,0),sr,True)
pt2=arcpy.PointGeometry(arcpy.Point(4,5,90),sr,True)
dis=pt1.angleAndDistanceTo(pt2)
pt3=pt1.pointFromAngleAndDistance(dis[0],dis[1])
print pt3.firstPoint

The result I get is:

4 5 NaN NaN

It only returns x and y values.

How can I get the Z- value too?


texts: Are political borders on Google Maps reliable?

texts:
I'm working on a shapefile viewer for Android, and the file I'm using to test shows the surface land management for the state of Arizona.  I zoomed all the way in on the north-west corner of the state and saw this.  (The dotted gray line is Google's state line; the green line is mine.)

It looks like my math is off, or the shapefile is wrong, but then I zoomed in on the north-east corner of the state and saw this:

Google's state lines miss the Four Corners Monument by quite a bit, while mine are almost dead on.  (I also checked in satellite mode; again, my state lines are almost perfectly on the Monument.)  This led me to read the whole backstory of the Four Corners Monument being in the "wrong" place.  But the final Supreme Court ruling was that the original survey markers are the legal border, not the description of where the borders were supposed to be.

Has Google drawn its state lines incorrectly?  Perhaps they used the nominal borders instead of the survey markers?


texts: Creating a new time column in PostGIS / CartoDB

texts: I received a CSV file that has year (just two digits, '70'), month (as 'JAN'), day and time (as '0310') as separate columns. Bringing in to CartoDB / PostGIS they all appear as strings. I'd like to convert them to a single column that is a real 'date' type so I can display with Torque. What's the best way to go about doing that?


texts: Autocad DXF to Qgis and exporting to KML for viewing in Google Maps

texts:
/ Background: I have a DWG trace made in AutoCad using Gauss-Kruger coordinate system ( a CRS used in Slovenia ). I'm using AutoDesk Map 3D to transform the coordinate system to WGS84.PseudoMercator EPSG:3857 which is required by Google to show traces on maps. /

I saved the drawing as DXF file, opened it in QGIS, set both project and  layers as WGS84.PseudoMercator EPSG:3857 CRS, enabled 'on the fly CRS transformation' and exported the layers to a KML file (Google format) and tried opening them in Google Maps.

However, the trace is shown in the ocean near Nigerian coast instead of Slovenia.

I also tried opening shapefiles in qgis in viewing the trace with OpenLayers plugin but it makes no difference.

I'm new to this software and CRS.


texts: How can I find small distance between polylines?

texts:
I`m trying to find lines which have small distance between them and I am using ArcGIS 10 for this.
I had tried to make points at the both ends of the lines and to detect the distance between points and lines with "Near" but there is no result. Points are lying on the lines and I get always 0 as a result.
What is my task: to have only four lines for every square and if the lines are close one to other, they must have a complete overlap.
I had tried with topology, buffers, near (with distance, angle) and etc. and I can not find a way to detect them.
Can someone help me with this ?

And now small description for the pictures:
1 - This is the general situation. Two squares which every one has to be made from four lines.
2 - The smaller has an error. He is made from 5 lines, which is not allowed. The fifth line is starting from point 1 towards the circled area (2).
3 - Zoom in of the circled area and the problem. The lines are not overlapped. The distance between them is about several millimeters.


texts: Removing pseudo-nodes using QGIS Sextante plugin?

texts:
I have a road network (vector) that i'm interested in carrying out some space syntax analysis on.

Predictably my GIS data is in a pretty terrible state. For example, there are a bunch of roads on the network that look like they join but do not. I've managed to fix some of this (on the surface at least) with the GRASS v.clean tool in QGIS.

My big problem though is polylines. In nearly every case I have long stretches of road that are stored as separate polylines. For example, a long trunk road, split into 5 or 6 separate lines even if they have all snapped together. Worse still in some cases these sections have been given different (or slightly different) road names and have no common identifying attribute. If I run the 'topology checker' in QGIS - it identifies these as 'pseudo-nodes'.

I now want to clean these pseudo-nodes and I was really hoping that there would be a (fairly) simple method of simplifying this network to allow me to run the space syntax analysis tools to actually give a meaningful output. Unfortunately the 'topology checker' only shows you where the errors are, I now need a tool to help me clean them. I had though v.clean from GRASS would do the trick but haven't had much luck. I have been looking at the suggested answers here: Find pseudo nodes in free GIS software

There are a lot of reference to the sextante plugin for QGIS - though I have a feeling this has been incorporated into QGIS's tool now as the 'This plugin is deprecated!' according to the repository. If so, I can't find  "Extract endpoints of lines" and "Join adjacent lines". Does anyone know where I can find these functions - or alternatives?

p.s. using QGIS 2.10 on Ubuntu 14.04 - but open to suggestions of other opensource software/solutions: R; Python, GRASS 7, Open Jump etc.


texts: Can I build a basic Data Dictionary for Trimble Terrasync without Pathfinder Office?

texts: I don't have Pathfinder Office.  Can I build a simple Data Dictionary for Terrasync without it?  I seem to remember a friend once telling me that there is a way to build a data dictionary within the Trimble.


texts: Why is Geoserver publishing my WMS layer in WGS84, when I have specified WyLAM?

texts: I am trying to publish a WMS layer from Geoserver, based on a table in a PostGIS dastabase.  The projection of the original source data is "WyLAM" (EPSG:8155).  It seems that no matter which options I select in the "Coordinate Reference Systems" boxes, the WMS gets published in WGS 84 (GCS_WGS_1984).  I have added the WyLAM projection to both Geoserver and PostGIS, since it's not included by default, but this did not seem to change anything.  For what it's worth, the layer displays in the correct location, so it seems that Geoserver does know that the source data are WyLAM.  How can I get the WMS to publish in WyLAM projection?


texts: Getting the number of records from a table in ArcGIS 10.3

texts: How do you output the number of records in a table in ModelBuilder? I can open the table and see how many records are there but I want to tabulate the number of records available from a number of tables.


texts:
im trying to create a popup in Openlayers that would open when someone clicks a marker. And the marker is setup like this: new OpenLayers.Feature.Vector() and not with the ordinary marker layer.

The marker is then placed inside another layer vesselLayer.addFeatures([markerLayer])

which is then added to the map: map.addLayer([vesselLayer])

Now, How can i get a popup to show when someone selects one of the markers?

Let it be noted that each vesselLayer has only 1 marker. But there can be a 100 dynamically added vesselLayers to the map. (Made it that way for easy access trough the LayerSwitcher).

I tried reading the docs:
http://openlayers.org/dev/examples/select-feature-openpopup.html

But i could not figgure it out.

The error that comes out of the console is saying that it cannot read the property "events" of null..
I double checked the object instantiation and it does get instantiated and the property is a part of the name.

Heres the object logged out right after instantiation.