texts: Zoom levels for WMS and WMTS layers are different

texts:
I use MapProxy to 'convert' a WMTS Layer to a WMS Layer. I need to do that because my client software can handle WMS only! I have a configuration which is working but there is a problem with the zoom-levels or the resolution. I have set up MapProxy with my configuration. I open then the demo site at localhost:8080/demo. Here I can view an OpenLayers map for the original WMTS Layer and the result of my conversion (MapProxy does a nice job here).

When I view the OpenLayers map of the original WMTS layer I can easily read all the labels for cities and streets. When I look at the map of the WMS layer it's hard to read these labels. The map doesn't appear sharp and the font size is too small. Furthermore it's obvious that the zoom levels for the OpenLayers WMS and -WMTS maps are different. It seems that the tiles for the WMS layer have shrunk.

Has someone a explanation for that behaviour? Is that a problem of MapProxy or a OpenLayers problem? Or is it just a constraint of the WMS specification?

Some additional informations:
The WMTS capabilities of the layer I try convert is available here: http://www.basemap.at/wmts/1.0.0/WMTSCapabilities.xml

A screenshot which shows the problem. On the left the MapProxy demo page for my WMS layer (labels are quite small and also the map extract is larger than in the original WMTS layer). On the right the demo page for the original WMTS layer (labels are large and sharp):
http://img42.com/VlzDo

The MapProxy configuration:

services:
  demo:
  wms:
    md:
    title: basemap.at WMS
  wmts:
    md:
    title: basemap.at WMTS

sources:
  basemap_wmts:
  type: tile
  url: http://maps.wien.gv.at/basemap/geolandbasemap/normal/google3857/%(z)s/%(y)s/%(x)s.jpeg
  transparent: false
  grid: webmercator
  coverage:
    bbox: [977650, 5838030, 1913530, 6281290]
    srs: 'EPSG:3857'

grids:
  webmercator:
  base: GLOBAL_WEBMERCATOR

caches:
  basemap_cache:
  format: image/jpeg
  request_format: image/jpeg
  grids: [webmercator]
  sources: [basemap_wmts]

layers:
  - name: basemap_layer
  title: Simple Basemaplayer
  sources: [basemap_cache]

OpenLayers for the WMS-layer map:

<script src="static/OpenLayers.js"></script>
<script type="text/javascript">
var map;
function init(){
    OpenLayers.Util.onImageLoadErrorColor = "transparent";

    map = new OpenLayers.Map('map', {
      maxResolution: 3655.78125,
      maxExtent: new OpenLayers.Bounds(977650, 5838030, 1913530, 6281290),
      projection: new OpenLayers.Projection("EPSG:3857")
    });
    var layer = new OpenLayers.Layer.WMS( "WMS basemap_layer",
        "../service?",
        {layers: "basemap_layer", format: "image/jpeg", srs:"EPSG:3857",
         exceptions: "application/vnd.ogc.se_inimage"},
        {singleTile: true, ratio: 1, isBaseLayer: true} );

    map.addLayer(layer);
    map.zoomToMaxExtent();
}
</script>

OpenLayers configuration for the WMTS-layer map:

<script src="static/OpenLayers.js"></script>
<script type="text/javascript">
var map;
function init(){
    var mapOptions = {
    projection: new OpenLayers.Projection('EPSG:3857'),
    resolutions: [156543.033928, 78271.516964, 39135.758482, 19567.879241, 9783.9396205, 4891.96981025, 2445.98490513, 1222.99245256, 611.496226281, 305.748113141, 152.87405657, 76.4370282852, 38.2185141426, 19.1092570713, 9.55462853565, 4.77731426782, 2.38865713391, 1.19432856696, 0.597164283478, 0.298582141739],
    units: 'm',
    maxExtent: new OpenLayers.Bounds(-20037508.3428, -20037508.3428, 20037508.3428, 20037508.3428)
};

map = new OpenLayers.Map('map', mapOptions);

var layer = new OpenLayers.Layer.WMTS({
 name: "WMTS basemap_layer",
 url: '../wmts/basemap_layer/{TileMatrixSet}/{TileMatrix}/{TileCol}/{TileRow}.jpeg',
 layer: 'basemap_layer',
 matrixSet: 'webmercator',
 format: 'image/jpeg',
 isBaseLayer: true,
 style: 'default',
 requestEncoding: 'REST'
});

map.addLayer(layer)
map.zoomToExtent(new OpenLayers.Bounds(977650.00, 5838030.00, 1913530.00, 6281290.00));
}
</script>


texts: Unexpected token in a sql query (Hibernate)

texts:
I have no more idea in how to run this query:

SELECT row_to_json(fc)::text
FROM (SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features
FROM (SELECT 'Feature' As type,
ST_AsGeoJSON(ST_Transform(lg.geom,4326))\\:\\:json As geometry,
(
select row_to_json(t)
from (select comune) t
)
As properties
FROM comuni_2015 As lg WHERE id >= 0 ) As f )  As fc LIMIT 10

This is what I'm trying to do:

return (Geomcomuni) getSession().createQuery(query);
                  "SELECT row_to_json(fc)::text \n" +
                  "FROM (SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features \n" +
                  "FROM (SELECT 'Feature' As type, \n" +
                  "ST_AsGeoJSON(ST_Transform(lg.geom,4326))::json As geometry, \n" +
                  "(\n"+
                  "select row_to_json(t) \n" +
                  "from (select comune) t \n"+
                  ")\n" +
                  "As properties \n" +
                  "FROM springgis As lg WHERE id >= 0 ) As f )  As fc LIMIT 10");

but it gives me this error:

  line 1:23: unexpected token: ":"

have you some suggestion to give me? or some other way..


texts: Iterating folder of shapefiles using ArcPy and exporting mxd to PDF

texts:
I'm new to Python and just created a hard coded script for the process I want, and it works. Right now it calls on the one of three shapefile I have in this map document. I want it to do the same thing, but for each shapefile in this folder. I have no idea how to go about this! I always use ModelBuilder, so I know that is should be an iteration.

Does anyone have any idea how to go about this, or how I am supposed to change this code?

import arcpy

#Declaring Variables
ID = arcpy.GetParameterAsText(0)
OutputLocation = arcpy.GetParameterAsText(1)

mxd = arcpy.mapping.MapDocument(r"C:\Users\m3rexkac\Desktop\Rename_Test\DOE_auto.mxd")

df = arcpy.mapping.ListDataFrames(mxd,"*")[0]

#Now to Create a Layer Object, the line of code below returns a list of all the layers in the dataframe, with this name.  The [0] returns the first layer with that name instead of the whole list
Layer = arcpy.mapping.ListLayers(mxd, "Mary_Creek_Co_Parcels", df)[0]

#This zooms the dataframe to the selected feature
lyr = arcpy.mapping.ListLayers(mxd,"Mary_Creek_Co_Parcels", df)[0]
extent = lyr.getExtent()

df.extent = extent
df.scale = scale

#Exporting map dataframe view to pdf, I have it set up to save
arcpy.mapping.ExportToPDF(mxd, OutputLocation + "2.pdf", "Page_Layout", 640, 480)

#Cleaning out memory by deleting the python objects
del mxd
del df
del Layer


texts: Can't clip vector layer from a raster layer I created

texts:
I've been trying to clip an area in a raster image which I have created from an OSM. However when I fill out the raster clipping form I get the following message.

Error 1: Output dataset. tif exists,  but some command line options were provided indicating a new dataset  should be created. Please delete existing dataset and run again.

I don't really get this because assuming I'm making a new layer with only the clipped part, I AM going to make a new dataset right?


texts: Load multiple csv files without defined x y headers

texts:
Initial Question Here Followed with Developments Below

I'm aware of the question posted here but my 120 csv files representing UK postcode centroids do not have x and y headers and dive straight into eastings and northings. Please see the code below which assumes coordinates are defined by columns 'x' and 'y'.

import glob, os

# Define path to  directory of your csv files
path_to_csv = "C:/Users/You/Desktop/csv folder/"

# Set current directory to path of csv files
os.chdir(path_to_csv)
# Find each .csv file and load them as vector layers
for fname in glob.glob("*.csv"):
    uri = "file:///" + path_to_csv + fname + "?    delimiter=%s&crs=epsg:4326&xField=%s&yField=%s" % (",", "x", "y")
    name = fname.replace('.csv', '')
    lyr = QgsVectorLayer(uri, name, 'delimitedtext')
    QgsMapLayerRegistry.instance().addMapLayer(lyr)

My data:

I'm only interested in the 1st, 3rd and 4th columns.
I've also attached a screenshot of what the data looks like in excel below the example code. How do I define columns 3 and 4 in python to replace 'x' and 'y' referred in this question?

AL1 1AG 10  515487  206498
AL1 1AJ 10  515491  206410
AL1 1AR 10  516270  205897
AL1 1AS 10  515005  206908
AL1 1AT 10  516131  206148
AL1 1AU 10  516141  206088
AL1 1BH 10  514836  206657
AL1 1BU 10  514643  206972
AL1 1BX 10  514619  206893
AL1 1BY 10  514582  206986
AL1 1BZ 10  514667  207047
AL1 1DG 10  514450  206754

Development 1.0

With a response to my initial question, I plugged the code in the editor but received the following error:

exec(open('C:/Users/JAMESB~1/AppData/Local/Temp/tmp3xtgxhmj.py'.encode('utf-8')).read())
Traceback (most recent call last):
File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\code.py", line 91, in runcode
exec(code, self.locals)
File "<input>", line 1, in <module>
File "<string>", line 4
SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 2-3: truncated \UXXXXXXXX escape

My Script:

import glob, os

# Define path to  directory of your csv files
path_to_csv = "C:\Users\>mynamehere<\Desktop\csv folder"

# Set current directory to path of csv files
os.chdir(path_to_csv)
# Find each .csv file and load them as vector layers
for fname in glob.glob("*.csv"):
   uri = "file:///" + path_to_csv + fname + "?useHeader=no&delimiter=%s&crs=epsg:27700&xField=3&yField=4" % (",")
   name = fname.replace('.csv', '')
   lyr = QgsVectorLayer(uri, name, 'delimitedtext')
   QgsMapLayerRegistry.instance().addMapLayer(lyr)

I have placed all the csv files into a folder called "csv folder" on my desktop. Can you help me solve my issue?
Sorry for my ignorance in python, this is a first for me!

Development 1.1

I corrected my file path with forward slashes and used the code below but still to no avail:

Script

import glob, os

# Define path to  directory of your csv files
path_to_csv = "C:/Users/myname/Desktop/csv folder/"

# Set current directory to path of csv files
os.chdir(path_to_csv)
# Find each .csv file and load them as vector layers for fname in
for fname in glob.glob("*.csv"):
    uri = "file:///" + path_to_csv + fname + "?useHeader=no&delimiter=%s&crs=epsg:27700&xField=3&yField=4" % (",")
    name = fname.replace('.csv', '')
    lyr = QgsVectorLayer(uri, name, 'delimitedtext')
    QgsMapLayerRegistry.instance().addMapLayer(lyr)

Error

exec(open('C:/Users/JAMESB~1/AppData/Local/Temp/tmp4ejnw3eo.py'.encode('utf-8')).read())
Traceback (most recent call last):
    File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\code.py", line 91, in runcode
        exec(code, self.locals)
    File "<input>", line 1, in <module>
    File "<string>", line 13, in <module>
