texts: Repeat a QgsLayoutItem on several pages without using atlas

texts:
I am creating a script to export an attribute table in pdf format with PyQGIS 3. The number of features is different regarding the farmer's number of fields and features are sorted by the field id.
Regarding the number of fields, the attribute table may need to be extended to next page(s). I would like to repeat the main informations about the farm on the next page(s). Like a header only starting from the second page.
I first defined as many QgsLayoutItemLabel as there were page (making if statements regarding the number of features) but it makes my script really heavy.
There is no map in the layout so I would like to know if it's possible to do it without using the atlas composer.
The following script allows to put my header on the 2nd page, and I would like to repeat it on the 3rd page without creating a new QgsLayoutItemLabel.
project = QgsProject.instance()
layout = QgsPrintLayout(project)

BD = project.mapLayersByName('Base de données client')[0]
prop_field = BD.fields().indexFromName('Propriétaire')
prop = BD.getFeature(idCode).attribute(prop_field)

header1 = QgsLayoutItemLabel(layout)
header1.setText("Liste des contraintes - " + prop)
header1.setFont(QFont("Ms Shell Dlg 2", 11))
header1.setHAlign(Qt.AlignLeft)
header1.setMinimumSize(QgsLayoutSize(133.376, 8.215, QgsUnitTypes.LayoutMillimeters))
header1.attemptMove(QgsLayoutPoint(6, 3.590),page=1)

EDIT
I tried the following code that allows to duplicate an item, but I don't understand why both items are at the same location, on the 1st page, even if a different page is called:
# Crédit
def generateLayoutItemCredit(layoutn, pagenb):
    credit_text = QgsLayoutItemLabel(layoutn)
    credit_text.setText("Plan réalisé par le Club agroenvironnemental de l'Estrie")
    credit_text.setFont(QFont("Ms Shell Dlg 2", 11))
    credit_text.setHAlign(Qt.AlignLeft)
    credit_text.setMarginX(3)
    credit_text.setMinimumSize(QgsLayoutSize(130, 8, QgsUnitTypes.LayoutMillimeters))
    credit_text.attemptMove(QgsLayoutPoint(6, 268),page=pagenb)
    layoutn.addLayoutItem(credit_text) #add layout item inside the function
    return credit_text

generateLayoutItemCredit(layout, 1)
generateLayoutItemCredit(layout, 2)

EDIT 2
Following the same logic as the answer provided, here is how a) I repeated a label only from the 2nd page and b) I repeated a label with the text changing for each page
a) Header repeated from the 2nd page
def generateLayoutItemHeader(layoutn, pagenb):
    header_text = QgsLayoutItemLabel(layoutn)
    header_text.setText("Liste des contraintes - " + prop)
    header_text.setFont(QFont("Ms Shell Dlg 2", 11))
    header_text.setHAlign(Qt.AlignLeft)
    header_text.setMinimumSize(QgsLayoutSize(120, 8, QgsUnitTypes.LayoutMillimeters))
    header_text.attemptMove(QgsLayoutPoint(9, 3),page=pagenb)
    return header_text

# additionnal number of pages
pages_sup=4
for i in range(0, pages_sup):
    page = QgsLayoutItemPage(layout)
    page.setPageSize('Letter', QgsLayoutItemPage.Portrait)
    layout.pageCollection().addPage(page)

for i in range(1, pages_sup + 1):
    header = generateLayoutItemHeader(layout, i)
    layout.addLayoutItem(header)

b) Updating the page number
def generateLayoutItemPage(layoutn, pagenb, label):
    page_text = QgsLayoutItemLabel(layout)
    page_text.setFont(QFont("Ms Shell Dlg 2", 10))
    page_text.attemptMove(QgsLayoutPoint(193, 3.590),page=pagenb)
    page_text.setHAlign(Qt.AlignRight)
    page_text.setMinimumSize(QgsLayoutSize(15, 8.215, QgsUnitTypes.LayoutMillimeters))
    page_text.setText(label)
    return page_text

currPage = 1
for j in range(1, pages_sup + 1):
    labelp = 'page' + str(currPage) + '/' + str(pages_sup + 1)
    pagel = generateLayoutItemPage(layout, j-1, labelp)
    currPage+=1
    layout.addLayoutItem(pagel)


texts: Geometry Generator Line to Polygon Feature

texts:
I had great luck using Carlos' answer here to make lines from features in a point-shapefile layer to their labels, but when I tried the same technique with a shapefile layer of polygons, it won't work.

When I use the same expression with a polygon layer, for unknown reasons it just draws a straight line and ignores the intermediate point. What could be the reason for this?


texts: Where can I find recent rainforest (change) land cover data?

texts:
I want to try to get a handle on what's becoming of the world's rainforests right now but not seeming to find recent (i.e. around 2010) data.

The stuff here:

http://www.glcf.umd.edu/data/

as far as I can tell goes up to about 2006, but no further.

I am principally interested in the Amazon and Borneo rainforests.


texts: Shapefile misforming when overlaid on map

texts:
I thought I was attempting something pretty simple - adding a shapefile over a terrain map. This is how the image looks when the shapefile is added to Google Earth:

However, when I process the data through R, I get the following image (this is just part of the picture, showing how the shapefile is different:

Some aspects of the shapefile are preserved, but some regions are inverted, and I'm not sure what's causing these problems. This is the code I'm using to process the file:

PPR -> readOGR("Wetland_Files/PPR_Extent/", "PPRboundary")
PPR -> spTransform(PPR, CRS("+proj=longlat +datum=WGS84"))
PPR -> fortify(PPR)
ggmap(myMap) + geom_polygon(aes(x = long, + y = lat, + group = group), + data = PPR, +  color = "white", + fill = "orange", + alpha = 0.5) + labs(x = "Longitude", + y = "Latitude")

Anyone have any suggestions as to what is going on here?


texts:
A couple times a month, I check out the various tips and tricks CW posts on Stack Overflow and always learn something new.
I went back and forth as to whether to pin this question to a specific program of a specific, major version of GIS software.  Given the amount of change to the ux, I think ArcGIS Desktop 10 should get its own Tips and Tricks question and will be a great resource for new ArcGIS users and ArcGIS 10 users upgrading from 9.3.
Hope to see some good tips that aren't already covered in the ESRI ArcGIS Desktop Tips and Shortcuts.
ArcMap

If the File menu is slow to open: open Customize window, navigate to the file menu, and delete the { Sign in } menu item.

Revert to old(classic) snapping environment: see Question on how to use the editor tool and snapping in ArcGIS 10


texts: Selecting surface constraint using ArcGIS Pro ModelBuilder?

texts:
I created a model using ModelBuilder in ArcGIS Pro to make a LAS ground only layer from LAS files. I first created LAS dataset using LAS files and then used make las dataset layer option. I am able to apply surface constraint to the main LAS dataset but the issue is that surface constraint is not applied in final LAS layer.

I have attached screenshot of my model and the issue faced.

The short version of my problem is The model is not able to select surface constraint automatically.


texts: How to display messages in the QGIS log from within a QTimer callback?

texts:
So I'm new to QGIS and Qt, and trying to write a simple toy python plugin. I'm trying to get the plugin to do "stuff" at regularly timed intervals, and to test, I'm just trying to get the script to write to the plugin's log, which I can view in QGIS.

According to this, I should be using the QTimer class to make an asynchronous thread.

This notes how to actually use the timer, and is largely the script I'm basing mine on.

And this notes how to write to the log.

So here's my code:

def runTick(self):
    QgsMessageLog.logMessage("called runtick", "myPlugin")
    timer = QTimer()
    timer.timeout.connect(self.tick)
    QgsMessageLog.logMessage("starting tick", "myPlugin")
    timer.start(1000)

def tick(self):
    QgsMessageLog.logMessage("tick", "superrefresh")

When I run it, I get:

called runtick
starting tick

and then nothing. This is rather annoying, since I'd like to be able to know if my script is actually running what I tell it to run, or if it has just decided that doing what I tell it to is not worth the trouble.


texts: How to adapt a PostGIS query to modify the input table instead of creating a new table?

texts:
I use the following SQL query to create a table that contains dissovled polygons of an existing table:

-- Schema: public
-- Table: t_union
-- Input table: input_table

CREATE TABLE public.t_union AS SELECT
    row_number() over() AS gid,
    sbqry.geom
FROM
(SELECT
    (ST_Dump(ST_Union(geom))).geom::geometry(Polygon, /*SRID*/) AS geom
FROM public.input_table
) AS sbqry;

Input polygon features:

Output polygon features:

Is there a way to modify the polygons of the input table instead of creating a new table?

Edit 1:

Meanwhile, I've created a SQL query to modify the input table:

INSERT INTO public.input_table
SELECT
    row_number() over() + 100000 AS gid,
    sbqry.geom
FROM
(SELECT
    (ST_Dump(ST_Union(geom))).geom::geometry(Polygon, /*SRID*/) AS geom
FROM public.input_table
) AS sbqry;

DELETE FROM public.input_table WHERE gid < 100000;
UPDATE public.input_table SET gid = gid - 100000;

row_number() over() + 100000 AS gid is a clumsy workaround to ignore the existing gid values and I don't know how to delete the input geometries in a simple way.

EDIT 2:

Perhaps this can be solved by using UPDATE and SET:

UPDATE public.input_table SET (geom) =
(SELECT COALESCE
    ((ST_Dump(ST_Union(geom))).geom::geometry(Polygon, /*SRID*/))
FROM public.input_table);

Unfortunately, an error occurs when executing the query:

set-valued function called in context that cannot accept a set


texts: Faster equivalent of arcpy distanceTo method

texts:
Unfortunately arcpy distanceTo method is painfully slow. When speed is crucial, I am still using Avenue to work around it, but it is getting harder to transfer ArcView 3 to newer computers.

As requested by @PolyGeo code snippets tested on polygon made of 3250 vertices.

------------------

Avenue:

theView = av.GetActiveDoc
Catchments=theView.GetActiveThemes.Get(0)
aTab=Catchments.getFtab
shpFld=aTab.findField("Shape")
feat=aTab.ReturnValue(shpFld,0)
outLine=feat.AsPolyLine
pC=feat.ReturnCenter

aTime0=date.now
for each i in (0..10000)
  d=outLine.Distance(pC)
end
msgbox.info ((Date.Now-aTime0).AsSeconds.asstring,"")

-------------------

Python:

import arcpy, time
infc=r'C:\Temp\theme1.shp'
with arcpy.da.SearchCursor(infc,("FID","Shape@")) as cursor:
    for i,feat in cursor:
        break
outline=feat.boundary()
pC=feat.centroid

start = time.time()
for i in range(10000):
    d=outline.distanceTo(pC)

end = time.time()
arcpy.AddMessage(end - start)

Executing: Script3 centres Start Time: Fri Jul 28 19:32:08 2017
Running script Script3...
83.7339999675751
Completed script Script3...

---------------------

Note that last 5 lines are doing exactly the same thing, calculatiog distance from polygon centroid to it's boundary.
As one can see arcpy function is far behind out of date Avenue.

So, does anybody aware of python module that works significantly faster, please


texts: Openlayers 2.12 - add selectFeature to button results in "Geometry Undefined"

texts:
I have modified the Openlayers mobile Vienna WMS example, and am trying to change the labelButton to allow selecting of features to open a popup. I think it is necessary to do this as on a touch screen, it is hard to differentiate between touch to move and touch to select.

However I get an error "Cannot read property 'geometry' of undefined" when I change the labelButton code (Pasted below). I have tried moving the code block to after the layer is added to the map, but still get the same message. Can anyone help resolve the problem?

   var labelButton = new OpenLayers.Control({
        type: OpenLayers.Control.TYPE_TOGGLE,
        displayClass: "labelButton",
        hover: "true",
        autoActivate: 'false',
        eventListeners: {
                  deactivate:
      function onFeatureUnselect(event) {
            var feature = event.feature;
            if(feature.popup) {
                map.removePopup(feature.popup);
                feature.popup.destroy();
                delete feature.popup;
            }
        },
        activate: function onFeatureSelect(event) {
            var feature = event.feature;
            var selectedFeature = feature;
            var popup = new OpenLayers.Popup.FramedCloud("chicken",
                feature.geometry.getBounds().getCenterLonLat(),
                new OpenLayers.Size(150,200),
                "<h2>"+feature.attributes.name + "</h2>" + feature.attributes.description,
                null, true, onPopupClose
            );
            feature.popup = popup;
            map.addPopup(popup);
        }

        }});


texts: How to get a handle on default MousePosition elements in OpenLayers?

texts:
I have an OpenLayers interface with two maps side-by-side.  They are synced so that any pan or zoom action on either map affects both of them.

I also have MousePosition controls on each map so that the user can see the lon/lat coords on mouseover.  I like having a separate MousePosition on each map rather than a single one in an external div.  I was able to sync the MousePosition controls to one another using the code below.

My question is this:  Is there a better way to get a handle on the MousePosition elements?  The getElement approach I'm using feels a bit hacky.  Thanks.

map1.events.register("mousemove", map1, function (e) {
  var position = map1.getLonLatFromViewPortPx(e.xy).transform(WGS84_SM, WGS84);
  var lat = position.lat.toFixed(5);
  var lon = position.lon.toFixed(5);
OpenLayers.Util.getElement("OpenLayers.Control.MousePosition_68").innerHTML =
                           "" + lon + " " + lat + "";
});

map2.events.register("mousemove", map2, function (e) {
  var position = map2.getLonLatFromViewPortPx(e.xy).transform(WGS84_SM, WGS84);
  var lat = position.lat.toFixed(5);
  var lon = position.lon.toFixed(5);
OpenLayers.Util.getElement("OpenLayers.Control.MousePosition_37").innerHTML =
                           "" + lon + " " + lat + "";
});


texts: Determining EPSG code for coordinate transformation using ArcPy?

texts:
I am trying to write a Python script which will identify the EPSG code for the coordinate transformation.

I am familiar with the dataFrame.geographicTransformations function, however this returns a String describing the name of the transformation.  What I really need is the EPSG code of this.

EPSG codes are possible in a SpatialReference Object, but do not seem to translate into the geographicTransformations.


texts: From text data create a FeatureClass with table with python

texts:
I'm searching how the output of one program (text file with some writings, but also with data as numbers, dates, times, coordinates, all separated with commas) can be transformed to a Feature Dataset.

1) .txt -> .csv -> gp.TabletodBase and then it's nearly done..
But get problems writing csv with Excel dialect - all row is writen in one cell. So, I create the delimiter ';' to get the data into cell by cell.