NameError: name 'QgsMapLayerRegistry' is not defined


texts: How to use GDAL to convert Int16 data to Byte

texts:
I have some data with Int16 bands in HDF format.

I want to convert this to Byte bands in GeoTIFF format.

The data range for the Int16 is -2000 to 12000, so I want to map this to the Byte range 0 to 255.

How do I do this with GDAL tools? (if it is possible)
Otherwise how can I write a Python script to do this.


texts: How to close a polyline in QGIS while digitizing?

texts:
We need to work on a line layer where some features are closed and some are open. The closed features need to be closed perfectly, so that the first and the last point of the polyline perfectly match. In Kosmo GIS this can be achieved with right-mouseclick->close trace, in GvSIG also with right-mouseclick->close polyline (or simply with a keyboard entry "C"). I could not find a comparable functionality in QGIS. E.g. snapping to itself would be OK, but it seems not to work on features that are in the process of being created. Only if a feature is already added to the dataset, it becomes "snappable". Or am I missing something obvious?

Note that to divide the features between polygon and line layers is not an option for us. We need them to be on the same layer with the geometry type = line.


texts: Drawing polygon and extracting land cover inside using Google Earth Engine?

texts:
A new question derived from my Seeking website to do geoprocessing online?

In Google Earth Engine, is it possible (if so, how?) to draw a polygon and extract the land cover inside?

I would guess so, since it says in the Google Earth Engine API
 page:

  The functionality of Earth Engine is exposed through an API available in both JavaScript and Python. The API supports complex geospatial analyses including overlay, map algebra, array operations, image processing, classification, change detection, time series analysis, joins, raster-vector conversions, vector-based extraction of image statistics and much more. Algorithms are constantly being added, enhanced and updated. Through the API, users are free to script more complex analyses and creatively recombine existing algorithms. Reporting of results is supported through charting, mapping and table or image export.


texts: Selecting features within specific area in polygon

texts:
I'm working on a water utility network. My municipality is in another multipolygon feature. The attribute field "area name" defines other cadastral municipalities.
I want to select parts of the water network within a specific cadastral municipality.

I'm working in QGIS 2.18.


texts: Pandas dataframe to Shapely LineString using GroupBy & SortBy

texts:
I have a pandas dataframe that contains information to construct (poly)lines, and I want to use shapely & geopandas tools to make a SHP.
In the example below, I have 3 lines differentiated by "myid" and the order of the vertices is in "myorder."
Making shapefile from Pandas dataframe? is a great explanation for making a point shapefile, but I am looking for a polyline SHP. Creating Shapely LineString from two Points let's me know I need to use from shapely.geometry import LineString to make the polylines, but I don't understand from the answer there (nor the shapely documentation) how to specify groupby("myid") and sortby("myorder").
How would I do this?
Using Windows 10, Python 3.7.6, Conda 4.6.14.
myid = [1, 1, 1, 2, 2, 3, 3]
myorder = [1, 2, 3, 1, 2, 1, 2]
lat = [36.42, 36.4, 36.4, 36.49, 36.48, 36.39, 36.39]
long = [-118.11, -118.12, -118.11, -118.09, -118.09, -118.10, -118.11]
df = pd.DataFrame(list(zip(myid, myorder, lat, long)), columns =['myid', 'myorder', 'lat', 'long'])
display(df)


texts: Saving image with current styling using GDAL?

texts:
I want to automate some processes in raster images. Thus, one of the modules that I use is gdal_merge module of GDAL library in order to create RGB and False color satellite images.
So I execute:

gdal_merge.py -separate -n 0 -a_nodata 0 -o output.tif  band_4.tif band_3.tif band_2.tif

which works fine.
After that I want to render and save it with the current styling. Just like the option "save as rendered image" of QGIS. How can I succeed that with GDAL ?


texts: Way to Shrink Alaska and Hawaii for "Stereotypical View"

texts: Is there anyway to shrink Alaska and Hawaii to create a map that looks like the stereotypical map with AK in the corner, etc.? I thought about trying to isolate AK and HI in their own layers and then move those independently in the Print Composer but thought I'd see if there was an easier way.


texts:
If I develop QGIS plugins for my company do I have to make them available to others? I use QT-Designer and think that it also depends on the question how you develop your plugins, right?

For example:
If I'd work for public authorities they could be totally fine with publishing my plugins as there is no competition between them and because you don't have to gain a high ROI.

But: If I work for a private business that has a lot of competitors than my CEO could have a problem if I'd tell him that I develop a plugin for QGIS that every rival firm could use without contributing to the expenses we had for developing the plugin.
I indeed think that it is the duty of companies who use QGIS to contribute time or money to the QGIS project itself in order to keep the project "healthy" but this must be made apart from sharing plugin sourcecode in many cases.


texts: Conditionally control which objects are returned from overpass API

texts:
In generic form my question is: can I conditionally control what is returned from overpass query:

[out:json];
(
   rel[/* filter /];
   // if there are relations, stop and return them
   // if not, try to find ways and then relations that include them
   way[/ filter */];
   rel(bw); // returns relations having these ways
);
out body;

This is more specific example finding relation with San Jose city boundary:

[out:json];
area(3600396501)->.county;
(
   rel(area.county)[admin_level][admin_level!=6];
   // if relation found, return it. For San Jose city boundary
   // coincide with county boundary, so it's not within area :(
   // Find ways that are part of that relation:
   way(area.county)[admin_level][admin_level!=6];
   // one way would be sufficient, but how can tell it give me first?
   rel(bw)[admin_level][admin_level!=6];
)->.all;
(
   rel.all[name="San Jose"];
   rel.all[alt_name="San Jose"];
   // in this case I do not need ways, but for some cities the
   // boundary is just one way, without relation. Ideally if there are
   // relations I would stop here and return only relations.
   // For cities with boundaries as ways:
   way.all[name="San Jose"];
   way.all[alt_name="San Jose"];
);
out body;


texts: Errors in QGIS (about clip, editing shape etc.)

texts:
1) A few days ago, I run "clip" tool on some shapefile.

Green shapefile was obviously full of the area, but after clip process, some regions are lose their area...

Why does this error occur?

2) I edit shapefile in QGIS, for example... In editing mode, I'm edited some object use by tool like "divide object" & "merge selected object"

Before editing, shapefile show me image like lower jpg.

But after editing... QGIS show me image like below jpg.

Result contain sliver...

I run same process in ArcMap, and ArcMap give me a correct result...

Is it just an error in QGIS?


texts: I have a dataset with a linestring containing coordinates like: 376753.328, 2042051.611 I've tried unpacking using a couple different coordinate mappings to get the GPS coordinates. I'm told the coordinates are in India and I have a projection of EPSG:3857. I'm using primarily geopandas but use QGIS too. What else is needed to figure how to translate these into a more commonly recognized coordinate mapping


texts: PROJ4 Library error in combination with MapServer 6 and CentOS 6

texts:
I have recently compiled MapServer 6.2 on a CentOS 6.3 machine. Since I compiled it --with-proj=/usr, all the PROJ4 related stuff resides in /usr/share/proj.

I've also installed proj4 with yum install proj, which installed PROJ4-Version 4.7.0.1.

I wanted to establish a simple WMS Server but ended up having these "no system list, errno: 2 " errors when issuing a GetMap-Request.

Weirdly, my PROJ4-Installation did not establish any epsg-File in /usr/share/proj, so I downloaded the current source of PROJ4 and extracted the epsg file from there and copied it to /usr/share/proj/epsg - after that I am getting similar errors, just with errno: 20 instead of errno: 2.

Here is my mapfile: wms.map

This is the request I am getting errors with:

http://*.*.ch/wms?map=wms.map&SERVICE=WMS&VERSION=1.3&REQUEST=GetMap&LAYERS=Axis_NS&STYLES=&SRS=EPSG:21781&WIDTH=400&HEIGHT=300&FORMAT=image/png&BBOX=591800,185500,615000,201000

How can I debug this, and why did my PROJ4-Package not have an EPSG File included? Notice there is also no /nad-Folder, just two files nad27 and nad83 as well as nad.lst.


texts:
Vectorized few lines using QGIS with Google Maps as a base.
If I do not Trabachar with maximum zoom does not respect the topology of the base.
It is possible to make snap using Google Maps based on QGIS?


texts: Couldn't Excute command line PostGIS pgsql2shp

texts:
I am using postgres 9.3 with postgis 2.1 on Ubuntu 14. When I tried to dump pgsql table using pgsql2shp, I get the following error.

Error: pg_wrapper: pgsql2shp was not found in /usr/lib/postgresql/9.3/bin

The command I used is the following.

pgsql2shp -f "/home/user/testshp.shp" -h localhost -u postgres -P postgres  dbname "SELECT * FROM table"

In  sudo nano /etc/environment I have added,:/usr/local/bin/pgsql2shp

I also also created a link using a guide from here.

sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/shp2pgsql
sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/pgsql2shp
sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/raster2pgsql

Is any the things I did causing the problem?


texts: Merging values from multiple fields based on condition in Attribute table?

texts:
In QGIS attribute table of a point shapefile I have several rows with species  occurrences. e.g. BA1404.

Whenever a species is present, the cell value = 1 otherwise 0

I created the new column "BA" and I converted the species code to species names: e.g. BA1404 = code

Acer negundo = species name

case
when "BA1404"  = 1 then 'Acer negundo'
when "BA3501"  = 1 then 'Aescolus sp.'
when "BA2007"  = 1 then 'Ailanthus altissima'
when "BA2009"  = 1 then 'Celtis sp.'
when "BA1304"  = 1 then 'Fraxinus pennsylvanica'
when "BA2010"  = 1 then 'Gleditsia triacanthos'
when "BA2002"  = 1 then 'Juglans nigra'
when "BA0404"  = 1 then 'Pinus strobus'
when "BA2008"  = 1 then 'Platanus sp.'
when "BA3305"  = 1 then 'Populus balsamifera'
when "BA0500"  = 1 then 'Pseudozuga menzisii'
when "BA1105"  = 1 then 'Quercus rubra'
when "BA1700"  = 1 then 'Robinia pseudacacia'
when "BA3304"  = 1 then 'Populus x canadensis'
when "BA1800"  = 1 then 'Different species'
else 'Error'
end

Now the present species are in column "BA".
BUT, in some cases there are TWO or more species in one row (red circles in figure 1)

Then the result should be:

'Pseudozuga menzisii, Robinia pseudacacia, Populus x canadensis'

How can I resolve this?

@she_weeds: I tried with all the steps you suggested:

made a CSV table with codes and species and dragged/dropped it into the QGIS project

now the CSV file is a table

select the point layer -> Field calculator -> Update existing field (case sentitive fields are adjusted). But still it does not give me any output as you can see.

What is weird: when I choose the map layer from the layer list; layer name is different from actual layer name...


texts:
I installed the plugin "tweeter 2qgis", I need the Python library "tweepy".

Downloading the plugin, a dialog window appeared, saying "we will download tweepy to your temp folder, but you must install it by yourself!".

Ok, before several minutes I've finally found this folder, but the extention's file is .py and I can't open anything. I have Windows, maybe is because of this?


texts: How to pass a token value to QGIS processing modeler?