f= open(rawtable, 'r')
f2 = open(outtable, 'wb')
csv.register_dialect('dotcoma', delimiter=';')
writer=csv.writer(f2, dialect='dotcoma')

for line in itertools.islice(f, 6, 7): #reading labels that are in 6'th row
    labels= line.strip('\n').split(",")
    writer.writerow(labels)

But well, then Arcgis does not recognize the table..
I have been trying to cast the data, like integers, floats, strings.. but that did not help..

  2) .txt -> .shp...   ----> THAT WORKS!!!

The shapefile does not have a type of time.. So was skiping it.

EDITED
However tried the Pyshp library. But after creating fields blindly I record the tuple of all data in a row and I'm getting an error: IndexError: tuple index out of range..

I simply write fields from the list of data, depending to which group they go:

for fieldname in labelkeys:
    if fieldname in t_int:
        w.field(fieldname, 'N', 8, 0) #1probl was here as I made '8' vs. 8 ! stpd
    ...
    else:
        w.field (fieldname,'C', 25)

Another problem with coordinates. I was writing floats... and it did not like it.. So I had to write integers.
Plus one cannot pass a tuple to w.record().
Thanks to StackExchange I've learned the benefits of asterix :) and one simply has to write w.record(*mylist).

3) .txt -> feature class...
Here I have problems writing the entire row at once..(like the list of values) into a NewRow. Since field names are made automatically I can not for each make like row.fieldname=somevariable. I have tried while looping to SetValues  row.setValue(labelkeys[unit], data[unit]) and then at the end of a list rows.InsertRow(row) , but it seemed that it sets just one value. Though I tried by hand and got an error in data casting..

The error comes that I read data like string '34', cast it to int(34) and try to write into field with type 'DOUBLE' (since I grouped some data and wanted to make more generic types).

if labelkeys[unit] in t_int:
    if data[unit]:
        row.setValue(labelkeys[unit], int(data[unit]))

So for now.. nothing is really working..

Maybe there is some other way I have not though about?

If not, could anyone help to get the things work?...


texts: Extracting data from a SpatialPointsDataFrame into a SpatialPolygonsDataFrame in case some polygons have no overlapping spatial points

texts:
Data used in this post:

Climate data (10') from the Climatic Research Unit, which I have changed into a SpatialPointsDataFrame https://crudata.uea.ac.uk/cru/data/hrg/tmc/
Shapefile from European NUTS 3 regions, in R a SpatialPolygonsDataFrame https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts

I am trying to get the average climate for each NUTS 3 region, by overlaying the two spatial objects and taking the average of the climate values which fall within each polygon:

NUTS_precipitation <- over(NUTS_map, CRU_pre, fn = mean)

However, I have many missing values because within some of the polygons fall no points (see the example of Belgium below, NAs in red). The solution would probably be to include also the surrounding climate values, but I have no idea which function to use for this.

I hope I am making myself clear. Can anyone help me?


texts: How to define leaders to a label in QGIS 2

texts:
In QGIS 1.8 and 1.9 I used a plugin called Easy Labels to develop leader lines that would be saved as a memory layer. I have not been able to replicate this with QGIS 2.0, I can't even find the Easy Label Plugin. Assuming it is no longer valid in QGIS 2.

Question:
How can I create a leader from a label moved outside a polygon to the centroid of the polygon?


texts: Combining data frame with bnd as spatial polygon data frame?

texts:
I want to use this data set (publicly available for R)

library(R2BayesX)
library(spBayesSurv)
data("LeukSurv")

Observations: 1,043
Variables: 9
$ time     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
$ cens     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
$ xcoord   <dbl> 0.20507166, 0.28555678, 0.17640573, 0.24476295,
$ ycoord   <dbl> 0.49724366, 0.84895259, 0.73649394, 0.21058434,
$ age      <int> 61, 76, 74, 79, 83, 81, 76, 87, 66, 78, 57, 87, 79, 84
$ sex      <dbl> 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1
$ wbc      <dbl> 13.3, 450.0, 154.0, 500.0, 160.0, 30.4, 41.3,
$ tpi      <dbl> -1.96, -3.39, -4.95, -1.40, -2.59, 0.03, 3.95, 1.91
$ district <int> 9, 7, 7, 24, 7, 11, 17, 21, 18, 7, 17, 9

and there is a boundry file for this data set

> nwengland <- read.bnd(system.file("otherdata/nwengland.bnd", package = "spBayesSurv"))
Read 25774 records
Note: map consists of 29 polygons
Note: map consists of 24 regions
Reading map ... finished

I want to combine these two as SpatialPolygonDataFrame so I can use errorsarlm() or lagsarlm() any other spatial R code which needs listw=()

I think we can use district to do this but I couldn't find how. I believe proj4string for north west england is "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"

BTW: Data frame to SpatialPolygonsDataFrame with multiple polygons doesn't work

@spacedman gave me an idea and tried to do this:

> pcens = SpatialPolygonsDataFrame(nwsp,
              aggregate(LeukSurv,
                        list(d=LeukSurv$district),
                        function(v){LeukSurv}))

What I am trying to do is that I want to add each information to comes from each county to spatial polygon and make it spatial polygon data frame. Is it possible?

PS: I couldn't make it working.


texts: How to update an existing GeoJSON with Fiona

texts:
With the code below I can create a new GeoJSON and add data to it:

data =[(LINESTRING (-98.98234 19.38234, -98.990341 19.342345, -98.992449 19.378621, -98.991299 19.378233, -98.992349999 19.377381), 1)]
schema = {'geometry': 'LineString','properties': {'id': 'int', 'name': 'str'}}
with fiona.open(file_name,'w','GeoJSON', schema) as output:
    for id, name in data:
        output.write({'geometry':mapping(line), 'properties':{'id':id, 'name':name}})

How can I write to an existing GeoJSON? Is that possible with the GeoJSON driver?


texts: How to create a PostGIS layer using QGIS API?

texts:
I want to create a new PostGIS layer and save changes made to it.

Presently, I am connecting Python to PostGIS using psycopg2 but I want to know whether there is any QGIS API function available.

Note: I am able to open an existing PostGIS layer, but I don't know how to create a new one.


texts: Shortest path problem from A to B on routing engines

texts:
I am in the process of creating a routing engine. I have encountered the following problem. Lets say the user gives a point A and point B and expects to get the A->B shortest path. I am using simple Dijkstra for now.

Let's say that I can somehow find the (latitudeA, longitudeA) and (latitudeB, longitudeB) coordinates, which are the closest coordinates to the A and B points that the user inputted. From those coordinates I could then also find the nodeA_ID and nodeB_ID on the graph. The problem is that for those nodes it is very likely that the A->B path doesn't exist at all. For example, if node A was only part of a one-way road that went to the opposite direction that the user wanted to go to.

However, a A' and B' must exist, very close to A and B respectively, so that the A'->B' path exists. So the routing engine should try and find that A'->B' path instead.

Also, that A'->B' path might not even be optimal. There could have been a A''->B'' path, where A'' and B'' were only a couple of meters away from A and B. So the routing engine should find the A''->B'' optimal path instead.

How do routing engines handle this situation ?


texts: Can GDAL run XML Composition File for GeoPDF from command line?

texts:
The documentation page for PDF via GDAL includes the description for using an XML Composition File:
https://gdal.org/drivers/raster/pdf.html#creation-of-pdf-file-from-a-xml-composition-file-gdal-3-0
But, the description only gives examples of the XML Composition file being used from C++ or Python code.
Is it possible to use the XML Composition file directly from GDAL commands on the command line?


texts: Is it possible to merge two polygons with a single touching point into one?

texts:
Is my understanding correct, namely that two polygons touching at a common point cannot be merged because the result would be a self-intersecting polygon?
I know there are dozens of Q&A threads on self-intersecting polygons and I’ve looked through them.
This post: Detecting invalid polygons (self-intersecting or touching) seems to indicate merging 2 polygons with a single common point into one would create a self-intersecting polygon, and would thus be invalid.
Likewise, the answer to this post Is it possible to have self intersecting polygons in shapefile? suggests such a geometry could not be stored in a shapefile.
However, this post seems to indicate it is possible to accomplish it: Combining corner-touching polygons into a single feature.

The reason I ask is that this post: Merging polygons that share one vertex (node neighbors) using ArcGIS Desktop? suggests the user is trying to do exactly that and encountering errors.
I’m posting it as a separate question here because my comment on that last post (poorly framed as a subsidiary question) was deemed invalid. Apologies in advance if I have not done the linking correctly in this question or I’ve deviated from protocol by referring to a deleted comment. I’m trying to improve my understanding of valid and invalid geometry.


texts: Route "buffer" algorithm - polygon around a route line

texts:
What is the specification for an algorithm which takes as input an array of coordinates (lat/lon) representing a route on map and a distance argument, outputting the coordinates of a polygon surrounding the route at the given distance?  I am interested in the math required to calculate the polygon points and how things such as line segments which join at a very acute angle are handled.  Related questions which do not produce the answer I am looking for (algorithm specification / equations):

Buffer along a route

https://stackoverflow.com/questions/7685633/length-buffer-around-polyline

How to create a Polygon around a Linestring as a GeoFence in C# .NET


texts: Modelbuilder Parameters - cannot select input files when model dialogue box is open

texts: I created a model in modelbuilder and set the inputs to be parameters. When I open the model and the tool dialogue box shows up, the input files that I wish to select are not shown. These inputs are used for the "Select layer by location". Why is that and how do I fix it? I have all the input data set up in a geodatabase, and I can see it when I navigate to the .gdb in ArcCatalog.


texts: How can I convert a tif or xyz file to HGT?

texts:
Due to system incompatibilities between projects, I have been asked to convert my DEM which is in geotif format to HGT.

I have no experience with HGT and cannot find any export options to this format.

Would anyone now if/how I can do this?

Yours.


texts: OSMOSIS: How to extract all railway stations in use?

texts:
I'm trying to retrieve a list of all railway stations of Belgium that are in use. I downloaded the country osm from http://download.geofabrik.de/europe/belgium-latest.osm.pbf. I tried the following command to extract all stations:

osmosis --read-pbf file=belgium-latest.osm.pbf --tf accept-nodes railway=station,halt --tf reject-nodes disused=*,abandoned=*,railway=disused,abandoned --tf reject-relations --tf reject-ways --write-xml be_railway_stations.osm

However, in my output I still have disused stations, e. g.

<node id="1313036421" version="6" timestamp="2012-07-06T13:22:20Z" uid="445028" user="Sebke" changeset="12130333" lat="50.8773652" lon="4.3551874">
    <tag k="abandoned" v="yes"/>
    <tag k="disused" v="yes"/>
    <tag k="name" v="Laken - Laeken"/>
    <tag k="name:fr" v="Laeken"/>
    <tag k="name:nl" v="Laken"/>
    <tag k="railway" v="station"/>
</node>

Where do I go wrong?


texts:
osm2pgrouting works without giving any errors but most of the roads have lost their geometry.

Result is ;

What I've done so far;

OSM Data downloaded from http://download.geofabrik.de/ in pbf format.
Data clipped and converted to osm format using the command below

osmconvert turkey-latest.pbf -B=ankara.poly -o=ankara.osm

Data is inserted into PostgreSQL using the command below;

osm2pgrouting --f ankara.osm --conf mapconfig.xml --dbname=workshop --host=localhost -p 5432 --username postgres -W 1 --clean


texts: Increase the size of overlapping lines in QGIS 3.6

texts:
In QGIS I have a layer that is a MultiLineString that has many overlapping lines, which represent driving routes around the city which are generated using OSRM.

My data looks something like

[
[[a,b],[c,d],[e,f],[g,h]],
[[y,z],[a,b],[c,d],[e,f],[g,h]],
[[e,f],[g,h],[i,j],[k,l]]
]

where some segments are only recorded once, some twice, and some many more times (up to thousands of times).

How can I increase the width (line thickness) of segments of the routes based on the number of times they are overlapped?

I have a working solution using the JavaScript libraries OpenLayers and Turf but the solution is inefficient and cannot scale. See below for an example of the desired effect.

Answers that involve re-formatting the data or working in Python are within my scope.


texts: Installing PostgreSQL with ArcGIS Desktop without breaking Python?

texts:
I am looking to learn PostgreSQL.

I have a current working version of ArcGIS 10.5 with Python version 2.7. I have heard of people having installed another version of Python and having trouble with ESRI Python working after that.

I want to install PostgreSQL (BigSQL installer). This installer includes a version of Python.

Will installing PostgreSQL cause Python in ArcGIS to stop working?

I am trying to understand what Postgresql installer BigSQl does with my current installed ArcGIS python version. I am hoping that PostgreSQL and ArcGIS python programs run independent of each other--then there is no issue. If there are issues, I'd like to know BEFORE I install PostgreSQL. There is documentation that ArcGIS program not "playing nice" with other Python versions
Installing ArcGIS 10 with Python(x,y)

In Researching, I have discovered that Postgresql Bigsql installer in windows version does change the Python executable in the Pythonpath environment. The source of my info can be found Here. Does the installer “overwrite” ArcGIS 10.5 Pythonpath OR does the Bigsql installer “add” to the Pythonpath and thus both co-exist without issues?


texts: Converting ESRI Raster Dataset to QGIS compatible Dataset?

texts: I am trying to convert an ESRI .gdb "File Geodatabase Raster Dataset" (Format: FGDBR) to an QGIS compatible File with or without using ArcMap.


texts: How do I create an atlas with multiple maps for each coverage layer with QGIS?

texts:
I am doing a NSW transport map that would like to segregate using Local Government Areas (LGAs). However, I specified that each map should be in the same scale, hence I would like to have many "same-scale" submaps per LGA maps.

Is it possible to do that?


texts: mapserver displays raster image in mirror

texts:
I've installed Mapserver and used a map file to set a raster layer and vector(shp) layers on top of it. However, although displayed in the expected order, the png file is displayed in mirror(compared to it while being built in OpenEV):

any ideas on what is the cause of this and a possible remedy are welcomed. Thanks!

from the .map file, the raster layer:

  LAYER
    NAME "level_2"
    TYPE RASTER
    STATUS DEFAULT
    DATA "lv2.png"
  END

gdalinfo on the .png file

Driver: PNG/Portable Network Graphics
Files: /var/www/mapsdemo/e/data/lv2.png
Size is 1591, 967
Coordinate System is `'
Metadata:
  Software=gnome-screenshot
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  967.0)
Upper Right ( 1591.0,    0.0)
Lower Right ( 1591.0,  967.0)
Center      (  795.5,  483.5)
Band 1 Block=1591x1 Type=Byte, ColorInterp=Red
Band 2 Block=1591x1 Type=Byte, ColorInterp=Green
Band 3 Block=1591x1 Type=Byte, ColorInterp=Blue


texts:
The data for this example can be found here:
https://github.com/davidbaniadam/maps

The problem is essentially the same with Alaska and mainland USA, i.e. I want the size of Alaska to be correct in relation to the other states.
So I am creating this plot:

But I want to keep the actual size of the "small island" with the text "Bornholms Hospital". How can I do that? I have also used the tmappackage but had problems with the labels (could not control them). The code is below. At the end I use print with viewport. Is there another way to combine two plots?

library(sf)
library(dplyr)
library(ggplot2)
library(stringr)
library(grid)

geo <-  st_read("ADM/KOMMUNE.shp")

# This is my plot
geo_RH_sj <- geo %>%
  filter(REGIONKODE=="1084")  %>%
  filter(!stringr::str_detect(KOMNAVN,"Christians") & !str_detect(KOMNAVN, "Bornholm"))

geo_RH_bh <-  geo %>%
  filter(KOMNAVN == "Bornholm")

df <- data.frame(place = c("Skt. Lukas Stiftelsen",
                           "Arresødal",
                           "Diakonissestiftelsen",
                           "Bornholms Hospital",
                           "Søndergaard"),
                 Y =   c(55.727444, 55.976848, 55.678323, 55.092019, 55.743852),
                 X =   c(12.555917, 12.031369, 12.515989, 14.707619, 12.323568),
                 Antal = c(24, 12, 16, 2, 16)
)

df1 <- st_as_sf(df, coords = c("X", "Y"),
                crs = 4326,
                remove = FALSE,
                agr = "constant")

mainland <- ggplot(data = geo_RH_sj) +
  geom_sf(fill= "#fecc5c", color = "#fecc5c") +
  geom_sf(data = df1) +
  geom_sf_label(data = df1, aes(label = place), nudge_x = c(2000, 0, 0, 0, 0),
                nudge_y = c(2000, 2000, -2000, 0,  2000)) +
  coord_sf(xlim = c(677299.3, 739128.7), ylim = c(6162392.4, 6225395.8), expand = FALSE) +
  theme_void() +
  theme(panel.grid.major = element_line(color = "white"))

a <- st_bbox(geo_RH_bh)

bornholm <- ggplot(data = geo_RH_bh) +
  geom_sf(fill= "#fecc5c", color = "#fecc5c") +
  geom_sf(data = df1) +
  geom_sf_label(data = df1, aes(label = place), nudge_x = c(0, 0, 0, 15000, 0),
                nudge_y = c(0, 0, 0, 2000,  0)) +
  coord_sf(xlim = a[c(1, 3)], ylim = a[c(2, 4)], expand = FALSE) +
  theme_void() +
  theme(panel.grid.major = element_line(color = "white"))

mainland
print(bornholm, vp = viewport(x = 0.82, y = 0.5, width = 0.3, height = 0.2))


texts: Assign Mask From Layer in SearchCursor Loop In ArcPy: Error Accessing Environment

texts:
I'm trying to run a loop based on a SearchCursor that selects layers from an input polygon and polyline. I want to select a feature from the polygon to use as a mask for running Euclidean Allocation on selected lines within the polygon feature. The goal is to have multiple areas of Euclidean Allocation run on a subset of polylines, each masked by a corresponding polygon.

Here's my original code. Note that this function is part of a much larger, otherwise functional script.

The input parameters for the three functions below are inline = a polyline feature class, inpoly = a polygon feature class, and valfield = a single field which corresponds to both inline and inpoly: i.e., each value of valfield = x in inline corresponds to lines that fall within inpoly features of the same value. All of these parameters are derived from in-memory feature classes fed in by earlier functions in the script--they're not located directly on disk.

def eucal_by_poly(inline, inpoly, valfield):
    polylist = []
    with arcpy.da.SearchCursor(inpoly, ("OID@", "SHAPE@", valfield)) as cursor:
        vallist = sorted({row[2] for row in cursor})
        for v in vallist:
            lyr = "lyr"
            selexp = "{0} = {1}".format(arcpy.AddFieldDelimiters(lyr, valfield), v)
            mask_lyr = arcpy.MakeFeatureLayer_management(inpoly, lyr, selexp)
            arcpy.env.mask = mask_lyr
            line_lyr = arcpy.MakeFeatureLayer_management(inline, lyr, selexp)
            eucalsub = arcpy.sa.EucAllocation(line_lyr, "", "", "", "LAT_ID")
            eucalsubpoly = arcpy.RasterToPolygon_conversion(eucalsub, "eucalsubpoly", "SIMPLIFY", "Value")
            polylist.append(eucalsubpoly)
    polyeucal = arcpy.Merge_management(polylist, "polyeucal")
    del cursor
    return polyeucal

When I run this code I get this error:

Traceback (most recent call last):
  File "G:\tgf_python\tools\py_scripts_working\transect_working_rasversion.py", line 406, in <module>
    polyeucal = eucal_by_poly(cliptrans, eucal_poly, IDfield)
  File "G:\tgf_python\tools\py_scripts_working\transect_working_rasversion.py", line 397, in eucal_by_poly
    arcpy.env.mask = mask_lyr
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 529, in set_
    self[env] = val
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 581, in __setitem__
    ret_ = setattr(self._gp, item, value)
RuntimeError: Object: Error in accessing environment <mask>

I've tried several variations, including using arcpy.env.mask = "SHAPE@", to no avail. I also tried saving the poly layer to a real file on disk, both as a feature class and a rasterized version. In all cases I get the same "Error in accessing environment" message when trying to apply the mask.

Does anyone have a solution?  I'm using ArcGIS 10.1 with an Advanced License and Python 2.7 on Windows 7. (Any other suggestions on improving code are welcome.)


texts: XY event source layer problem (converts into feature class after running python script ) in ArcMap 10.1

texts:
I am using ESRI 10.1 suit (ArcGIS Desktop & ArcGIS Server, No ArcSDE) at my end and after running python script n number of times event layer converts into feature class.

Following steps I am following :

-I am using query layer in ArcMap & fetching records from MS sql server, table type is plane table.

-Created XY event source layer from query layer and displaying points on the map.

-layer is classified using unique value symbology.

-My query layer is updating five times in day and using python script (example2) to automate the unique value symbology (using windows scheduler) five times in day.

-Overwriting mxd on ArcGIS Server 10.1  using ArcPy in specific time interval.

-After 10 times my XY event source is automatically converts to feature class then I am not able to use ArcPy script to update the legend at my end

So is there any bug related to XY event source layer or am I doing something wrong??

Researching on Google and ArcGIS web help not getting any reason/limitation for that.

Any help will be great !! thanks in advance !!

EDIT :

Legend update python script

Map service overwrite script


texts: Generate N polygon centroids based on attribute in QGIS

texts:
I have a number of polygons with an integer value in a field. I would like to create N number of centroid points for each polygon where N corresponds to the value in the polygon's field.

So sort of like Random Points in Polygon, but centroids instead of Random.

I've done it before by generating random points, adding the centroid X and Y as a attributes to the source polygon and then rejoining the points with the source polygon's fields with Join By Location. And then replacing the geometry in the random points with the source polygon's centroid X and Y.

But this feels suuuper clunky.

So is there some easy way to accomplish the same thing? Maybe some succinct SQL for Virtual Layers? Or clever Geometry Generator snippet for the polygon layer (generating new data is not necessary)?


texts:
I would like to load a web service as a background in Trimble Pathfinder Office 5.30.  But I always get an error indicating that either it couldn't understand the address or received data in an unexpected format.  I have not been able to find any reference or example of a working address.  (The help documentation assumes a working address without any specifications about what is acceptable.)
As an example, can anyone tell me the correct address format to use ESRI's World Imagery layer? (https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer)
Or can you tell me an address that does work?  (I would prefer OpenGIS over ArcIMS, which are the only two types I can select.)


texts: Counting points within polygon using PyQGIS?

texts:
I have one polygon and onde shp of points (some within and others without the polygon).
I have tried the following code:

layers = QgsProject.instance().mapLayers()

for feature in layers['poly_teste_1'].getFeatures():
    geom_poly=feature.geometry()

for feature2 in layers['points_test2_2'].getFeatures():
    geom_point=feature2.geometry()
i=0
for points in layers['pointo_test2_2'].getFeatures():
    if geom_poly.contains(geom_point):
        i=i+1

print(i)

The problem is that the code a printing the total number of points, and not only the ones that are within the polygon.


texts:
I am using QGIS 2.18. Upon making updates to my map, the update is not showing in the export to map web page. I have updated it several times and restarted QGIS but it is not replicating my label styling.

I updated my map to the first image using the size assistant function for the proportional circle. However, when I export to map using the QGIS2web plugin, it still gives me the old representation called capture 1.

1: https://i.stack.imgur.com/yYL8o.png This is the result, it is still not updating in qgis2web.


texts: Why is the Google Maps circle shape not matching to the distance calculated with PostgreSQL?

texts:
I am running a PostgreSQL 9.1 Database with PostGIS 1.5 extension. Several coordinates are stored in there with type GEOMETRY (SRID 4326).
I'm building a web-application with Google Maps Javascript API v3 to display those points. For any point on the map and a definable distance, I now want to highlight markers if they are within the radius of that point.

This query returns the IDs of all those points:

WITH chosen AS (SELECT * FROM ST_SetSRID (ST_MakePoint ($1, $2), 4326))
SELECT mypoints.id FROM mypoints, chosen
WHERE ST_DWithin (mypoints.location::geography, chosen.ST_SetSRID::geography, $3);

The ::geography casts are needed here to use meters as measurement unit instead of degrees. $1,$2 is lon,lat, $3 is the chosen distance in meters.

Clientside, only the markers resulting from the query receive full opacity and I also display a google.maps.Circle with the chosen distance as radius. As you can see on the picture below, some markers that should be highlighted if you trust the red circle are not returned by my query. The difference also only seems to occur on the X-axis.

Now what is wrong here?

Is the circle shape from the GMaps JS API inaccurate for my region (Germany)?

Or am I missing something in my query?

So far, i really think that martin_f is right about the distortion. There are multiple coordinate systems involved here: https://stackoverflow.com/questions/15015946/google-map-api-v3-projection

Can anybody confirm that every Shape I draw in Google Maps Javascript API is in SRID 3857 and every Marker I add with Lat/Lon Values in SRID 4326 are internally converted to 3857 before they are shown on the map?


texts: Comparing ArcGIS Desktop and ArcGIS Server?

texts: What is the primary difference between ArcGIS Desktop and ArcGIS Server?


texts: Plotting netCDF4 data with the right projection

texts:
I guess this will sound simple to you but I just cannot get the right plot within R. I have a netcdf4 file and would like to visualize it within R just like I would do in QGIS and ArcGIS.

sample data

I tried to use insights from similar existing posts but did not succeed. The image get squeezed as shown below:

library(raster)
mycrs <- CRS("+proj=lcc +lat_1=35 +lat_2=51 +lat_0=39 +lon_0=14 +k=0.684241 +units=m +datum=WGS84 +no_defs")
pr <- raster(file.nc4, varname="precipitationCal")
projection(pr) <- mycrs
pr

  pr

class       : RasterLayer
dimensions  : 3600, 1800, 6480000  (nrow, ncol, ncell)
resolution  : 0.1, 0.1  (x, y)
extent      : -90, 90, -180, 180  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=35 +lat_2=51 +lat_0=39 +lon_0=14 +k=0.684241 +units=m +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\3B-HHR.MS.MRG.3IMERG.20150405-S170000-E172959.1020.V03D.HDF5.dap.nc4
names       : precipitationCal
zvar        : precipitationCal

plot(pr)

NB: my region of interest is lat=c(41,84) and lon=c(-142,-53)

file.nc4 attributes:

Grid.GridHeader: BinMethod=ARITHMETIC_MEAN;
Registration=CENTER;
LatitudeResolution=0.1;
LongitudeResolution=0.1;
NorthBoundingCoordinate=90;
SouthBoundingCoordinate=-90;
EastBoundingCoordinate=180;
WestBoundingCoordinate=-180;
Origin=SOUTHWEST;

EDIT:  Using @mdsummer's suggestion below:

r <- '3B-HHR.MS.MRG.3IMERG.20150405-S170000-E172959.1020.V03D.HDF5.dap.nc4'
pr <- raster(r, varname="precipitationCal")
library(maptools);
data(wrld_simpl)
projection(wrld_simpl)
proj4string(pr)=proj4string(wrld_simpl)
plot(t(pr))
plot(wrld_simpl, add = TRUE)

I got this image but it is still not on the right latlon projection?.


texts: Finding geometry of cluster from points collection using PostGIS?

texts:
I have a regular grid, each point in the grid belongs to some cluster, and I need to find geometry of this clusters.

My tables look like so (only necessary fields):

grid_nodes(
    *id uuid,
     geom GEOMETRY(Point)
)

clusters(
    *id uuid
)

grid_node_cluster(
    grid_node_id uuid,
    cluster_id uuid
)

I have not experience with Postgis and GIS, so I've been googling for couple of hours, and found a few similar (but not the duplicates) questions:

polygon from exterior points
finding the outermost border of a set of geomertries (circles)
Using postgis ST_ClusterWithin() on a table

where I found some clues, but not found the answer.

At the end I wrote SQL query like this, but unsuccessfully:

SELECT
    ST_ConcaveHull(geom, 0.8, true) as geom,
    st_geometrytype(ST_ConcaveHull(geom, 0.9, true)),
    uuid_generate_v4() as id
FROM (
    SELECT unnest(ST_ClusterWithin(gn.geom, 0.025)) as geom
    FROM ds_forecast_territory.grid_nodes as gn
    JOIN ds_forecast_objects.grid_node_clusters as gnc
        ON gn.id = gnc.grid_node_id
        AND gnc.cluster_id = 'ef2e3314- d0c9-4751-9b0f-cd5b7cf3c4e5') as clusters
WHERE st_geometrytype(ST_ConcaveHull(geom, 0.9, true)) = 'ST_Polygon'

where clause in the query and uuid are just dirty hacks for QGIS. This query somehow works, but resulting geometry is messy and inaccurate.

This is example of the points, and borders which I want to find:

And this is the result I was able to get with my query:

Can I solve my issue using Postgis, or is it impossible and I need to use some spatial libraries for Python?

I've found it's probably called alpha shape, but I still can't solve the issue.

With help from @Spacedman and @dbaston, the answer;

DO
$$
DECLARE
  _cluster_id UUID;
  _cluster_geometry GEOMETRY;
BEGIN
  FOR _cluster_id IN SELECT id FROM ds_forecast_objects.clusters
  LOOP
    SELECT st_union(st_expand(gn.geom, 0.0041, 0.0023)) INTO _cluster_geometry
    FROM grid_nodes as gn
      JOIN grid_node_clusters as gnc
          ON gn.id = gnc.grid_node_id AND gnc.cluster_id = _cluster_id;

    INSERT INTO clusters_geometry(cluster_id, geom)
      VALUES (_cluster_id, _cluster_geometry);
  END LOOP;
END
$$

And the result in QGIS:


texts: GeoTools UnsatisfiedLinkError in display raster layer of RPF data

texts:
Background: Used this tutorial as base

System: Eclipse IDE in Windows 10 VM

So I have some RPF data that I'd like to display within a JFrame. Right now I'm trying to figure out a way of just displaying a single file (.i42) within a layer and then my next object is to display a image mosaic populated by multiple files of the same type.

It looks like I'll need to use RPFTOCFormat class to get a Reader to build the layer. Using Maven to import the library, my pom.xml (relevant section only) is the following:

<properties>
    <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
    <geotools.version>19-SNAPSHOT</geotools.version>
</properties>

<dependencies>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-imageio-ext-gdal</artifactId>
        <version>${geotools.version}</version>
    </dependency>
</dependencies>

Within my code (very similar to the tutorial referenced above) I try and use this Reader to build a GridReaderLayer to populate my map with.

private void getLayersAndDisplay() throws Exception {
    List<Parameter<?>> list = new ArrayList<>();
    list.add(new Parameter<>("image", File.class, "Image", //element id and label before file path box
            "GeoTiff or World+Image to display as basemap", //description below file path box
            new KVP(Parameter.EXT, "tif", Parameter.EXT, "jpg"))); //expected file extensions
    //list.add(new Parameter<>("shape", File.class, "Shapefile", "Shapefile contents to display",
    //      new KVP(Parameter.EXT, "shp")));

    JParameterListWizard wizard = new JParameterListWizard("Image Lab",
            "Fill in the following layers", list);

    int finish = wizard.showModalDialog();

    if(finish != JWizard.FINISH) {
        System.exit(0);
    }

    File imageFile = (File) wizard.getConnectionParameters().get("image");  //get path from associated element id
    //File shapeFile = (File) wizard.getConnectionParameters().get("shape");
    //displayLayers(imageFile, shapeFile); //display those layers from file objs
    displayLayers(imageFile);
}

private void displayLayers(File rasterFile) throws Exception {
    AbstractGridFormat format = new RPFTOCFormat();

    reader = format.getReader(rasterFile);

    final MapContent map = new MapContent();
    map.setTitle("Jonathan's attempt");

    Style rasterStyle = createRGBStyle();
    Layer rasterLayer = new GridReaderLayer(reader, rasterStyle);

    map.addLayer(rasterLayer);

    frame = new JMapFrame(map);
    frame.setSize(800, 600);
    frame.enableStatusBar(true);
    frame.enableToolBar(true);

    JMenuBar menuBar = new JMenuBar();
    frame.setJMenuBar(menuBar);
    JMenu menu = new JMenu("Raster");
    menuBar.add(menu);

    frame.setVisible(true);
}

However when I run the class, I get the following error:

it.geosolutions.imageio.gdalframework.GDALUtilities loadGDAL
WARNING: Failed to load the GDAL native libs. This is not a problem unless you need to use the GDAL plugins: they won't be enabled.
java.lang.UnsatisfiedLinkError: no gdaljni in java.library.path

I thought I imported the lib via that dependency specified within the maven pom but apparently I'm mistaken.


texts: Spliting layers from multispectral images

texts: I have an image with 4 bands and I wanted to De-layer stack to have each bands separately. How and with which software I can do that?


texts: Using QGIS on network with heavy firewalls?

texts:
I'm having a problem keeping QGIS 2.18 running on a government network computer with heavy firewalls.

The software downloaded with no apparent issues, however I do not have access to any web layers, web tools, or plugins. And the software is really unstable.

I used a local raster and pulled a local .txt table and it couldn't handle this without timing out and forcing a shut down.

Could this be an Internet problem even though this app runs from my local drive?


texts: What could be a proper and simple invalid WKT/WKB geometry?

texts:
We are converting a geometry table to another one and we would like to convert every invalid geometry to a single one easily identifiable by users in WKB or WKT form and not breaking the processing or display.

What could it be?

We do not want to use 'POLYGON EMPTY' because it could be something that exists in the original table. It could also be the result of some geoprocessing done on the converted table and then these ones would not be distinguiable from the original invalid ones...

Any better idea?


texts: What is the problem with the polygon? Getting a *"found non-noded intersection"* exception error


texts: Geoserver Dynamic Colormap using CQL won't validate

texts:
I'm trying to create a dynamic heatmap using GeoServer 2.8.3, but when it comes to validating the colormap-part of my SLD i keep getting the same error:

line 66: cvc-datatype-valid.1.2.1: '${env('high',1)}' ist kein gültiger Wert für 'double'.
line 66: cvc-attribute.3: Wert '${env('high',1)}' des Attributs 'quantity' bei Element 'ColorMapEntry' hat keinen gültigen Typ 'double'.

My code looks as followed

<ColorMap type="ramp">
                <ColorMapEntry color="#FFFFFF" quantity="0" opacity="1.0"/>
                <ColorMapEntry color="#4444FF" quantity="0.3" opacity="1.0"/>
                <ColorMapEntry color="#FF0000" quantity="0.7" opacity="1.0"/>
                <ColorMapEntry color="#FFFF00" quantity="${env('high',1)}" opacity="1.0"/>
              </ColorMap>

While the first three entrys work just fine, the fourth is causing the error to appear. I am pretty confused since i took the code right from the documentation (http://docs.geoserver.org/2.8.x/en/user/styling/sld-reference/rastersymbolizer.html#cql-expressions). Since i want the colormap to be fully customizable by the getmap-method, i would like to make every colormap-entry (color aswell as quantity) dynamic using the env-function.

Putting in ('high',1.0) instead of ('high',1) wont work.

I also found another question regarding the same matter, but it didnt help either (How do I use env parameter variables in an SLD colormap?).

Did anyone experience something similar?

The dynamic quantitys do work, but once i apply the same to the colors, it won't display the map at all.

<ColorMap type="ramp">
                <ColorMapEntry color="${env('c1',#FFFFFF)}" quantity="${env('zero',0)}" opacity="1.0"/>
                <ColorMapEntry color="${env('c2',#4444FF)}" quantity="${env('low',0.3)}" opacity="1.0"/>
                <ColorMapEntry color="${env('c3',#FF0000)}" quantity="${env('med',0.7)}" opacity="1.0"/>
                <ColorMapEntry color="${env('c4',#FFFF00)}" quantity="${env('high',1)}" opacity="1.0"/>
              </ColorMap>

Do i have to do something special with the color-codes? I am not getting any error outputs for these.


texts: How to calculate the hill area in Google Earth

texts:
I want to calculate the area of this hill (the 3D area not flat area) located in Google Maps.
https://www.google.com/maps/place/Oujda,+Morocco/@34.4788585,-2.0182252,434m/data=!3m1!1e3!4m5!3m4!1s0xd7864984106d6ef:0x1d86b33244c4650!8m2!3d34.681962!4d-1.900155
Here its photo :


texts:
I have three rasters which represent distance (in metres) from certain things (bus stops, parks, libraries).  I would like to reclassify these rasters so that they get integer values based on the following criteria:

between 1 - 1000 = 1
between 1000 - 2000 = 2
between 2000 - 3000 = 3
between 3000 - 4000 = 4
between 4000 - 5000 = 5
between 5000 - 6000 = 6
between 6000 - 7000 = 7
between 7000 - 8000 = 8
between 8000 - 9000 = 9
9000 and above = 10

I know how to do this in Spatial Analyst, but I don't have a license and I'm trying to use Open Source tools.  I've attempted to use r.reclass and followed the tips in this question, but the input raster needs to be an integer raster.  At the moment, my distance rasters are floats.

I've also tried using the reclassify tool in the Geospatial Modelling Environment but the input raster also needs to be an integer raster.

I've also tried using the Raster Calculator in QGIS by following this example, but it appears that that only works with a single criteria not 10.

I've been pulling my hair out trying to convert the rasters to integer, but can't figure it out.

Does anyone know how to either reclassify a float raster or convert a float raster to an integer raster?

Note, the end goal of this will be to perform a weighted overlay.  Eventually, I may have Spatial Analyst in ArcGIS to do this, but I'd like to know how to do it using Open Source tools.

Thanks.


texts: Why won't this mapbox map load into a DIV container?

texts:
All,

For the life of me I can't figure out why I can't load a map into a CSS template I'm trying to customize.

I've placed a fiddle here:

http://jsfiddle.net/mapbaker/Du7zq/10/

Basically I'm trying to get a map to load into the red DIV Container (called mapBox)

I hope its something small and/or user error!

Thank you!

-m


texts: PostGIS topology: creation of topogeometry polygons for simplification

texts:
I am using PostGIS topology to manage a geologic mapping project.

My basic topology is defined by linework. The interstital polygons are meaningful, but can be constructed simply from lines and attribute data stored in points. I then generate the polygons "on the fly" from my linework, spatially joined to point attribute data. This allows me to keep the polygons in sync with my linework.

This works rather well, but it means that the polygons are Geometry objects, not Topogeometries. So I can't simplify the layer in a topologically-aware manner using St_Simplify.
I would also like to make larger polygons that are supersets of individual polygons based on attribute data, to allow generalized "dissolve" of map units.

Basically, it seems like I am looking for the equivalent of ST_GetFaceGeometry that returns a topologically-defined element. Perhaps topology.AddFace is what I need, but going from line topogeometries through polygonal geometries and back to a different topological representation of polygons seems confusing.

Plus I am not sure how to translate these face_ids to topogeometries and what happens if the linework is edited and edge_ids change – can I easily delete bad faces?

Any guidance from someone who has dealt with similar situations?


texts:
I received some raster data in .hts format from one of our clients. They use this data format in a tool called ATOLL.

Is there a way to convert this data to some standard format such as .bil or .img?


texts: Is there a FOSS solution for legend parsing in a web mapping environment?

texts:
Our web mapping requirements are fairly standard - some basic querying and the ability to print out the resulting map (or save the image). Typically we've used ESRI products in the past (right now we're using GeoCortex) but we still haven't been able to solve the problem of how to parse the legend information to match what is actually being shown on the map. For example, if there are four polygons showing up on the map, there should only be four legend items for the polygons. Is there an open source solution for this problem?

EDIT: Let's see if I can clarify this more. Currently, we have a geology layer with 33,000+ polygons. The legend for this layer is a single image and contains several hundred legend items (ie. it's the legend for the entire geology layer, stored in a JPEG/PNG).

Typically, a user will zoom in on an area (such as a single NTS map sheet, or even smaller) and so a subset of geology polygons will be shown on screen. They will then save the map as a print document (PDF or image), which is generated by GeoCortex/ArcGIS Server. The print document will contain a header, an image of the map and a legend. However, the legend for the geology layer will be the legend for the entire geology layer, not what is actually shown on the map (ie. a much smaller subset).

So, I'm looking to see if there is a FOSS solution to enable me to clip out or generate on the fly an subset of the geology legend which reflects the actual items on-screen/on the map, to go with the print document. I hope that clarifies things; I apologize if it doesn't!


texts:
I've just installed QGIS 2.0.1 on Ubuntu and the above error message appears.

Does anyone know how to resolve this problem?  I think I need to change the permissions of the relevant folder, but I'm not certain.


texts: How can I group Points to make polygon via Python

texts:
I have an array of coordinates

Point= [lat,long]

ListPoint= [[38.203655,-76.113281],[38.10,-76.12],..]

I wanna make a polygon that contains all these points ? There is a project developed by Flicker - The shape of Alpha but the project is not open.

Do you have any ideas how can I do that.


texts:
I want to calculate a cost distance model between a source of raw material and a group of sites. Raster and vector data are in the same CRS (EPSG 32615). I was able to use r.walk and get my cost and direction rasters. Whenever I want to use r.drain, I get:

Converting outputs
Loading resulting layers
The following layers were not correctly generated.
Least cost path
You can check the log messages to find more information about the execution of the algorithm

This is the log:

GRASS GIS 7 execution commands
            g.proj -c proj4="+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs"
            r.external input="C:/Users/gejst/Desktop/COst/rwalk/rwalk.tif" band=1 output=tmp1544491016043 --overwrite -o
            r.external input="C:/Users/gejst/Desktop/COst/rwalk/dir.tif" band=1 output=tmp1544491016044 --overwrite -o
            g.region n=2103521.41379 s=1879836.37719 e=394754.028731 w=180548.446796 res=30.0935068748
            r.drain  input="tmp1544491016043" direction="tmp1544491016044" start_coordinates=312000,1961000 output="outputef66d40f601f408cba4515b177ffa77b" drain="drainef66d40f601f408cba4515b177ffa77b" --overwrite
            g.region raster=outputef66d40f601f408cba4515b177ffa77b
            r.out.gdal --overwrite -c createopt="TFW=YES,COMPRESS=LZW" input=outputef66d40f601f408cba4515b177ffa77b output="C:\Users\gejst\AppData\Local\Temp\processing111bde5bc84d4563a5b91f7eeea042c8\de27ee9605134eaf96c3b82d9c8de29f\output.tif"
            v.out.ogr -s -e input=drainef66d40f601f408cba4515b177ffa77b type=auto output="C:\Users\gejst\AppData\Local\Temp\processing111bde5bc84d4563a5b91f7eeea042c8\2d50226dc6254d79b65c557fce716a27" format=ESRI_Shapefile output_layer=drain --overwrite

I am aware of these questions:

No output generated with "r.drain" GRASS module

How to get r.drain to run in QGIS/GRASS?

and I think I have not make mistakes in the coordinates of the  points (brackets, 0s)

In short I do not know what I am missing.I tried to run r.drain in Qgis 2.18.26 (stand allone install), 3.4 (osgeo4w install) running in win10 64.

Any hint is welcome


texts: Merging two different polygons to single polygon

texts:
I have a raw data with number of polygons and sub polygons in it. Attached Snap is an example.

I can simply merge two polygons by selecting both. But i can't do that to the huge data. Is there any way to make it dissolve the smallest one in largest one.


texts: Extracting random point coordinates using ArcGIS 10 model builder

texts:
I have a question about using model builder in ArcGIS 10.1.
In my model, I am creating a processing extent, and within that I am creating a random point to be used in further processing.
This model is iterative, and for each iteration, I am creating a new random point (effectively overwriting the values for the point).
What I would like to be able to do is to write the x-y coordinates of these random points off to a file before the next iteration is performed (either all together or as separate files for each point).
In a way, what I want is the inverse of the "Get Cell Value" tool and to be able to write that value to a file.
I guess it wouldn't be too difficult in python, but I was hoping to be able to do this in model builder.


texts: Converting string from GetParameterAsText to uppercase?

texts:
I am trying to make a tool in ArcToolbox using a Python script so users can input a street and lot number and it will zoom to the correct location/s in the cadastre. It's working but I have to type the name of the street in upper case in the toolbox and I'd like it work for uppercase or lower.

I assume I should use something like .upper() but I don't know where to put it. I've tried and failed in a few places.

Here's my snippet:

RoadName = "'" + arcpy.GetParameterAsText(0) + "'"


texts: ArcGIS runtime SDK 10.1.1 for java does not display anything

texts:
I am currently trying to display a map (using OpenStreetMapLayer) and nothing shows up but the esri logo and the OpenStreetMap copyright. I am not getting any error message even if the event mapReady never occurs.

I am running the application on a distant server (via Remote Desktop) on Windows Server 2008 R2 Enterprise (SP1) and using the jre 1.7.0_25 x86. Although the SDK version is 10.1.1, an ArcGIS 10.2 application is installed on the server (ArcMap and others).

Another strange thing is that even the SDK samples seem to encounter the same issue when run in the same conditions (I cannot try on another PC). See the screenshot below:

I have been searching the web for two days but found nothing and, as there are no error messages, I ran out of ideas to google.


texts: Create difference TIFF image from two TIFFs

texts:
I am trying to use gdal_calc to calculate the differences between two almost identical images. The problem is that the images are TIFFs with three bands (RGB) and the result is difficult to interpret. I use --calc but I don't really understand the syntax to get a mas(?) with only the differences.

I used this expression:

 --calc=(A-B) -> get a gray image

I would like to have White in the diferences and NoData in the equal values but I can't use "!=", for example:

--calc=((A-B)!=0) --calc=255*((A-B)==0) --NoDataValue=0


texts: ArcGIS JSAPI Toggle Identify (AMD style)

texts:
I am trying to convert the Toggle Identify example to AMD style. However, I am stuck with an empty map. Debugging in Firebug doesn't help too much. Could anyone help me to find where I am wrong? The code is in here:

http://jsfiddle.net/neWnE/

I am trying to modify Blord Castillo's code from this thread:

Identify with Layer Toggling

Any support appreciated!!!


texts: Memory considerations for cloud based GeoStack

texts:
Some time ago I asked about hardware considerations of cloud based GeoStack implementation. As a follow up to this question

There are some good pointers in related question already, but they mainly refer to PostGIS.

So I wanted to ask you for some advice regarding memory allocation for a simple Web mapping Ubuntu box, based on PostGIS, GeoServer and OpenLayers/GeoEXT.

By simple in this case  I mean serving maximum 10-20 WMS layers, caching, getFeatureInfo, and getting legends. Underlying dataset will consists of ~3000 polygons of administrative areas. It will be hosted on cloud VPS - either AWS or rackspace.

I have three questions regarding this setup:

Does it make sense to keep PostgreSQL and Geoserver on two separate instances? Or is it better to focus on one, more powerful server?
Cloud servers are scalable. So assuming 100% is the total amount of memory available - what would be your suggestions of setting memory parameters in PostgreSQL? And how about java memory allocation for GeoServer?
Any other hints for tweaking memory?

Thanks for help.

PS I asked this question on GeoServer Users mailing lists but haven't recieved any suggestions so far.


texts:
I'm using the QGIS ImportPhotos plugin to import photos taken on an iPhone and the coordinates of each photo displays as points based on the lat/long as expected.

However what I would like to do is display the direction each photo was taken in using a simple arrow. I believe that the direction information is recorded by the iPhone geotagging based on this response: https://stackoverflow.com/questions/3930266/what-information-is-stored-in-efix-jpeg-photos-taken-on-the-iphone-with-geotaggi

However when I look at the data the only direction attributes that come through is the Azimuth - is there anyway this can easily be converted to North or is the direction data not included in the data that the plugin collects from each photo?


texts: How to rename PostgreSQL/PostGIS columns while importing Shapefiles?

texts:
I discovered shp2psql to successfully import geospatial data from ESRI Shapefiles into PostgreSQL/PostGIS.
Now I would like to automatically rename the target column names. Is there any mechanism where I can define the original column names and their corresponding target column names in something such as a look-up table? Or would you do the renaming within the database using triggers or rules? Or do you recommend another import tool?


texts: How to set up ArcGIS Server 10.1 test environment on Mac laptop with VMware?

texts:
In my production environment I have ArcGIS Server 10.0 Enterprise, ArcSDE 9.3, and ArcIMS. I want to shut down ArcIMS as soon as possible for obvious reasons, and upgrade ArcSDE 9.3 to 10.1 and ArcGIS for Server 10.0 to 10.1. As I am new to ArcGIS for Server 10.1, I would like to set up a test environment just to get to know the software architecture but will not be able to do this on my network at work. I would therefore like to do this on my laptop with VMware.

Could anyone guide me as to how I should do this on my Mac with VMware?


texts:
I have GDAL 2.4.2 plus ECW plugin.

gdalinfo --formats
>>ECW -raster- (rw): ERDAS Compressed Wavelets (SDK 3.x)

And gdalinfo shows ecw file info.
But when i tried to use gdal2tiles i got error:

"C:\Program Files\GDAL\gdal2tiles.py" -r near -a 0,0,0 -z 15 D:\OFP\input.ecw D:\dev\nginx-1.12.1\www\OFP

ERROR 4: `D:\OFP\input.ecw' not recognized as a supported file f
ormat.
Usage: gdal2tiles.py [options] input_file [output]

gdal2tiles.py: error: It is not possible to open the input file 'D:\OFP\input.ecw'.

What can be wrong here?

UPDATE

I found that gdal module do not know about this format:

from osgeo import gdal

driver_list = []
for i in range(gdal.GetDriverCount()):
    driver = gdal.GetDriver(i)
    driver_list.append(driver.GetDescription())

print(str(driver_list))

['VRT', 'DERIVED', 'GTiff', 'NITF', 'RPFTOC', 'ECRGTOC', 'HFA', 'SAR_CEOS', 'CEO
S', 'JAXAPALSAR', 'GFF', 'ELAS', 'AIG', 'AAIGrid', 'GRASSASCIIGrid', 'SDTS', 'DT
ED', 'PNG', 'JPEG', 'MEM', 'JDEM', 'GIF', 'BIGGIF', 'ESAT', 'BSB', 'XPM', 'BMP',
 'DIMAP', 'AirSAR', 'RS2', 'SAFE', 'PCIDSK', 'PCRaster', 'ILWIS', 'SGI', 'SRTMHG
T', 'Leveller', 'Terragen', 'ISIS3', 'ISIS2', 'PDS', 'PDS4', 'VICAR', 'TIL', 'ER
S', 'JP2OpenJPEG', 'L1B', 'FIT', 'GRIB', 'RMF', 'WCS', 'WMS', 'MSGN', 'RST', 'IN
GR', 'GSAG', 'GSBG', 'GS7BG', 'COSAR', 'TSX', 'COASP', 'R', 'MAP', 'KMLSUPEROVER
LAY', 'PDF', 'Rasterlite', 'MBTiles', 'PLMOSAIC', 'CALS', 'WMTS', 'SENTINEL2', '
MRF', 'PNM', 'DOQ1', 'DOQ2', 'PAux', 'MFF', 'MFF2', 'FujiBAS', 'GSC', 'FAST', 'B
T', 'LAN', 'CPG', 'IDA', 'NDF', 'EIR', 'DIPEx', 'LCP', 'GTX', 'LOSLAS', 'NTv1',
'NTv2', 'CTable2', 'ACE2', 'SNODAS', 'KRO', 'ROI_PAC', 'RRASTER', 'BYN', 'ARG',
'RIK', 'USGSDEM', 'GXF', 'NWT_GRD', 'NWT_GRC', 'ADRG', 'SRP', 'BLX', 'PostGISRas
ter', 'SAGA', 'IGNFHeightASCIIGrid', 'XYZ', 'HF2', 'OZI', 'CTG', 'E00GRID', 'ZMa
p', 'NGSGEOID', 'IRIS', 'PRF', 'RDA', 'EEDAI', 'EEDA', 'SIGDEM', 'GNMFile', 'GNM
Database', 'DB2ODBC', 'ESRI Shapefile', 'MapInfo File', 'UK .NTF', 'OGR_SDTS', '
S57', 'DGN', 'OGR_VRT', 'REC', 'Memory', 'BNA', 'CSV', 'NAS', 'GML', 'GPX', 'LIB
KML', 'KML', 'GeoJSON', 'GeoJSONSeq', 'ESRIJSON', 'TopoJSON', 'Interlis 1', 'Int
erlis 2', 'OGR_GMT', 'GPKG', 'SQLite', 'ODBC', 'WAsP', 'PGeo', 'MSSQLSpatial', '
PostgreSQL', 'MySQL', 'OpenFileGDB', 'XPlane', 'DXF', 'CAD', 'Geoconcept', 'GeoR
SS', 'GPSTrackMaker', 'VFK', 'PGDUMP', 'OSM', 'GPSBabel', 'SUA', 'OpenAir', 'OGR
_PDS', 'WFS', 'WFS3', 'HTF', 'AeronavFAA', 'Geomedia', 'EDIGEO', 'GFT', 'SVG', '
CouchDB', 'Cloudant', 'Idrisi', 'ARCGEN', 'SEGUKOOA', 'SEGY', 'ODS', 'XLSX', 'El
asticSearch', 'Walk', 'Carto', 'SXF', 'Selafin', 'JML', 'PLSCENES', 'CSW', 'VDV'
, 'GMLAS', 'MVT', 'TIGER', 'AVCBin', 'AVCE00', 'NGW', 'GenBin', 'ENVI', 'EHdr',
'ISCE', 'HTTP']


texts: Get OpenLayers.Control.MousePosition in Toolbar

texts:
I have created an OpenLayers.Control.MousePosition:

var coord = new OpenLayers.Control.MousePosition({
            prefix: 'You are here: ',
            separator: ' | ',
            numDigits: 5,
            emptyString: 'Mouse is not over map.',
            displayProjection: "EPSG:4326",
            });

        map.addControl(coord);

And a Toolbar:

  var toolBar = new Ext.Toolbar({
            items: [
                 {
                  xtype: "label",
                  text: "Scale = 1:  ",
                 },
                 {
                    xtype: 'displayfield',
                 }
                ]
            });

This Toolbar is the tbar in my mapPanel:

var mapPanel = new GeoExt.MapPanel({
          region: "center",
          map: map,
          tbar:toolBar
          });

Is there any posibillity to get the Value of the Mouse Position as a displayfield in my Toolbar?
I already tried this: http://trac.geoext.org/attachment/ticket/195/195.1.patch,
but my map is not working with "  new GeoExt.Toolbar.ControlDisplay". Is there anything i can use instead?

I am using: OpenLayers 2, GeoExt and ExtJS.


texts:
I try to create a common extent for rasters in QGIS 2.99. Here is my script:

import processing

processing.runalg("script:unifyextentandresolution",
                  "/qgis_data/rasters/Image2.tif;
                  /qgis_data/rasters/Image1.tif",
                 -9999,"/qgis_data/rasters",True)

However, processing.runalg is being used for QGIS 2.99. I need it for QGIS 2.99


texts:
I've got a ImageMosaic (GeoServer) created from 4 band *.sid files. Images become transparent around the outer edges of *.sid files . When 4-th band is removed with gdal_translate -b 1 -b 2 -b 3 -of /X/Y.sid /X/Y_2.sid (gdal of course says Warning: The target file has a 'sid' extension, which is normally used by the MrSID driver, but the requested output driver is GTiff. Is it really what you want? link to a question about this)

When gdal translates sid to geotiff the images are displayed as expected but overhead storage required for geotiff files is unwanted.


texts: Are there any more specific tutorial about the -t_srs parameter in gdalwarp?

texts:
I want to use python to operate a seriel remote sensing images.Though the gdal official site http://www.gdal.org/gdalwarp.html show a explaination, it's hard to apply because no grammar tutorial for each command.

For instance, this demo:

  "gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84 " -overwrite raw_spot.tif utm11.tif".

What's the means of code "+proj=utm +zone=11 +datum=WGS84"? Are there any practical examples?


texts: Calculating Latitude/Longitude X miles from point?

texts:
I am wanting to find a latitude and longitude point given a bearing, a distance, and a starting latitude and longitude.

This appears to be the opposite of this question (Distance between lat/long points).

I have already looked into the haversine formula and think it's approximation of the world is probably close enough.

I am assuming that I need to solve the haversine formula for my unknown lat/long, is this correct? Are there any good websites that talk about this sort of thing? It seems like it would be common, but my googling has only turned up questions similar to the one above.

What I am really looking for is just a formula for this. I'd like to give it a starting lat/lng, a bearing, and a distance (miles or kilometers) and I would like to get out of it a lat/lng pair that represent where one would have ended up had they traveled along that route.


texts: Create buffer in meters in WGS84 layer using PyQGIS

texts:
In my QGIS 3 plugin I have this code that returns the list of all features intersected with a buffer calculated from two input values: the ID of the feature in the center of the buffer and a radius length. The problem is that the layer is in CRS WGS84, that uses degrees as unit. I need to calculate the buffer using meters as unit.

Is there a way to calculate this buffer in meters without change layer's CRS?

expr = QgsExpression("id = '"+InputID+"'")
it = layer.getFeatures(QgsFeatureRequest(expr))
ids = [i.id() for i in it]
layer.selectByIds(ids)
selection = layer.selectedFeatures()
for feat in selection:
    geom = feat.geometry().centroid()
    buff = geom.buffer(inputDistance,5)
intersectedFeaturesList=[]
for feature in layer.getFeatures():
    if feature.geometry().intersects(buff):
        intersectedFeaturesList.append(str(feature.attributes()[0]))
return intersectedFeaturesList

EDIT: I solved in this way, but I'm open to new suggestions.

import processing
from processing.core.Processing import Processing

Processing.initialize()
parameter = {'INPUT': layer, 'TARGET_CRS': 'EPSG:3857','OUTPUT': 'memory:Reprojected'}
result = processing.run('native:reprojectlayer', parameter)
layer3857 = result['OUTPUT']
expr = QgsExpression("id = '"+InputID+"'")
it = layer3857.getFeatures(QgsFeatureRequest(expr))
ids = [i.id() for i in it]
layer3857.selectByIds(ids)
selection = layer3857.selectedFeatures()
for feat in selection:
    geom = feat.geometry().centroid()
    buff = geom.buffer(inputDistance,5)
intersectedFeaturesList=[]
for feature in layer3857.getFeatures():
    if feature.geometry().intersects(buff):
        intersectedFeaturesList.append(str(feature.attributes()[0]))
return intersectedFeaturesList


texts:
I cannot get arcpy.Integrate_management to work.  I can run as a standalone process in ArcMap just fine, and if I copy the result as a python snippet I get:

arcpy.Integrate_management(in_features="SFMTA_BikeLanes #", cluster_tolerance="30 Feet")

Now, I try to modify this to use in my Python script, so that the process runs on a feature class that exists in my geodatabase.  I've tried:

#Set tolerance
ftTolerance = "30 Feet"
#arcpy.Integrate_management("G:\Shared drives\Bike Facilities Data\SF_BikeLanes.gdb\SFMTA_BikeLanes #", ftTolerance)

I've also tried:

#arcpy.Integrate_management("G:\Shared drives\Bike Facilities Data\SF_BikeLanes.gdb\SFMTA_BikeLanes", ftTolerance)

And, also tried:

# Define bike lanes feature class
bikeLanes = r'G:\Shared drives\Bike Facilities Data\SF_BikeLanes.gdb\SFMTA_BikeLanes'

arcpy.Integrate_management(bikeLanes,ftTolerance)

With all these, the script fails with the following error:

 Traceback (most recent call last):
  File "G:\My Drive\GIS\Bike_Aerial_Map_Automation\Scripts\Automation_Test.py", line 184, in <module>
    arcpy.Integrate_management(bikeLanes, ftTolerance)
  File "c:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\management.py", line 2133, in Integrate
    raise e
RuntimeError: Object: Error in executing tool

Failed to execute (TestAutomation).

Can anyone tell me what I'm doing wrong?


texts: Adding an empty dynamic layer makes my shapefiles display correctly - without it they are incorrect

texts:
I have a completely empty dynamic layer

public class CustDynamicLayer2 : BaseDynamicLayer
{
public override void DrawDynamicLayer(esriDynamicDrawPhase DynamicDrawPhase, IDisplay iDisplay, IDynamicDisplay iDynamicDisplay) { }
}

I add a shape file with these four calls:

gcs = srFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);
mapControl.SpatialReference = gcs;
mapControl.Map.AddLayer(new CustDynamicLayer2()); //without this my shapefile is badly distorted
mapControl.AddShapeFile(fileDir, fileName);