texts:
There are various places in the Processing modeler where I would like to pass the value of a variable.
We're iterating through a bunch of features in an input layer. After the processing, we will have a new shapefile per input feature.
I want to name each of the new shapefiles based on a value in the new attribute table...this gives me a unique name per new shapefile.
For example, using the Field Calculator, We run a simple math operation and store the result in a new result field.
We can create a new OutputVector for the calculated result, and I want to use an attribute table value to create the name, like:
myResultVector-%SOMECOLUMNLABEL

Where SOMECOLUMNLABEL is the column that contains the attribute of interest, and myResultVector is a nice namespace.
If I could name each of the output vectors uniquely, I could save them all in a single folder for future processing.


texts: Layer Order panel duplicating layers in QGIS?

texts:
I'd like to just confirm that a problem I am having in QGIS 2.10 is not unique to me.

I am using the Layer Order panel to control the rendering order of layers in my projects, it's very useful when working with lots of layers. However in QGIS 2.10, when I move a layer it does not change the layer's rendering order, it simply duplicates the layer (see below, notice all the duplicated layers).

The Value Tool also sees the duplicated layers (if it's a raster). The layers are not duplicated in the 'Layers' panel. It only happens when I tick the 'Control rendering order' tick box at the bottom of the Layer Order panel.

If others are affected I'll submit a bug report.


texts: Export Canvas extent inside QGIS at a 1m resolution GeoTIFF format

texts: I am working with the "XYZ Tiles Google Earth" inside QGIS and I have an area of interest that is roughly a 2 nautical mile radius that I would like to export as a GeoTIFF with a 1m resolution to use as a reference in another program. Currently in order to get it to the resolution I desire I have to zoom in to the scale I want and then export the canvas extent. The problem with doing it this way is that I would have to somehow keep track of where I am at and end up with a lot of small exports that I would have to stitch together later on. I was hoping if someone knew of a method in which I could export the entire area of interest I have at 1 meter resolution size scale. Or if anyone knows of an entire different method?


texts:
I'm attempting to follow along with this tutorial on using web images as markers within QGIS.  I'm using static web-linked images, not webcams as in the tutorial.

The code seems to be working correctly, except the svg that is generated does not render within QGIS.

This is the code i'm using in the symbology function editor, taken from the above tutorial:

import requests
import base64

@qgsfunction(args='auto', group='Custom')
def show_camera(feed, feature, parent):
    svg = """
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<svg version="1.1" baseProfile="full" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
  <g>
    <image xlink:href="data:image/jpeg;base64,{0}" />
  </g>
</svg>
"""
    data = requests.get(feed, stream=True).content
    name = feed[-16:]
    b64response = base64.b64encode(data)
    newsvg = svg.format(b64response).replace('\n','')
    path = r"C:\test\{0}.svg".format(name)
    with open(path, 'w') as f:
        f.write(newsvg)
    return path.replace("\\", "/")

When I attempt to open the generated svg file directly in a browser, it also does not render, just showing me a blank page.

Any suggestions?  Are there any better methods to use images as markers within QGIS in the 2 years since this tutorial was posted?


texts: How to join attributes by location with command line?

texts:
I would like to join the attributes of two shapefiles by location but with a command line (Linux).
This feature is easy in QGIS but it is VERY slow.

Is there a way of doing this with ogr or PyQGIS?


texts:
I am working on an XP machine with FME'13 and I have a Address field and I have to find the spaces(middle spaces and not the leading and trailing ) in each and every row.

St._Petters_Road is the valid space(it need not be flagged); but
St.__Petters_Road is invalid (as it has two consecutive whitespaces).

It is assumed that the words are correct. But there should not be two two consecutive spaces.

and I need to flag the invalid attributes in a new column as 0 and 1
Can we do that with transformers or in SQL?


texts: QgsProject.layerTreeRoot() contains no layers

texts:
I'd like to automate a layer rasterization process using PyQGIS but the first step---obtaining a list of the layers to rasterize---is failing in QGIS 2.18.9 on Windows 10.  What am I missing?  Regardless of which project it's pointed to or whether Python is running from an IDE or the OSGeo4W shell directly the test script below always indicates zero layers even though project.read() always returns true.

import sys
sys.path.append('C:/Program Files/QGIS 2.18/apps/qgis/python')

from qgis.core import *
from PyQt4.QtCore import QFileInfo
qgis = QgsApplication([], False)
# EDIT: setPrefixPath() resolves the issue per Germán Carrillo's answer below
qgis.setPrefixPath('C:/Program Files/QGIS 2.18/apps/qgis', True)
qgis.initQgis()

project = QgsProject.instance()
project.read(QFileInfo('C:/Users/<username>/Documents/<path>/<project>.qgs'))

layers = project.layerTreeRoot()
print len(layers.children()), "layers loaded"

qgis.exitQgis()

I'm aware of this question and the code above conforms to a number of other similar examples which are supposed to work.  Using QgsApplication([], True) or various other flavours of the above calls makes no difference.  QgsMapLayerRegistry.instance().mapLayers() is also empty.

If the same projects are opened in the regular QGIS UI querying:

QgsProject.instance().layerTreeRoot().children()

from the Python console contains all layers as expected.  Having the project open in QGIS proper versus not having QGIS running has no effect on the external Python python behavior.  So this doesn't look like contention between QGIS instances.


texts: How to add a label of a point to a polygon

texts:
i have several polygons with in each polygon a point with a label (number).

I want to add the label in the polygon to the attribute-table of the polygon it belongs to.


texts: Struggling to install QGIS from source in RHEL 6.4

texts:
I have installed QGIS from source in Ubuntu using the this link. In ubuntu I compiled qgis from source without any problem.

but now I moved to RHEL 6.4, not able to run QGIS from source. The steps I have followed to install QGIS from source are described below.

First, I have added all repository which required by QGIS. They are as follows:

devtools-1.1.repo
elgis-testing.repo
elgis.repo
epel-testing.repo
epel.repo
remi-php70.repo
remi-safe.repo
remi.repo
rhel-source.repo
rpmforge.repo
slc6-extras.repo
slc6-os.repo
slc6-updates.repo

and I also installed Python 2.7, QScintilla2, GLIBC 2.14 from source because of version mismatching.

Then,I followed this link to build the QGIS, till this stage everything is fine. Build completed without any error.

After launching the QGIS application, It is showing an error saying that
"Import Error:cannot import name Qgsfeature". this error is related to python versions, As I have two python versions installed my machine(Python 2.66 and Python 2.7). then I set the PYTHONPATH, But it is not affecting to the Application, So still displaying the same error.


texts:
I'm not a API person but a GIS user and I would like to export points from a category from Wikimapia to KML or any other useful GIS formats.

Has anyone developed that?

It could be interesting to add the analysis that we need to the data from wikimapia.


texts: Unzip a file using python caller in FME 2011

texts: The idea is to download a zipped file from a website and unzip it using FME 2011. I have attached a screenshot of FME workbench which is not working.


texts: jquery communication with arcgis services

texts: I am working with the ArcGIS JavaScript API for my application.  My app users are stored in a table published in a map service. Whenever a user wants to connect to the app, he needs to provide a username and password. How can I retrieve data from my table and compare them to those entered by the user? I've found somewhere that I need to use JQuery but i don't know how. If you could help, thanks in advance:)