Basically if I don't add that empty class CustDynamicLayer2, the shapefile is distorted.  This makes no sense to me because the class doesn't do anything.


texts: Leaflet JS: Tiles misaligned using WMS resource

texts:
in Leaflet JS I instantiated a WMS tile layer object, given a base URL of the WMS service.
The wms service projection system is EPSG:32633 (WGS 84 / UTM zone 33N) and I had to use proj4js to project it properly on Leaflet object.

WMS url: webapps.sit.puglia.it/arcgis/services/BaseMaps/Ortofoto2013/…?

var crs = L.CRS.proj4js('EPSG:32633', '+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs');
var ortoOpt = {
    layers: ["0"],
    format: "image/png",
    transparent: true,
    request: "GetMap",
    version: "1.3.0",
    service: "WMS",
    crs: crs
  };
var ortofotoLayer = L.tileLayer.wms(wms_baseurl, ortoOpt);

Sadly, especially with higher zoom, the tiles appear misaligned, as in the picture.

It's possible to fix this?


texts: Selecting parcels in vector map from raster layer using ArcGIS for Desktop?

texts:
I am relatively new to Arcgis and I am stuck. I have a parcel map (vector) and a land cover map(raster). My land cover map is only 1's and 0's, flood and no flood. I want to be able to highlight or select all of the parcels from the vector layer that are majority flood (>50%).

I have scanned the forums for help, and I have tried at least half a dozen different options, raster to vector conversion, zonal statistics ect...But I can't seem to move forward. I think it's a query issue, right?


texts: CartoCSS Bubble map Jitter/Non-overlapping points

texts:
I am new to CartoCSS and web mapping in general. I am working on a project with different sized points that represent populations in different locations around cities. What I want to do is have it such that the points do not overlap, especially at changing zoom levels. I played around a bit with the marker-overlap true/false, but that doesn't seem to help. I've used Jitter functions to do something similar in R and other data visualization things, but from my cursory search I did not see anything about this. Do I need to write instructions for different zoom levels in relation to this? I am hoping not to use a cluster function because the data in each dot represents different ethnic groups, something I want preserved at a larger scale.

EDIT: The problem is not points that share the same location - rather, the points are close by but "smush" together into one blob at higher zoom levels. I also want to retain the bubble size effect. See images : https://imgur.com/gallery/m9dHg

EDIT: how to combine chloropleth, bubble and cluster designs in CARTO? <- this is a very similair problem, but again they have points that are on the same point rather than located close together.

Sample:

   /** Nepal */

#diasporaoutreachinfo [ _1_country = "Nepal"]{
  marker-fill-opacity: 0.9;
  marker-line-color: #FFF;
  marker-line-width: 1;
  marker-line-opacity: 1;
  marker-placement: point;
  marker-multi-policy: largest;
  marker-type: ellipse;
  marker-fill: #41006D;
  marker-allow-overlap: true;
  marker-clip: false;
}
#diasporaoutreachinfo [ _1_country = "Nepal"][ _3_count <= 5700] {
  marker-width: 25.0;