texts:
I'm currently trying to translate the tasks I performed using Qgis (3.16.5) into a Python (Spyder) script, such that I will be able to perform these tasks on multiple files more easily.
However, I'm experiencing difficulties importing Qgis into Spyder. First, I received an error. After looking on internet, I found many topics related to this and I tried most of them.
Most of them were related to changing the PYTHONPATH and PATHs. I ended up with the script below:
'''
@echo off
  SET OSGEO4W_ROOT=C:\OSGeo4W64
  call "%OSGEO4W_ROOT%"\bin\o4w_env.bat
  call "%OSGEO4W_ROOT%"\apps\grass\grass78\etc\env.bat
  @echo off
  path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
  path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass78\lib

  set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python;
  set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\Python27\Lib\site-packages
  cmd.exe

However, I now get a new error: "DLL load failed while importing _core: Cannot find procedure" and still cannot import qgis.core. I tried this code both with OSGeo4W64 and with a Qgis version downloaded without OSGeo4W64.
How can I import the Qgis.core?


texts:
I have a physical copy of a map with railways that I want to import directly into ArcGIS (the network of railways), but I have no data on coordinates of any kind (it's a really old map but the lines are clear to see).

What would be the best way to go about this?

I'm a beginner.


texts: Is there a way to tie line width to map scale in QGIS?  All I found was an outside limit for display as the scale decreases.  I'm using an outline line style for roads, and the roads are too big when I zoom out.  What about converting them to a raster layer?


texts:
the standard method requires to run ArcPy script in command line, i've struggled with it to run but without a result, i can import arcpy, set the workspace to the geodatabase path and describe the env variable using arcpy.Describe(arcpy.env.workspace).name it gives me the geodatabase name, but when i try to get featureclasses using arcpy.ListFeatureClasses() i get [], i have made sure that all the environement variables in PATH and PYTHONPATH are ok .

i would like to know if there's an alternative way to run ArcPy script, like it could be run in ArcMap console, is there a way to send scripts to ArcMap console using ArcObjects, what's the difference between ArcMap console and standard command line console?

so any way to run an ArcPy script using ArcObjects 10 or Arcgis Engine 10 will be very welcome.


texts: How to make slope map from elevation map in Lat/Long?

texts: I have an elevation raster map of the world, in Lat/Long. When I make a slope map, the results of the raster cells will be either 0, or values close to 90. The problem seems to be that the Slope function does not understand the horizontal units (which is obviously degrees, but should be meters). How can I solve this? I am using ArcGIS.


texts: Matching d3.js geo projection with pyproj Proj4

texts:
I am having trouble figuring out the correct notation for Proj4 projections with pyproj. I got the following projection in d3.js and I want to match it with pyproj on my backend to send a background bitmap for the map.

var projection = d3.geo.albers()
    .center([0, 55.4])
    .rotate([4.4, 0])
    .parallels([50, 60])
    .scale(4800)
    .translate([width / 2, height / 2]);

Here is what it looks like:

The points I need to project onto this are lat, lng.

 toP = pyproj.Proj("+proj=aea ... ? ")
 x, y = pyproj.transform(EPSG4326, toP, x, y)


texts: Layers reference from Geoprocessing using ArcObjects?

texts:
I am using ArcObjects (ArcGIS 10.1) in .NET with some geoprocessing. It works great. There is an option to add the resulting feature class to the map. That works great too. However how do I get a reference to this layer to move it to a group layer?

Dim GP As GeoProcessor = New GeoProcessor
GP.AddOutputsToMap = True 'True if we want to add outputs to map, False if we dont.

Dim parameters As IVariantArray = New VarArray
parameters.Add(inDatasets)
parameters.Add(toPointFullName)
parameters.Add("")

Dim res As IGeoProcessorResult = GP.Execute("Merge_management", parameters, Nothing)
'todo: how do I get a reference to ILayer of the resulting layer added to map?


texts: Writing script to append Multiple GPX Files using ArcPy

texts:
I am using ArcMap 10.7.

How do I write a script to insert multiple GPX Files as a Feature Class?

What I was thinking is batching the files together in to one file and then running the code I already have for that one file.

import arcpy
arcpy.GPXtoFeatures_conversion('c:\\GPX_Files\\Hike.gpx', 'c:\\gisData\\Hike.shp')


texts: How to stick plugin widgets for ever to its own position?

texts: Every time I'm starting Qgis my custom widget disappearing. Is there a way to stick plugin to keep it all the time at the same place in the screen layout ?


texts: Buffer analysis of points, counting intersects of resulting polygons

texts:
I am using ArcGIS pro and this is the first time I've done something like this. I don't have a "proper" GIS background and have just been plotting point maps and choropleth maps.
I've been asked to plot some points on a map, draw a buffer of a fixed distance around them, then determine the number of overlaps for each buffer, and then determine which points have the most overlaps.

I think I do something like this:

1) draw buffer of 400m around points
2) self join this buffer layer to itself -
   I think this is a one to many join (as I want to see all the buffers a single
   buffer overlaps and I think I set this to have a match option of Boundary
   touches)

3) dissolve this joined layer so I can see which point buffer intersects with which point buffers and how many. each point has a unique ID. So if I dissolve by the field Target_FID (which is generated by the join) and the unique ID and sum the Join_count field, I should generate a layer where each point has a sum join count.

I'm having problems visualising it.

If I've not explained anything clearly, I'm really sorry, I've only really got a vague idea of what words to use to explain what I want.
What I want to do is
Draw some buffers on some points - can do this

Determine their overlaps with the other buffers on the layer - can do this, but do I use intersect, boundary touches, what set up? and I know I need one to many so that I can count ALL the overlaps

Calculate the number of overlapping polygons - this is probably the field Sum Join_count (but do I subtract one from this value as the polygon probably overlaps with itself?) can do this

Map this to determine the "catchment area" (is that the right term?) for each point, so that I can show for each point the number of overlaps and determine if there is any clustering or if there are any hotspots. Can sort of do this, I can change the symbology to graduated colours but at the moment they're complete circles, what I think I want is more intensity at the greatest point of overlap and more transparency where the overlaps are least. (if that makes sense?)

This document
https://www.newham.gov.uk/Documents/Environment%20and%20planning/EB01.%20Evidence%20Base%20-%20Cumulative%20Impact%20V2.pdf

Explains exactly what I want to do:
assessing the degree of spatial proximity of each point to other equivalent
points  by looking at the number of others within 400m (5 minute walk).

Is there a specific ArcGIS Pro tool to do this, is my workaround of self join, sum, display the right way?
I'm spinning my wheels generating "stuff" but not the "thing" I need.


texts: When using Raster Image Fill in QGIS  I have no problem getting the image there but is there any way to make it stay one size within the shapefile rather than changing depending on how far you zoom in and out?


texts:
I published a workspace on FME Server. The workspace converts Geodatabase Files to MapInfo TAB.  I'd like to call the generic workspace with FME Server to convert other geodatabase.

Any idea?


texts: What are these terms in Marine Navigation ETA, DTG and TTG?

texts:
Consider a scenario where a ship has to navigate and reach a specific destination. I have a route plan comprising of multiple route way-points through which my ship has to go through to reach its destination point. Throughout the navigation, I'm required to calculate the ETA (Estimated Time of Arrival), DTG (Distance To Go) and TTG (Time To Go).

What exactly are these terms and can you provide mathematical expressions for calculating each?


texts: Using Spatial Join to join geologic features to polylines

texts: I want to spatially join polygons (geologic features, major aquifers) to where the NHD polylines (rivers/streams) fall. I want to be able to click on a polyline and see what geologic feature or major aquifer the NHD polyline flows through it. Do I use the spatial join for this? I have joined other attributes like Stream Order, Velocity, Watershed area because I can join them based on a ComID because all information came from NHDplus. I am using ArcMap 10.4.


texts: How use DataFrame Geographic Transformation User Interface with ArcObjects?

texts:
I am bringing data from SQL Server into ArcMap with ArcObjects as a Query Layer.  If I do this manually, a 'Geographic Coordinate System Transformations' UI pops up.  I would like the same thing to happen when I use my code OR to be able to set the transformation automatically.  But I'm having trouble finding how to handle on the fly transformations with ArcObjects.  Can someone point me in the right direction?

Edit: This is a Python Addin that I am developing.  I'm using comtypes to get at the ArcObjects.  Please tell me if I'm wrong but it seems that the simplest way to set the transformations is with arcpy.mapping.DataFrame.geographicTransformations.  Then I can append to the list the transformations I want the dataframe to use, depending on which Query Layers are being asked for.

I would still like to have the same dialog that is used when manually setting the DataFrame's Geographic Transformations pop up after the data is brought into the map so the user is aware of what transformations are being used.  Does someone know if this is possible?

Edit2:  Maybe if I could trigger the event that warns the user that a featureclass with a differing projection than the dataframe has been added to the map.  Where would I go to start looking into this?


texts: Changing default settings for the Save Vector Layer as... dialogue box - QGIS

texts:
I would like to change the defaults selected in the Save Vector Layer as... dialogue box.  I use export to kml and Shapefile daily and would like it to default to the settings I use, rather than changing them every time.  Is there a way to do this?
I am using QGIS 3.14.1 Pi for Mac and QGIS 3.10.10 A Coruna for Windows.


texts: Using cursor to make all LONG fields with NULL values 0 in ArcPy?

texts:
I have a table that I am joining to calculate several long integer fields, when the join is executed all of the values are joined as NULL that don't have values.

How can I use an Update Cursor to calculate every NULL value to 0 that is in a field that has a long integer type?

How can I expand upon this example to iterate through each field that is not a string(double, long) and change all NULL values to 0?

This only needs to be done to one field; however the script below does not update my NULL values to 0.

import arcpy

distfc = 'aFC'

with arcpy.da.UpdateCursor(distfc, ["FREQUENCY"]) as cursor:
    for row in cursor:
        if row[0] ==  NONE:
            row[0] = 0
            cursor.updateRow(row)


texts: ArcGIS: Create a line layer from a point layer and csv data file

texts:
UPDATED DECEMBER 20

I have a Regular point layer and a csv file that has the same ids as the point layer. The csv file is organized as

ID     From    To     Weight
--     ----    --     ------
01     1234    1235   2
02     1234    1236   4
03     1234    1237   1
04     1235    1234   4
05     1236    1235   3

// etc

the point layer is organized in the following way:

MY_ID
-----
1234
1235
1236
1237

the From / To fields correspond to a my_ID in the point layer. I'd like to create a new line layer, such that the From, To, Weight correspond to both the point and csv files


texts: "Wrong or missing parameter value: Input layers"

texts:
I get this error message when I try p. ex. merge layers (native or SAGA algorithm) or package layers (processing) into  geopackage. When I set files from a directory I don't have this problem.

Am I missing some configuration? I've tried short paths, layernames without space...

Linux Mint 19.2 - QGIS 3.10.1


texts: Styling each feature of a layer individually

texts:
I’ve a PostGIS table which is published as a Geoserver layer. Every feature in the table has some attributes that should style that feature (e.g color, size, opacity, external graphic url, ... etc).

I can create an SLD document and let it read its styling attributes from the table using <PropertyName>, but the problem is more complex than this, for example: some points features need to be styled as a simple circles, while some other need to be styled as external-graphic images. And from the SLD Reference, the SLD PointSymbolizer must have either <Mark> or <ExternalGraphic> element, for that reason I should create two rules one for the simple point symbolizer and the other for the external-graphic point symbolizer, and should add a filter to each rule so that it deals only with the appropriate features, something like:

<FeatureTypeStyle>
    <Rule>
        <ogc:Filter>
            <ogc:PropertyIsEqual>
                <ogc:PropertyName>is-external-graphic</ogc:PropertyName>
                <ogc:Literal>0</ogc:Literal>
            </ogc:PropertyIsEqual>
        </ogc:Filter>
        <PointSymbolizer>
            <Graphic>
                <Mark>
                    <WellKnownName>circle</WellKnownName>
                    <Fill>
                        <CssParameter name="fill">
                            <ogc:PropertyName>fill-color</ogc:PropertyName>
                        </CssParameter>
                    </Fill>
                </Mark>
                <Size>6</Size>
            </Graphic>
        </PointSymbolizer>
    </Rule>
    <Rule>
        <ogc:Filter>
            <ogc:PropertyIsEqual>
                <ogc:PropertyName>is-external-graphic</ogc:PropertyName>
                <ogc:Literal>1</ogc:Literal>
            </ogc:PropertyIsEqual>
        </ogc:Filter>
        <PointSymbolizer>
            <Graphic>
                <ExternalGraphic>
                    <OnlineResource
               xlink:type="simple"
               xlink:href="smileyface.png" />
                    <Format>image/png</Format>
                </ExternalGraphic>
                <Size>6</Size>
            </Graphic>
        </PointSymbolizer>
    </Rule>
</FeatureTypeStyle>

My question is: Is there any better solution that is more dynamic?


texts:
Can someone guide me for applying thresholding technique to input raster image using Python.

I intend to extract impervious surfaces from a high resolution raster such as roads, sideways and rooftops etc.


texts: Display multiple identical polygons (with different attributes) on map

texts:
I need to create maps of multiple species in a certain area. However sometimes multiple species occurr in the exact same area. This results in two identical polygons on top of each other. Many of these polygons are for example a UTM1 square of 1 by 1 km.

When making a map I should be able to make it clear that there are multiple species that occur in that polygon. However through the syles menu I can only make the colors blend while ideally when 3 polygons are on top of each other I should get some sort of pie diagramm or a checkerboard pattern containing the colors of both species.

Does anyone know a plugin or elegant way to do this? So far I have only found stuff like the v.clean option, but I do not want to remove the overlapping polygons, quite the opposite, I want to visualise the overlapping polygons as well as the non-overlapping ones with a clear legend.


texts: What's the difference between using PostGIS functions and operators?

texts:
PostGIS offers a list of operators implementing spatial relations. Some of them have a direct function equivalent, for example <-> and ST_Distance.

EXPLAIN
SELECT
st_distance(st_setsrid(st_makepoint(0, 0), 4326),
             st_setsrid(st_makepoint(10, 10), 4326))

and

EXPLAIN
SELECT
st_setsrid(st_makepoint (0, 0), 4326) <->
st_setsrid(st_makepoint (10, 10), 4326)

yield the exact same query plan. Are there any scenarios in which one method is preferable to the other?


texts: Which interface can replace IArcPadExtension?

texts: I trying to migrate a code written for arcgis desktop and I have there a IArcPadExtension, but here (http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#/type_changes_between_9_3_and_10/000100000408000000/) we can read that this Interface removed. Can someone tell me which interface can be used instead of IArcPadExtension?


texts: Selecting largest attribute of repeat attributes using ArcGIS Desktop?

texts:
I am have a feature class of thousands of cities, however sometimes the cities are not a single polygon and instead have small fragments separate, with their own object ID, but the same name field.

I was trying to label the cities so that it would only pull the largest one (Shape.STArea() is the name for volume in my database) of each city name (NAME is the attribute for city name).

I don't have any formal programing experience, and got as far as creating a query to only label the single largest attribute by doing:

Shape.STArea() =(SELECT MAX( Shape.STArea() ) FROM DBO.CityArea)

However I couldn't figure out how to do this for each "NAME" I was tinkering with loops but I couldn't seem to get anything to work.

Note: none of the cities have the same name other than the fragments, and unfortunately I can't export or change the data, as it is not mine, and I am making layer files for our DPW who all pull from this database.


texts: I am trying to make a map in QGIS composer.  The final image needs to be transparent and I cannot get rid of the unwanted lines by making the bounding box smaller.  Is there a way to remove/erase lines that I do not want in the final map (circled in red)?


texts: Geoserver GetLegendGraphic not working for single rule

texts:
<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
    xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
    xmlns="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>Parcel</Name>
    <UserStyle>
      <Title>Parcels</Title>
      <FeatureTypeStyle>
        <Rule>
          <MinScaleDenominator>100</MinScaleDenominator>
          <MaxScaleDenominator>6000</MaxScaleDenominator>
          <PointSymbolizer>
             <Graphic>
               <Mark>
                 <WellKnownName>circle</WellKnownName>
                 <Fill>
                   <CssParameter name="fill">#CC3300</CssParameter>
                 </Fill>
               </Mark>
               <Size>6</Size>
             </Graphic>
           </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

This is my style file. I want to get legend graphic wit a query http://localhost:8080/geoserver/PGIS/wms?REQUEST=GetLegendGraphic

This is returning a circle but not include text:

Bu it should be like this:

If I add multiple rules, text label is appearing.


texts: Exporting multiple municipality Quarter Sections

texts:
I am very new to python and am attempting to select a municipalities pipe data by their QS number(in attribute table), dissolve the feature class on objected and qs# then export to a .gdb. I can do this one at a time but it requires me to reenter the QS number in parts of the script. I would also like the exported feature class to have a name like "QS_01_27", "QS_1-31" etc depending on the qs name. So I guess in short I would like to:

Select qs from dataset
dissolve qs on objectID, qs#
export to a .gdb with a unique name

This needs to be repeated for 689 quarter sections. below is what I have been stuck on and what I have started.

#Import system modules
import arcpy
import os

#Set workspace and mxd
workspace = "R:/desktop/this.gdb
mxd = "R:/desktop/this.mxd"

#set Variables
fc = "ProjectSewerMains"
outfc = "QS01_27"

#Select layer by QS
arcpy.Select_analysis("fc","outfc", """ QSNO = '01-27' """)

#Dissolve
arcpy.Dissolve_management(outfc, " ",     ["OBJECTID", "QSNO"], " ", "MULTI_PART", "DISSOLVE_LINES")


texts: Determine which points are contained inside a given polygon

texts:
I have these database tables defined:

CREATE TABLE regions (
  id SERIAL PRIMARY KEY,
  name VARCHAR(150) NOT NULL
)

CREATE TABLE region_points (
  id SERIAL PRIMARY KEY,
  region_id SERIAL REFERENCES regions (id) ON DELETE CASCADE,
  index INTEGER NOT NULL,
  point GEOMETRY(Point, 4326) NOT NULL
)

CREATE TABLE stations (
  id SERIAL PRIMARY KEY,
  name VARCHAR(150) NOT NULL,
  point GEOMETRY(Point, 4326) NOT NULL
)

A region is a polygon, its points being defined in the region_points table. A station may reside at a single point somewhere inside a region polygon.

I'd like to write a single query that gets all stations inside a particular region polygon. At the moment, I don't know how, so I do it in multiple steps:

# 1. Get all points that comprise the region polygon.
SELECT point
FROM region_points
WHERE region_id = 1234
ORDER BY index

# 2. Programmatically (non-SQL) construct a linestring from the points. Not shown here.

# 3. Get all stations inside the polygon.
SELECT id, name, point
FROM (
    SELECT ST_GeomFromText($1, 4326) AS polygon
) AS query1, (
    SELECT id, name, point
    FROM stations
) AS query2
WHERE ST_Contains(query1.polygon, query2.point)

As you can see, I have to construct the linestring programmatically and then pass it into the query as an argument to ST_GeomFromText, because I don't know how to do it in SQL.

I don't want to do anything programmatically. How can I re-write all of this as one SQL query?


texts:
I am working with Google Maps. Specifically, I am working with http://hpneo.github.io/gmaps for a gps tracker simulation, here I need to update the position on each point.
At the gmaps script, markers are stored on an array, and then I was thinking of taking each point and using the setPoint method to update the position, however, the gmaps script doesn't have implemented this method. I was thinking on implement this method, but my question is then:
After having an array of markers on the gmaps script, how could I identify each one in order to update position on the correct marker.
Or probably I must store it on an external array, an associative array which help me to identifies each marker, but I am thinking that when I update them, the array on the gmap script also will be keeping at the same position without been updated.

I am attaching my code.

thanks

/** Positions and map statuses **/
var isLoaded = false;
/** **/
var refreshIntervalId;
var vehicles;
var map;
function setVehicleAsCenter (registration) {
    alert("Hola");
}
function loadPoints (positions) {
    console.debug('loading markers');
    var lttd;
    var lgtd;
    for (var i=0; i < positions.length; i++) {
        lttd = positions[i].latitude;
        lgtd = positions[i].longitude;
        marker = map.addMarker({
            lat: lttd,
            lng: lgtd,
        });
        markers[positions[i].registration] = marker;
    };
    map.fitZoom();
    isLoaded = true
}
function updatePoints (positions) {
    /**
    how could be the alorithm here
    */
    console.debug('updating markers');
}
function requestPoints() {
    $.ajax({
        url:'{% url 'gpstracking.ajax.request_tracks' %}',
        type: 'get',
        dataType: 'json',
        data: {
            vehicles: vehicles
        },
        success: function (positions) {
            if (isLoaded == false) {
                loadPoints (positions);
            } else {
                updatePoints (positions);
            }
        }
    });
}
$(document).ready(function() {
    /** Buttons Click Event Set **/
    $('.map-mode').click(function(){
        vehicles = '';
        $("#jstree").jstree("get_checked",null,true).find('a[rel="vehicle"]').each(function(){
            vehicles =  vehicles + $.trim(this.text) + "|";
        });
        if (vehicles == '') {
            console.debug('No vehicles to display');
            return;
        }
        option = $(this).attr('rel');
        if (option == 'show') {
            console.debug('Ordering to show');
            clearInterval(refreshIntervalId);
            requestPoints();
        }
        if (option == 'listen') {
            console.debug('Listening');
            requestPoints();
            refreshIntervalId = setInterval("requestPoints()", 10000);
        }
        if (option == 'clear') {
            console.debug('Clearing');
            clearInterval(refreshIntervalId);
            markers = new Object();
            map.removeMarkers();
            isLoaded = false;
        }
    });
    /** Map loading **/
    map = new GMaps({
        div: '#map-canvas',
        lat: -16.4237766667,
        lng: -71.54262,
    });
});


texts: How can I calculate degree of distortion on an Azimuthal Equidistant projection?

texts:
The Tissot's indicatrix of an Azimuthal Equidistant projection shows ellipses being stretched along latitudes.

How can I calculate the degree (amount) of this distortion?  That is, what is the ratio of major axis to minor axis of a Tissot ellipse as a function of latitude (assuming the center of projection is the north pole)?

Edit:
Found numerical values on the GIS help page.

Degrees from center     Percentage of scale distortion along parallels
15                      1.2
30                      4.7
45                      11.1
60                      20.9
90                      57

I can interpolate this well and get what I need.  But would still be good to know what the analytical formula is.  Feels like a function of cos(latitude) since an Azimuthal Equidistant projects latitudes with linearly increasing distance from center instead of cosine.


texts: Popups in Leaflet through Geoserver wms active layers only

texts:
I have popups working well in Leaflet as pulled from Geoserver utilizing getfeatureinfo and drawing layers using wms. Unfortunately I cannot figure out how to only request a popup information for layers currently drawn on the map. Right now I request my getfeatureinfo with this code:

map.addEventListener('click', onMapClick);
popup = new L.Popup({maxWidth: 1000});
function onMapClick(e) {
var latlngStr = '(' + e.latlng.lat.toFixed(3) + ', ' + e.latlng.lng.toFixed(3) + ')';
var BBOX = map.getBounds()._southWest.lng+","+map.getBounds()._southWest.lat+","+map.getBounds()._northEast.lng+","
+map.getBounds()._northEast.lat;
var WIDTH= map.getSize().x;
var HEIGHT = map.getSize().y;
var X = map.layerPointToContainerPoint(e.layerPoint).x;
var Y = map.layerPointToContainerPoint(e.layerPoint).y;
var URL = 'http://xx.xx.xx.xx/geoserver/workspace/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&LAYERS=layer1,layer2&QUERY_LAYERS=layer1,layer2,layer3&STYLES=&BBOX='+BBOX+'&FEATURE_COUNT=5&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&FORMAT=image%2Fpng&INFO_FORMAT=text%2fhtml&SRS=EPSG%3A4326&X='+X+'&Y='+Y;
popup.setLatLng(e.latlng);
popup.setContent("<iframe src='"+URL+"' width='450' height='500' frameborder='0'></iframe>");
map.openPopup(popup);
}

How can I send the getfeatureinfo request url with only the layers checked as active in my TOC in leaflet? (i.e. only layer 1). Right now I get a popup with information from and of the three layers that have a feature where clicked whether they a selected and drawn or not.


texts: Fastest way to copy feature dataset from ArcSDE to File Geodatabase in ModelBuilder?

texts:
What geoprocessing tool is the fastest for copying feature datasets from ArcSDE to a file geodatabase?

I want to make a model in ModelBuilder that copies an entire featuredataset from my SDE database to a file geodatabase. I don't know what the fastest tool is. I tried "Copy element" but that is painstakingly slow compared to if you copy-pasting in ArcCatalog. Featureclass to Geodatabase, only copies the featureclasses in the dataset, and places them in the root of the fgdb. Featureclass to featureclass is a lot of drag and dropping, and I don't know about the speed.


texts: How to adjust the tolerance when showing vertices with the V key (ArcGIS)

texts:
When in the editing mode in ArcMap, it is possible to approach the mouse cursor close enough to a feature and press the V key. This will show the vertices of the feature.

However, you need to move the cursor quite close (I believe the default tolerance is 10 pixels) to be able to show the vertices. I want to be able to quickly see vertices of features without moving my mouse cursor too much.

Is it possible to adjust this tolerance value? Editing the snapping tolerance doesn't affect the "vertex show tolerance".


texts: R geographically weighted regression GWModel

texts:
I am exploring the use of GWmodel to run some GWR regressions. I set it up and tried to run some sample regressions using the gwr.basic function, but ran into the following error,

  Error in t(X * W.i) : dims [product 4659] do not match the length of object [0]

This is what I did,

gwr.basic(Y~X+I(X^2), data = spdf,bw = bw.gwr.1, kernel = "gausssian", adaptive = TRUE, F123.test = TRUE)

And here is the summary of my input data set, spdf

**summary(spdf)
    Object of class SpatialPointsDataFrame
    Coordinates:
              min       max
    Long 95.71586 103.70657
    Lat  40.12491  42.89963
    Is projected: NA
    proj4string : [NA]
    Number of points: 1553
    Data attributes:
         Y             X
     Min.   : 51.33   Min.   :100
     1st Qu.:196.19   1st Qu.:115
     Median :224.74   Median :123
     Mean   :216.66   Mean   :123
     3rd Qu.:245.70   3rd Qu.:131
     Max.   :346.72   Max.   :150**

What may be happening?


texts: Affine Tranformation Matrix; shapely asks 6 coefficients, but rasterio delivers 9; how to translate

texts:
Some lines have been extracted via scikit-image, now tho goal is to vectorize them within the origin coordinate reference system
#some geotif opened with rasterio
src = rasterio.open(raster_path)
transform = src.tranform
print("transform",transform)
>> | 0.10, 0.00, 350482.83|
| 0.00,-0.10, 5617200.95|
| 0.00, 0.00, 1.00|
print("type(transform)",type(transform))
>> <class 'affine.Affine'>
#some line extract from the geotiff with scikit-image
line = ((849, 2363), (1626, 1711))
#extract start- and end-pointof line
p0, p1 = line
# construct points in shapely
A = Point(p0[0],p0[1])
B = Point(p1[0],p1[1])
AB = LineString([A, B])
print("AB",AB)
>> LINESTRING (849 2363, 1626 1711)
AB = affine_transform(AB, transform)
>>ValueError: 'matrix' expects either 6 or 12 coefficients

How does one get a 6 coefficient Matrix or a 12 coefficient Matrix from from a geotiff with rasterio?
or
How does one translate an 9 coefficient affine transformation to a 6 coefficient one?


texts: How to toggle custom geolocation marker in OpenLayers 3?

texts:
I've managed to to add a geolocation marker to an OpenLayers 3 map that reflects the actual position and accuracy of geolocation.

I got a custom button on the map to start the geolocation, but how can I stop geolocation and remove the marker from the map? That means use the same button to toggle my function on and off.

Here my code until now:

// new Position Control

    var button = document.createElement('button');
    button.innerHTML = 'P';

    var handleGetPosition = function(e) {
        getPosition();
    };

    button.addEventListener('click', handleGetPosition, false);
    button.addEventListener('touchstart', handleGetPosition, false);

    var element = document.createElement('div');
    element.className = 'get-position ol-unselectable ol-control';
    element.appendChild(button);

    var GetPositionControl = new ol.control.Control({
        element: element
    });

    map.addControl(GetPositionControl);

// set up the geolocation api to track our position
function getPosition() {
      var geolocation = new ol.Geolocation({
        projection: map.getView().getProjection(),
        tracking: true,
        trackingOptions: {
          enableHighAccuracy: true,
          maximumAge: 5000
        }
      });

     var accuracyFeature = new ol.Feature();
      geolocation.on('change:accuracyGeometry', function() {
        accuracyFeature.setGeometry(geolocation.getAccuracyGeometry());
      });

      var positionFeature = new ol.Feature();
      positionFeature.setStyle(new ol.style.Style({
        image: new ol.style.Circle({
          radius: 4,
          fill: new ol.style.Fill({
            color: '#3399CC'
          }),
          stroke: new ol.style.Stroke({
            color: '#fff',
            width: 1
          })
        })
      }));

      geolocation.on('change:position', function() {
        var pos = geolocation.getPosition();
        positionFeature.setGeometry(pos ?
            new ol.geom.Point(pos) : null);
        view.setCenter(pos);
        view.setZoom(19);
      });

      new ol.layer.Vector({
        map: map,
        source: new ol.source.Vector({
          features: [accuracyFeature, positionFeature]
        })
      });
};


texts: FME Server Automation: Watch URL/Published data in the internet

texts:
I'd like to build an automation in FME Server consisting of one trigger (watch URL/xlsx-sheet) and an email-action.
The problem is about the trigger: There is a federal environment agency that publishes an excel-file with spatial data every few months. Every time the file gets an update (or they publish a new file) I'd like to recive an email.
The URL I'd like to watch is of this syntax:
https://www.environment-agency.de/fileadmin/documents/abc/xyz/DataSheet_01.xlsx
Is this possible? I tried the directory watch trigger, but it says, the path in unvalid.


texts: Opening shapefile linked custom metadata textfile within QGIS

texts:
I used a text file as a simple metadata for my shapefiles. It has basic information such as Source of raw data, Who can edit the file, Accuracy level, what those dbf filed means, Cost aquiring data, etc.
The text file is saved with the same name and path with my Shapefile e.g.

Conservation Area.shp
Conservation area.txt     ----> the file

How do I open this text file when its shapefile is an active layer?
I know, the perfect solution would be to use standard metadata format (shp.xml) or learn to write a plugin.
Temporarily, I'm doing it as follows.
Create field named "meta" in the shapefile and populate EVERY row of the field with path and filename of the textfile.Then Create action > type : "Windows", action: notepad "[% "meta" %]
To View,  I select the new action Icon and click on any vector.
I am using win 7, x64, QGIS 2.10.


texts: Can I use "Web Map Service Module" and "Web Feature Service Module" of Geoserver in my project?

texts:
I'm developing a GIS application with java Spring Framework and GeoTools library. I have some vector layer. I want to get this layers with wms and wfs protocols.
I've downloaded geoserver code. It includes module as bellow:

Main Module
Core Platform Module
Web Map Service Module
Web Feature Service Module
and etc.

Can I use Web Map Service Module and Web Feature Service Module of Geoserver in my project?
If yes, How?
I want a sample code to use these modules apart from geoserver. Help me please.


texts: Arcpy - Convert point in WSG84 to EPSG:2056

texts:
Currently, i have some points in WSG84 coordinate and i want to convert them to EPSG:2056 using this snippet:

coordList = [[46.34188, 8.02521], [46.34191, 8.02525], [46.34195,8.02516], [46.34188, 8.02521]]
polygonArray = arcpy.Array()
for p in coordList:
  point = arcpy.Point(p[1], p[0])
  ptgeo = arcpy.PointGeometry(point, in_spatial)
  ptgeo1 = ptgeo.projectAs(out_spatial)
  pt1 = ptgeo1.lastPoint
  polygonArray.add(pt1)

This snippet will produce this result:

2645081.068384 - 1132306.244629
2645084.122581 - 1132309.602093
2645077.161386 - 1132313.996292
2645081.068384 - 1132306.244629

But if i use this website: https://mygeodata.cloud/cs2cs/ to convert my coordinate it give me different result:

2645159.23603 - 1132445.21286
2645162.29062 - 1132448.57073
2645155.32853 - 1132452.96545
2645159.23603 - 1132445.21286

And i realize that the result from the website above is the correct one. So what is the problem with my code? Why there is the different between my code and the result from mygeodata?


texts:
I have this code in which I would like to perform a for loop in which every time I get one raster as input (the one called L, different for every month) and divide it always by the same raster (the one called U).
The code seems to work, but it always gives me as result an invalid raster since all the pixels have value 3.40282e+38, BUT if I do the same manually with the Raster Calculator (so not with a python script) it gives me the right result.
Can someone tell me where is the mistake in my code?
import qgis
import gdal
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

U90 = QgsRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/rasterUfinali/U_1990.tif', 'U90')

entries = []

u = QgsRasterCalculatorEntry()
u.ref = 'U@1'
u.raster = U90
u.bandNumber = 1

entries.append( u )

inputpath= r'/Users/macbook/Desktop/TESI/PROGETTO/L/'

suffix_input=['Ljan', 'Lfeb', 'Lmar', 'Lapr', 'Lmay', 'Ljun', 'Ljul', 'Laug', 'Lsep', 'Loct', 'Lnov', 'Ldec']

outputpath= r'/Users/macbook/Desktop/TESI/PROGETTO/lamb/'

suffix_output=['lambdajan','lambdafeb','lambdamar','lambdaapr','lambdamay','lambdajun','lambdajul','lambdaaug','lambdasep','lambdaoct','lambdanov','lambdadec']

layers_path = []

layers = [ QgsRasterLayer(inputpath + suffix_input[i] + '.tif', suffix_input[i] ) for i in range(len(suffix_input)) ]

raster_ref = [ '{}@1'.format(suffix) for suffix in suffix_input ]

for i in range(12):

    layer = QgsRasterCalculatorEntry()
    layer.ref = raster_ref[i]
    layer.raster = layers[i]
    layer.bandNumber = 1
    entries.append( layer )

    new_path = outputpath + suffix_output[i] + '.tif'

    l = layer.ref  +'/'+ u.ref
    print(l)

    lamb = QgsRasterCalculator(l, new_path, 'GTiff', layers[i].extent(), layers[i].width(), layers[i].height(), entries )

    lamb.processCalculation()
    iface.addRasterLayer(new_path)


texts: GRASS plugin missing from QGIS 2.14.3 Essen in Ubuntu 14.04

texts:
I had installed QGIS 2.14.0 (or .1 not sure) and I did some work with the GRASS plugin a couple of months ago, so it was there.

I upgraded a couple of weeks ago to QGIS 2.14.3 and I was planning in using the GRASS plugin is now missing.

GRASS is installed and works from the Processing tools, but the plugin is not to be found anywhere... (and yes the I have installed the plugin with sudo apt-get install qgis-plugin-grass)

I have been looking for a solution and none of the ones I found seem to work. The closest one was this one but none the answers changed the situation.

(System: GRASS70, QGIS 2.14.3 Essen in Ubuntu 14.04)

Update: I found out that the grass-plugins that I have in my machine are "libgrassplugin6.so, libgrassprovider6.so and libgrassrasterprovider6.so". So the problem seems to be that I have GRASS 7 installed but the plugin is for GRASS 6, right?

Update2: I found this question in the users list (sorry I did not look earlier) and seems to be the same/similar problem... from the answers there you don't quite find the solution but I am working on it...


texts: How to save a a CartoDB buffer to its own data table/layer

texts:
This seems like a dumb question, but I want to create (or add to) test_table2 a buffer feature from test_table1 point layer.
INSERT INTO test_table2
(SELECT ST_Buffer(the_geom_webmercator, 1000) AS the_geom_webmercator, cartodb_id FROM test_table1)

In this example, test_table2 already exists and shares the same columns as test_table1
Update
In the CartoDB query window online in the Data View/Map View, this statement works:
SELECT ST_Buffer(the_geom_webmercator,500) AS the_geom_webmercator, cartodb_id FROM test_table1

So I basically need to INSERT this into a new or already created table.
Update2
So I got the following statement to work:
INSERT INTO testcreate (the_geom) VALUES (ST_Buffer(CDB_LatLng(40.7127, -74.0059),0.01));

I don't think this really solved my problem, unless I can get ST_Intersects, and ST_Transform (as I need this in a projected system - the result of this is unprojected) to fit inside the VALUES() function.
Advice provided from @schwanksta on twitter:
schwanksta: @nygeog Hm, maybe misunderstanding, but can't you UPDATE t1 SET buffer = (SELECT ST_Buffer(t2.geom,100) FROM t2 WHERE t1.id=t2.id)
@nygeog t1 would be the table you're trying to stuff the value into, t2 would be the source pt layer. sorry, reversed your naming scheme!
@nygeog You'd need to create a column with AddGeometryColumn first
Tried this:
UPDATE t2 SET buffer = (SELECT ST_Buffer(t1.the_geom,100) FROM t1 WHERE t1.cartodb_id=t2.cartodb_id)

This didn't seem to work for me, I could be translating to meet my data wrong (I did create two tables t1 and t2, but getting urllib2.HTTPError: HTTP Error 400: Bad Request from this.


texts: How to extract a new linestring vector from a existing multipolygon in QGIS?

texts:
What I want to do: create a new line (named linestring in Postgis?) from existing multipolygon layer in Postgis.

Purpose: once these lines are created, I'll run some distance calcs from specific points.

What I'm using: Postgres 9.3 with Postgis enabled. Quantim GIS 2.2 Valmiera.

The existing multipolygon layer are neighborhoods in Rio de Janeiro city. The lines I want to create are clean beaches that some neighborhoods have.

What I'm doing: Like in the image below, I select the layer in QGIS, click on the button "Toggle Editing", click on "Node Tool", select the neighborhood to show its points with the Node Tool, select the points based unpon a rectangle I draw, the poinst become highlighted in blue, and then I don't know how to save them as a new line.

UPDATE: I've found this tutorial to manually create a new feature. It will be my backup technique if there is no other.


texts: OpenLayers 3 - Multi feature vector layer - how to select all features

texts:
I have multiple vector layers which is constructed of 3 parts: start point, line, end point. Each of those features has its own style. start point is a green circle, line is grey and end point is a red circle.

I also enabled selection and it is working, but how do I programmatically select the 2 other features when selecting one of the parts?

My code is below:

var featureLine = new ol.Feature({
        geometry: new ol.geom.LineString(points),
        type: "line",
        id: id1,
        startdate: startdate1,
        starttime: starttime1,
        enddate: enddate1,
        endtime: endtime1,
        soaktime: soaktime1,
        haulnumber: haulnumber1,
        area: area1,
        startlat: startlat1,
        startlon: startlon1,
        endlat: endlat1,
        endlon: endlon1
    });
    var featureStart = new ol.Feature({
        geometry: new ol.geom.Point(ol.proj.transform([startlon1, startlat1], 'EPSG:4326', 'EPSG:3857')),
        type: "start",
        id: id1,
        startdate: startdate1,
        starttime: starttime1,
        enddate: enddate1,
        endtime: endtime1,
        soaktime: soaktime1,
        haulnumber: haulnumber1,
        area: area1,
        startlat: startlat1,
        startlon: startlon1,
        endlat: endlat1,
        endlon: endlon1

    });
    var featureEnd = new ol.Feature({
        geometry: new ol.geom.Point(ol.proj.transform([endlon1, endlat1], 'EPSG:4326', 'EPSG:3857')),
        type: "end",
        id: id1,
        startdate: startdate1,
        starttime: starttime1,
        enddate: enddate1,
        endtime: endtime1,
        soaktime: soaktime1,
        haulnumber: haulnumber1,
        area: area1,
        startlat: startlat1,
        startlon: startlon1,
        endlat: endlat1,
        endlon: endlon1
    });

    var vectorSource = new ol.source.Vector({
        projection: 'EPSG:4326',
        features: [featureStart, featureEnd, featureLine]
    });

    var styleStart = [
        new ol.style.Style({
            image: new ol.style.Circle({
                radius: 7,
                fill: new ol.style.Fill({
                    color: '#32CD32'
                })
            })
        })
    ];

    var styleLine = [
        new ol.style.Style({
            fill: new ol.style.Fill({
                color: 'rgba(255, 255, 255, 0.2)',
                weight: 4
            }),
            stroke: new ol.style.Stroke({
                color: '#808080',
                width: 4
            })
        })
    ];

    var styleEnd = [
        new ol.style.Style({
            image: new ol.style.Circle({
                radius: 7,
                fill: new ol.style.Fill({
                    color: '#FF4500',
                    opacity: 0.8
                })
            })
        })
    ];

    var styles = [styleStart, styleLine, styleEnd];

    var vectorLayer = new ol.layer.Vector({
        source: vectorSource,
        style: function (feature, resolution) {
            if (feature.get('type') === "line") {
                return styleLine;
            }
            if (feature.get('type') === "start") {
                return styleStart;
            }
            if (feature.get('type') === "end") {
                return styleEnd;
            }

            return styles;
        }
    });

    map.addLayer(vectorLayer);

And my select event on initialization is as follows:

        var selectClick = new ol.interaction.Select({
        condition: ol.events.condition.click
    });
    map.addInteraction(selectClick);


texts:
There are several answers for this question on GIS stackexchange but none using open/free software.

How can someone without access to the Esri suite remove interior holes from a polygon?

Example data is available here


texts: Removing buffers of overlapping buffers using ArcMap

texts:
I work with ArcMap 10.2 and would like to remove buffers of the overlapping buffers, but not all overlapping buffers. So for example:

I would like to remove the buffer in the middle and remain the outer two. I have many buffers though which overlap, I would do it for all of them, so I do not have to throw away all the overlapping buffers. I would like to use the ModelBuilder because I do not know much about Python.


texts: Difference (in decimals) in the area before and after union

texts:
I'm using ArcGIS 10.4.1

DATA

I have two layers

5_catchments: polygon shapefile for 5 catchments
nloss_clip_final.shp: polygon shapefile for Nloss

You can find and download the data from this link.

I want to get the attributes of the two layers into one layer. I can do this using intersect or union. I used union to do so.

After using union, I exported the attribute table as .txt file and imported it into R to check the total area for the 5 catchments. I found difference in decimals of the total area in km2

Area (in km2) of 5 catchments before union

Area (in km2) of 5 catchments after union

Although the difference in the area (in km2) of the 5 catchments, before and after union, is in decimals (mainly in the most downstream catchment Rangitikei_McKelvies where the difference is around 0.047 km2), I wonder why there is a difference? Is it a rounding issue? or something else.

UPDATE

All the input and output layers are in NZGD2000 (New Zealand Transverse Mercator 2000). Resolution and tolerance values are below


texts: Is it possible to merge multiple polylines into one

texts:
I have vehicle GPS coordinates of the route. In my case vehicles can drive through one road multiple times, so overall it can look like round trips. Is there any way to merge two traces into one, so at the end there will be only one red line on the map, using tools like PostGIS, OSRM or something else?

Also, currently I don't see a way to filter traces programmatically before rendering them on map, by checking if trace already present near the road.

I'm new to this site and unsure how to classify this problem.

After some searching, looks like PostGIS function ST_RemoveRepeatedPoints can help to remove close coordinates, it has parameter tolerance so the coordinates within the tolerance value will be considered as duplicates.


texts: How to save IDocumentInfo attributes to a map document programmatically?

texts:
We have been successfully creating mxd files based on MapDocument, adding layers, setting spatial information and saving them.  Our code has been working, literally, for years and has made the transition from 9.x to 10.x without trouble.

The other day we decided that we wanted to incorporate field based links that required us to use the HyperlinkBase attribute of the DocumentInfo class.  We have tried the 5 or 6 most obvious ways to do that and have not been able to accomplish it.  No matter what we insert into the IDocumentInfo or IDocumentInfo2 interface, it does not get saved with the mxd.

Furthermore, if you open an mxd that already has attributes (that you added with ArcMap, for example) twiddle them, and then save the mxd again ... all the attributes are cleared. This is intensely frustrating of course. Our code is of this form ...

IMapDocument mapDoc = null;

try
{
    mapDoc = ArcObjects.CartoUtil.OpenMxd(mxdPath, null);

    AddMapLayersToMxd( "our hidden, sanitized stuff goes here" );

    IDocumentInfo2 infoDoc = (IDocumentInfo2) mapDoc;
    infoDoc.HyperlinkBase = BlahBlah.ServerBaseUrl;

    mapDoc.Save(true, false);
    mapDoc.Close();
}
finally
{
    ComUtil.Release(mapDoc);
}

Where AddMapLayersToMxd() does just what it says. This works fine and the mxd has the layers and other things we add to it.  Our desired DocumentInfo ... nowhere to be found.

If you crack the resulting mxd, you will find only the DocumentInfo added by the system ... paths, dates, RelativePaths, etc.  None of the user-defined items.


texts:
I have an equidistant point grid (0.5m) and created a raster from it in QGIS and ArcMap, with different results:

QGIS:

Raster -> Conversion -> Rasterize
Result: raster, with points in the
top left corner of each "pixel"

ArcMap:

Toolbox: Point To Raster
Result: raster, with points centered

Which is the better approach? Can I emulate ArcMap's behavior with QGIS?


texts: How do I install Orfeo toolbox and Qgis 1.9?

texts:
I'm having a lot of trouble getting Orfeo toolbox and Qgis 2.0 working.
I followed the qgis wiki instructios, added the ppa, installed qgis 1.9 (it's beautiful, great work!).
But when i tried to install the otblib from the ubuntugis repositories, it forces an uninstall of Qgis.
Any ideas ?

My system is linux mint 15.


texts: Facing Issue on Creating MVVM Pattern From An ESRI Runtime 100.2 Existing App

texts:
Following This Sample I am hoping to be able a real MVVM pattern from the tutorial but based on my understanding the application is missing the Model and View classes!

I have the MapViewModel.cs like this

 public class MapViewModel
    {
        public MapViewModel(){ }

        private Map _map = new Map(Basemap.CreateStreets());
        public Map Map
        {
            get { return _map; }
            set { _map = value;  }
        }

    }

and the MainWindow.xaml

<Window.Resources>
    <local:MapViewModel x:Key="MapViewModel" />
</Window.Resources>
<Grid>
    <esri:MapView Map="{Binding Map, Source={StaticResource MapViewModel}}" />
</Grid>

but whre are the "MapView and "MapModel classes? Can you please help me to extracts and create those classes from the MapViewModel and create a real MVVM model?


texts: I installed CartoDB on a local machine and it worked fine. However, when I tried to go offline I was able to login, but I wasn't able to do anything. It looks like the page is looking for online data. Is it possible to completely work with CartoDB in an offline mode?


texts: Getting table Name from IObjectClass using ArcObjects?

texts:
I'm using IRelationshipClass.DestinationClass to get an IObjectClass.  From here, I want to be able to get the table name, but can only access the AliasName.

Any ideas?

The IObjectClass can be cast to a IDataset, which has the Name property, which is not the same as AliasName.  I should say thanks to Kirk, as it was some code he posted elsewhere that pointed me in the right direction.

For reference,

((IDataset)relClass.DestinationClass).Name


texts:
Following the guidance at,
http://docs.qgis.org/html/en/docs/user_manual/working_with_gps/plugins_gps.html
I checked my Plugins Manager and verified that GPS Tools is enabled. However, unlike what the reference site indicated, I only have one button on my tool bar--the one for "Create New GPX Layer"; not the one for "GPS Tools", which is used for loading GPS data from a GPS unit.

Any idea why I'm missing that button? Or how I should load data from my GPS unit? (I'm using QGIS 1.8.0 and Windows 7.)

Thanks
TVZ


texts: Raster output symbology differs in tool vs. Python

texts:
I have a bunch old DOQQs from the 90's that are projected in NAD27, giving them the black angled collar you see in the image. I am trying to set the collar to null by running SetNull on Value = 0. If I run the Set Null tool from the SA toolbox it turns out perfectly. If I run SetNull in the python window however, it changes the symbology from stretched to unique values, requiring me to manually change it back on every DOQQ. Is there a way to keep the symbology set to stretched in the save output, or to change it back in Python?

from arcpy.sa import *
arcpy.env.workspace = r'F:\DOQQ1991\Tiles'
rastiles = arcpy.ListRasters()
for tile in rastiles:
    nullras = SetNull(tile, tile, "VALUE = 0")
    nullras.save("{0}\{1}".format(r'F:\DOQQ1991\NullTiles', tile))


texts: Lines length plugin - QgsVectorLayer error

texts:
Several similar questions have already been made, but none of the is helping me.

I have a layer of polygons and I want to known the length of each line of each polygon. And based on the length I want to create another layer.
So I convert polygons > line and apply the explode lines tool. But then when I try to get the length, this is the error: AttributeError: QgsVectorLayer object has no attribute geometry.

My code is this:

def run(self):
   # show the dialog
    self.dlg.show()
    self.dlg.layerCombo.clear()
    layers = QgsMapLayerRegistry.instance().mapLayers().values()

    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.layerCombo.addItem( layer.name(), layer )
    result = self.dlg.exec_()
    # See if OK was pressed
    if result == 1:
        index = self.dlg.layerCombo.currentIndex()
        layer = self.dlg.layerCombo.itemData(index)
        #polygons to lines
        processing.runalg ('saga:convertpolygonstolines',layer,layer_line)

I know I need to use somenthing like this to have a new field with length values:

for f in features:
        geom = f.geometry()
        leng = geom.length()
        res = layer.dataProvider().addAttributes([QgsField("Length", QVariant.Int)])
        layer.updateFields()
        fieldIndex = layer.dataProvider().fieldNameIndex( "Length" )

I know that my problem is using QgsMapLayer.VectorLayer...But I don't know
whatelse to use. I already read PyQGIS Developer Cookbook, the answer to this question (Calculate length of selected features in PyQGIS), this question (Adding new field with expression in pyqgis)  ...

What should I use instead of QgsVectorLayer?

Ok. I made the changes suggested by Michael.
My code is now like this:

            res = layer.dataProvider().addAttributes([QgsField("Length", QVariant.Int)])
            layer_output3.updateFields()
            fieldIndex = layer.dataProvider().fieldNameIndex( "Length" )
            features = layer.getFeatures()
            layer.startEditing()
            for f in features:
                geom = f.geometry()
                leng = geom.length()
                layer.changeAttributeValue(f.id())
                layer.commitChanges()

I am missing something..

So I corrected the itens mentioned by Michael and David. The part where the new field is add and the length is calculated looks like this:

            features = layer.getFeatures()
            coluna = QgsField('Length', QVariant.Double)
            res = layer.dataProvider().addAttributes([coluna])
            leng = layer.dataProvider().fieldNameIndex
            layer.updateFields()
            fieldIndex = layer.dataProvider().fieldNameIndex( "Length" )
            layer.startEditing()
            for f in features:
                geom = f.geometry()
                leng = geom.length()
                layer.changeAttributeValue(f.id(),fieldIndex, leng)
            layer.commitChanges()
            layer.updateFields()


texts:
does anyone know how I can load local MBTiles from the device storage using the MapBox SDK?

The vector Tiles come from this Website: https://openmaptiles.com/downloads/planet/

The same page offers Android Source Code which contains a solution for my problem but I don't want to pay over 500$.

https://openmaptiles.com/mobile-app/

Quote: "We have extended the Mapbox SDKs (C++) with the ability to load the map tiles from local MBTiles files that are either bundled with the app or downloaded by the user. "

I don't want to use a local web server for this.


texts: OpenLayers 4 change position on load and input change

texts:
I have map which should show position depends on input, which I control with switch. Long and Lat I am getting from API (depends on my input) which return to me coordinates. How can I show position from my coordinates? For me all time shows same position. Here is my code: Where am I wrong?

if(smth !='' && smth!=0){

    $.ajax({
        type: "POST",
        url: 'myApiUrl',
        data: {
            'region' : data_from switch,
            'my_key': 'my_key_here'
    },
        success: function(e){
            console.log(e);
            if(e.names.length > 0){

            center = ol.proj.transform([e.names[0].X, e.names[0].Y], 'EPSG:4326', 'EPSG:3857')

            var lonlat =ol.proj.transform([e.names[0].X,e.names[0].Y], 'EPSG:4326', 'EPSG:3857');
            var feature =  new ol.Feature({
            geometry: new ol.geom.Point([Math.floor(lonlat[0]),Math.floor(lonlat[1])])
                });
                    vectorLayer.getSource().addFeature(feature);
                    map.getView().setCenter(center);
                    map.getView().setZoom(14);

                console.log(lonlat);
                loadMarkers(e.names[0].X,e.names[0].Y);
            }
        }
    });
}

Also I need to change position on change of same input, which will also use same switch and return coordinates. I have create crated function, which normally should show my position on change of input, it does not work also. Please help me to find mistake. Here is the code:

function loadMarkers(x, y){
    //console.log(x, y);
    var lonlat =ol.proj.transform([x,y], 'EPSG:4326',   'EPSG:3857');

    var feature =  new ol.Feature({
        geometry: new ol.geom.Point([Math.floor(lonlat[0]),Math.floor(lonlat[1])])
    });

    vectorLayer.getSource().addFeature(feature);
    map.getView().setCenter(lonlat, 9);
}

This function I call on input change where inserting Long as 'Y' and Lat as 'X', but it does not show me anything also, also no errors.


texts:
I have a clipped DEM of my study area in ArcGIS Desktop 9.3 which I want to use for SWAT project set-up, but I need to mask it.

How do I go about it?

I'm new to the use of ArcGIS.

I don't want to use the extract by mask option.


texts:
I have problem with CRS Bosnia and Herzegovina?

MGI 1901 / Balkans zone 6  EPSG:3908
MGI 31276 / Balkans zone 6  EPSG:31276

Which of these two is OK?
When I open google sattelite image with openlayers plugin, MGI 1901 match an image, when I change to MGI 31276 vector layer mismatch about 300m.
When I import my layers to Global Mapper, datum of vector is "MGI Austria"

I need this projection:


texts: r-Mapview non reponsive with large simple feature dataframe?

texts:
I am trying to visualize a large sf dataframe with around 1000 rows containing linestrings with mapview::mapview() in R. Unfortunately a warning comes up when I perform the command:

mapview(data.frame)

Warning message:
In leaflet_sfc(sf::st_geometry(x), map = map, zcol = zcol, color = clrs,  :
  the supplied feature layer has more points/vertices than the set threshold.
  using special rendering function, hence things may not behave as expected from a standard leaflet map,
  e.g. you will likely need to zoom in to popup-query features

  to see the number of points/vertices of the layer use 'npts(x)'
  to see the threshold for the feature type use 'mapview:::getMaxFeatures(x)'
  to adjust the threshold use argument 'maxpoints'

And then I try to export as html and I get:

R code execution error

And the saved file is 0 kb.

Is there a fix or a workaround with some other function from the package?


texts: Pre-populating list in Survey123 using SQL Server column

texts:
What I currently have is a Survey which enters the data collected from the survey in a SQL Server table, I'd like for one of those columns in that table to become the choices in another survey.

If this is not possible, are there any other tools I can use to accomplish this using ArcGIS?


texts: Display millions of points in the browser

texts:
I have a web application which will track he movement of user, and save the history location to database.

And user may want to see his/her movement trace over a period of time, we have to display them in the browser(contain IE6).

Generally we can use the overlay API like Polyline in Google Map or Leaflet JS, however it is unacceptable for our use case since there maybe too many points for a time.

Now the user will upload their location every 60 seconds, and suppose they just upload 10 hours one day, then there will be 10 * 24 * 60 = 14400 records per user one day.

So how about user trying to check the trace for yesterday or even last month? Technologically we prefer to load just a part of the data and then load more ,but user like to see All of the trace.

That's means we can not draw the points in the client, I have though to simplify the points before returned to the client, but this can not make sure the number of the points to an acceptable level, so I wonder if someone have the similar experience for this kind of situation?


texts: Making an ellipse based on an ellipse-like polygon in QGIS

texts:
I'm using QGIS. I have a layer that contains ellipse-like polygons as in the image. I need to convert them to ellipses. They come from a different data sources. These are all actually regular shapes that originate from ellipses. Therefore, I think, somehow, an ellipse containing all vertices can be formed. But an approximate ellipse also works for me.

I use Smooth tool using different parameter values. But it doesn't give expected ellipses (brown polygons in the next image). And they are not bounding ellipses. If all ellipses had the same size, the Smooth + Buffer tools combination or "outer (positive) then inner(negative) buffer" could be work.

How can I  convert them to proper (bounding) ellipses?


texts:
How can I check the page size (and maybe modify)  of the GiST index in PostGIS? From the PostGIS documentation, I quote:

  "... GiST indexes support the concept of "lossiness" which is
  important when dealing with GIS objects larger than the PostgreSQL 8K
  page size. Lossiness allows PostgreSQL to store only the "important"
  part of an object in an index -- in the case of GIS objects, just the
  bounding box. GIS objects larger than 8K will cause R-Tree indexes to
  fail in the process of being built."

Not only am I having some trouble understanding this whole concept but I am a bit confused what this implies for my original question. What does this mean that the page size is? Is it necessarily 8K?


texts: Need to export a feature dataset into a new geodatabase

texts:
I am having trouble trying to export a feature dataset into a .gdb that I have created. I am working on data in ArcMap 10.2 that is located in my d drive (D:). I have created a folder and a (for example)test.gdb in my p drive where I want my test.gdb located. However, when I right click the feature dataset in my table of contents (because that is the only thing I want to export) and click on "export," I am confused on what to do next.

For my output I choose the path of the feature dataset I want to export. My input, I choose the folder in which I want my dataset to go too (p:\test.gdb). There is another "input" section that I am not sure what to insert as the path.


texts: How to classify vegetation in 500x500 m polygons in a very coarse way in QGIS

texts:
For a simple species distribution map for the Dutch Wadden island of Terschelling, I want to give each 500x500 m polygon (reference grid = EPSG 28992) a specific color that coarsly corresponds with the underlying vegetation/habitat type.

In the pictures you can see what I mean: The first map is one that I want to use as reference (of course I can also use sattelite photos) for vegetation (the dots represent the dispersion of a species and the polgons are vegetation types). The first picture was my initial idea, but since the maps will be displayed on a small scale in a book, i reconsidered and went on to the second picture. Here i added a vector grid of 500 x 500 meters and used the TomBio addon to convert the point shapes to 500x500 points corresponding with the grid. In the second picture I also added colours for vegetation types (e.g. yellow is beach vegetation, purple is heather vegetation, green is dune meadow etc.). I did this by just eyeballing the underlying map and give each polyon the colour I thought was right.

My Question
Is there a way of automatically generating pre-determined colours to polygons based on the vegetation in that polygon. i.e. a polygon consists of 80% forest, 15% dune meadow and 5 % heather, the polygon will be coloured dark green because most of the polygon consists of forest. This seems to me a better way then just eyeballing a polygon, but I can't figure out how to do it.

Note: I really want it to be a very simple map. So no subdivisions in forest or dunemeadow types etc. The classes i want are:
-Forest (dark green)
-Beach and dune (yellow)
-urban (gray)
-Polder (red)
-heather (purple)
-dune meadow (light green)
-Salt marsh (Brown)

So more of a habitclassification based on the (in some cases absence of) vegetation.