texts: How to determine in which grid cell each coordinate of an x,y series is located?

texts:
I have a circular arena which I have divided into a number of equal-sized 'cells' using a grid-like arrangement of concentric circles and line segments (in grey below) and also a series of x,y coordinates (plotted as a line in black below).

I was wondering what would be the best method (preferably using R) to find out which 'cell' each point is in ? I guess I could test to see which circle each point lies in by its distance from the centrepoint, and then which sector the point lies in by its angle from the centrepoint, but this seems rather laborious and inelegant..

This question bears some similarity to this one, except that I'm using circular lines to demarcate grid cells rather than linear lines.


texts:
I have created a network dataset from a road segment file of a city available through their data platform.
If I use the default settings and the only cost available is Length then the Service Area analysis works as expected. The rings created measure distance along the network for the set amounts of cutoffs at 5, 10, 15 meters respectively.
However, if I add Travel_Time (minutes) as a cost, then the service area created for 5, 10, 15 minutes are too small.
For calculating the Travel_Time, I have created a column by the same name in the road segments layer (before creating and building the network dataset) and used the Shape_Length (in meters) and SPEED (in kmph) attributes as follows:
Travel_Time = ((!Shape_Length!/1000)/!SPEED!)*60

After creating the network dataset, I have added new cost using this field and rebuild the network.
I am unsure why adding length as a cost makes the analysis not work as expected. If I increase the the cutoffs to 50, 100, 150, then the size of the Service Area increases. So, I am guessing it is not an issue with network connectivity.
I have already gone through the answer to Service Area Analysis and Time cost calculations in ArcGIS Pro.


texts: Offline Maps API for C# Desktop Application?

texts:
I have a C# WPF application which has browser control to display a Google Map.

Google Maps API v3 is being used in the application.

I use data received from a service to plot users current location and history. So a lot of markers and polylines are in use.

But to use Google Maps API, we need to be connected to the internet for sending requests to the server. And I believe that the API wont work offline even with the Enterprise Edition.

I need an alternative to the Google Maps API (preferably open source) to work completely offline, which provide similar services like marking locations, drawing lines, zoom, etc., I don't have lots of time to finish this. So instead of trying out the different available APIs, I was hoping somebody must have done a similar thing so I can finish my task on time. I am currently trying out ArcGIS. Is it good enough for my purpose? Is there an open source alternative?


texts: Create shapefile from data PostGIS or KML using Java

texts: I'm working on a geographic information system project, whose products form shapefiles. I am using the Java programming language and a PostGIS database. Is there a way to generate a shapefile from PostGIS or KML using the Java language?


texts: VB.Net sample for Google APIs Client Library for .NET

texts:
I'm just getting started on a project to replace Google Map's depreciated dynamic KML support.

I'm think Fusion Tables is the way to go. I have data that is updated daily - a map with several regions and I need to update the Color column with a daily forecast.

We currently have a six day's of forecast HTML w/KML sent to our site. In my new scheme I'm thinking I will still have six web pages but they will point to six Fusion Tables. Daily (almost, weekends excluded) we'll run the code to update the color values on the Fusion Tables.

Here's our current site: http://slocleanair.org/air/AirForcasting_map3.php  The page displays the next day forecast at the top of the page.

My normal toolset is VS 2013 Express and VB.Net.

I found the Google APIs Client Library for .NET project page and downloaded the binaries for the API. The sample link is currently throwing a 403 access error but my guess is the project is C# and overly complex and not for Fusion Tables. So I was wondering if someone can point me to/post a VB.Net sample for updating Fusion Tables from VB.Net?


texts:
I'm looking for a way to implement visibility (aka viewshed). I found a previous GIS SE question that suggested a few GIS applications (not what I need - I'm trying to embed the calculations), and also SAGA. SAGA looks more like what I'm looking for (C++, Java, C# classes), but I can't understand enough of the API documentation to figure out what classes I even need to instantiate (or subclass).

I'm not tied to the idea of SAGA - I'm just looking for library / engine that can do intervisibility calcs.

The constraints are that it needs to be:

reasonably "light weight" (since I'd like to be able to cover embedded / mobile)
reasonably open source (since my application will be open source, although I'm flexible on exactly which license as long as its (L)GPL compatible).

Can anyone recommend such a library or engine, and provide a description of or link to which bit of the library or engine I need to use to viewshed calculations?

Alternatively / in-addition, I'd also appreciate references to papers or tutorials that explain how to do these calculations in an efficient way (as applicable to embedded / mobile devices, so a GPU based implementation may not be as useful as something that is more general, and handles low power devices)

My preference is an existing implementation, rather than creating something myself (which is a a backup option).

Edit: C++ isn't a firm requirement - anything reasonably portable (C#, C++, Java) will do. I'm trying for a library or templates rather than embedding a large application.


texts: I get "The following layers were not correctly generated." error (Raster -> Extraction -> Contour) QGIS

texts: No matter what parameters I enter for the contouring, I always get the response in the log of "The following layers were not correctly generated." I'm working QGIS 3.6. I am not familiar with this program at all. What can I do to solve this?


texts:
is it possible to create pyramid from large rasters and add them with folders to google drive and mount google drive as local drive to geoserver ? i tried to search in internet to find about it and found syncdrive and grive for linux.Can these tools be used in linux servers(like centos,fedora,...)?
geoserver rest api has this capability to read raster pyramids from google drive?


texts: GeoServer + MongoDB + Boundless: CQL not working

texts:
After adding a collection of jsonArray data to my mongodb database I wanted to filter it with CQL, but it doesn't seem to work. My query goes like this:

user.name like '%username%'

When I apply this filter all my points disappear.
The only CQL filter that seems to work is when I try an int value on a root property:

id_str = 123

Versions I am using:

mongodb: 3.0.2
geoserver: 2.7.0
geoserver-exts: 2.7-SNAPSHOT


texts: Creating flat end type polylineline buffer in model without ArcGIS advanced

texts:
I have a polyline split into many short segments and would like to create a flat end buffer for each individual segment (not "dissolve" them all into one feature).

Is there any way to do this in an ArcGIS 10.2 model without purchasing ArcGIS Advanced licence?  Arcpy script, paid 3rd party extension, or open source extension etc?

Thanks!

Update: I ended up modifying the code found in this Q/A to create a custom script tool which:

creates an empty feature class for output
Iterates though every line of a different input feature class with many short polylines
at both the start and end of each line create a set of points which form a perpendicular line
connect those two "lines" (2 sets of 2 points) and feed them into the "shape" attribute of polygon
write those polygons to the empty feature class


texts: As per i know that there is an API called Openlayers.MapServer in openlayers 2, which is used to show MapServer file in browser in openlayers. How can i add an MapServer file in Openlayers 3, as i read somewhere that Openlayers 3 is currently under development for an API of MapServer, but please help me to add MapServer file in Openlayers 3.


texts:
Is it possible to specify the extent when creating an empty raster in GDAL?

I looked everywhere, but can't seem to find anything. The portion of my script where I write the raster would look something like this:

tifpath = "path/to/output/tiff.tif"
drvtif = gdal.GetDriverByName("GTiff")
rast = drvtif.Create(tifpath, Ncolumns, Nrows, Nbands, datatype)
rast.SetProjection(projection) # projection being a variable representing the desired projection
rast.GetRasterBand(1).WriteArray(array)

My array is an array of size Ncolumns x Nrows, filled with 0's (I want an empty raster). The number of pixels is predetermined, as is the cell size, and by multiplying npixels*cellsize, I've determined what my output extents need to be such that each pixel is a whole pixel (i.e. doesn't get cut off by the extent).

But, I can't find any options that would allow me to specify the desired extent.

It's also just occurred to me that I don't see any option to specify cell size.


texts: Attribute Table ArcMap / get specific values

texts:
I have an Attribute table with 151823 rows. I want to cut off the upper and lower 5% of the data so I sorted the Value field ascendingly but now I am not able to select the first and the last 7591 entries which would be the 5%, as the FID is not in the right order anymore. I tried it with a new field (NewID) with the intention to just have a new sequential numbering but this did not work. I used the script presented here

Sorting feature class to calculate sequential ID field using ArcGIS Field Calculator?

Does anyone know a simple solution?


texts: Single country coastlines in Python Matplotlib Cartopy

texts:
In the image below I am using cartopy to plot coastlines. But, I only wish to have the coastlines of Australia, and not the surrounding countries. Is there a way to remove the other countries.

Here is some sample code with the key steps in it for this problem:

fig  = plt.figure(figsize = (16, 16))
gs = gridspec.GridSpec(4,2, hspace = 0, wspace = 0)

for i in range(1,9):
    ax = fig.add_subplot(gs[i - 1], projection = ccrs.PlateCarree())
    ax.coastlines(resolution = '50m')
    ax.set_extent([110, 155, -23,-8.5])

.


texts:
I have an application that is able to import .lyr files into a Cesium map if they are bundled with all the required associated files in a zip file.  The application treats the bundled files in the .zip file as a .shp file with its associated files.  I want the application to import the .lyr file without having to bundle the .lyr file and other required files in a zip file.  Simply the user selects the .lyr file and the application will retrieve the associated files and plot items on the map.

Is there an easy solution to uploading .lyr files this way into Cesium?


texts: Where can I get land-use data for Nigeria?

texts: Anybody knows websites of land use for the Federal Republic of Nigeria? Free or paid. I am looking for the states of Kogi and Anambra.


texts: Is it possible to create a GPS track in QGIS and if so, how? (beginner question)

texts:
I'd like to use QGis as a trip planner. Once there's a topographic map on the screen - from whatever source - I'd like to point-and-click a track over that map, then export it in a common GPS format (GPX for instance).
Is that possible and if so, how?

Thanks.


texts: Removing Border from Layer List Symbology

texts: When I launched my webapp on our server I noticed these little black boundary boxes appeared around the symbology in my layer list. This is my first webapp deployed using WebAppDev so I'm fairly new with where to look to update things in the code. I've been able to update a few other things (color of text) by inspecting the code in the webapplication on a browser and tracing it back to the source from there but I'm not getting a straightforward answer on this one for some reason. I've messed around in the widgets->LayerList->css->style.css file (.jimu-widget-layerList .col-showLegend and .jimu-widget-layerList .showLegend-div) and tried adding a color and background-color to both of those but those attempts didn't render a succesful output.


texts: Creating counting grid from set of points using numpy?

texts:
I have a dataset of longitudes/latitudes as follows:

id,spp,lon,lat
1a,sp1,1,9
1b,sp1,3,11
1c,sp1,6,12
2a,sp2,1,9
2b,sp2,1,10
2c,sp2,3,10
2d,sp2,4,11
2e,sp2,5,12
2f,sp2,6,12
3a,sp3,4,13
3b,sp3,5,11
3c,sp3,8,8
4a,sp4,4,12
4b,sp4,6,11
4c,sp4,7,8
5a,sp5,8,8
5b,sp5,7,6
5c,sp5,8,2
6a,sp6,8,8
6b,sp6,7,5
6c,sp6,8,3

From such data, I want to generate a grid like this:

0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 1 1 2 0 0 0 0
0 0 0 0 0 1 1 1 1 0 0 0 0
0 0 0 1 0 1 0 0 0 0 0 0 0
0 0 0 2 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 3 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0

which gives the number of data records in each cell of the grid, using variable "spp" as a categorical (grouping) factor.

From this grid, I then want to create a heat map, superimposed on a geographical map, so that I end up with something like the figure below.

I can see how to plot a heatmap on a Matplotlib/Basemap, but I could not figure out how to generate the grid from the point data. Also, it is important that I am able to choose the grid size, so that several different resolutions can be evaluated.

Is there any way of doing this using Numpy meshgrid (or alternatively Scipy griddata)?


texts: GeoServer REST - Applying Style to Multiple layers

texts: I have a bunch of rasters in GeoServer that need to have the same style applied to them as their default style. Is there any way to do this?


texts: Assigning Raster Cell data to Point data

texts:
I have 23 sets of raster data for different types of habitat (which I have performed Block Statistics on in Arc10.2, in order to create a 1km grid). I also have 5 different Sets of point data showing presence data for a specific amphibian species. What I need to do is Assign the raster cell data to the point data which overlap. My final aim is to produce a table which will show the percentage of each habitat for each amphibian sighting

For example using great crested newt as one of my 5 species (I will end up with 5 tables) my Table headings could look like this,

GCN Sighting | % woodland | % grassland | % open water | ...% other 20 variables

From this I will produce box plots to illustrate habitat variability for each amphibian species, and hopefully will be able to draw some conclusions from it.

I am quite new to ArcGIS so please phrase it as if you were explaining this to a 5 year old!


texts: Find intersecting cell values for multiple rasters

texts: I am trying to find where 2 raster cells values intersect. I have a ground DEM raster and a water elevation raster.  I need to find where are values meeting and export the cells.  In the graphic below, the red circle is the cell that I want to find and export in ArcGIS 10.5 using raster calculator if possible.