texts: How to use ggplot with Processing in QGIS 2.0?

texts:
I'm trying to use R scripts within QGIS thanks to Processing. There are some example scripts available and I achieved to plots different columns of the attribute table (with the plot built-in package).

What I cannot do is using the ggplot package. This is one of the simple code I've written:

##Vector processing=group
##showplots
##Layer=vector
##Field1=Field Layer
##Field2=Field Layer
library("ggplot2")
ggplot(Layer)+
geom_point(aes([[Field1]], [[Field2]]))

But none plot is shown.

Log output says:

R execution commands
options("repos"="http://cran.at.r-project.org/")
tryCatch(find.package("ggplot2"), error=function(e) install.packages("ggplot2", dependencies=TRUE))
tryCatch(find.package("rgdal"), error=function(e) install.packages("rgdal", dependencies=TRUE))
tryCatch(find.package("raster"), error=function(e) install.packages("raster", dependencies=TRUE))
library("raster")
library("rgdal")
Layer = readOGR("/media/data/Universita/brema/gis/file_tesi/piezo_e_logger",layer="piezo_tutto_finale")
Field1="Br"
Field2="Ca"
png("/tmp/processing/f02d43598eb44b63aec89c02de38af54/RPLOTS.html.png")
library("ggplot2")
ggplot(Layer) + geom_point(aes(Layer[[Field1]], Layer[[Field2]]))

dev.off()

And that's a piece of the dbf table (trying to plot Cl vs Br):

ID                     Cl   Br
2818/22/0834    125.70  0.20
2818/22/0851    219.12  0.29
2818/22/0853    298.55  0.26
2818/22/0886    132.32  0.09
B0F                520.60   1.01
B12T           3435.49  3.83
B1T                2322.61  2.54
B21T           2800.90  3.18
B22F            366.87  0.64
B26T           1891.63  2.10
B33F            216.97  0.62
B4F                 383.98  0.44
P1                   82.01  0.14
P21                 365.56  0.74
P24                 199.98  0.73
P25                 308.71  0.90
P27                 217.58  0.58
P28                 436.09  0.32
P29                 711.77  1.08


texts:
I am trying to create a table in PostGIS for some LiDAR data and I need to set the the_geom column to my data's SRID.  My data is in EPSG 102743, but when I try to use 102743 as my SRID as the_geom I get this error:
ERROR:  AddGeometryColumn() - invalid SRID
spatialreference.org and epsg.io both say that it should be 102743, but that isn't working.  This is my PSQL statement:
ERROR:  AddGeometryColumn() - invalid SRID

What is the correct SRID for EPSG 102743?


texts: Difference between bounding box, envelope, extent, bounds?

texts:
Every time when I create a script to manipulate vector (point, line and polygon) in python I meet this doubt.

What is the difference between bounding box (bbox), envelope, extent, bounds?


texts:
I imported the OSM data into my Postgresql 9.5 database, I also installed the Postgis. So when I use these SQL-Syntax:

SELECT
ST_Y(ST_Transform(point.way, 4326)) AS lat,
ST_X(ST_Transform(point.way, 4326)) AS long,
poly.name as place,
point.name as location
FROM planet_osm_point point left join planet_osm_polygon poly
ON st_intersects(point.way, poly.way)
WHERE poly.admin_level='8'
AND point.sport='tennis'
AND point.name IS NOT NULL;

I get following table:

   lat        |       long       |            place            |   location

-----------------------+----------------------------+---------------+-----------------------------------
 48.1122856691474 | 11.4279876860241 | Gräfelfing                  | Tennis Grün-Weiß
48.180304681448  |  11.263177913097 | Fürstenfeldbruck            | Tennisclub Fürstenfeldbruck e.V.

So next i want to search a radius with the specified location (lat,lon):

SELECT
ST_Y(ST_Transform(point.way, 4326)) AS lat,
ST_X(ST_Transform(point.way, 4326)) AS long,
poly.name as place,
point.name as location
FROM planet_osm_point point left join planet_osm_polygon poly
ON st_intersects(point.way, poly.way)
WHERE poly.admin_level='8'
AND point.sport='tennis'
AND ST_Intersects(poly.way,ST_Transform(ST_Buffer(ST_SetSRID(ST_Point(48.1122856691474,11.4279876860241),4326)::geography,500)::geometry,3857))
AND point.name IS NOT NULL;

But now, I received no results, I also used the location of the first result (48.1122856691474 | 11.4279876860241). Why can postgres find any results?

Using a tip from @MappaGnosis, my SQL syntax now looks like this:

SELECT
name,
ST_Y(ST_Transform(way, 4326)) AS lat,
ST_X(ST_Transform(way, 4326)) AS long
FROM planet_osm_point
WHERE ST_DWITHIN(way, ST_TRANSFORM(ST_SETSRID(ST_MAKEPOINT(11.4279876860241,48.1122856691474),4326),900913), 500);

And it works! Now together with my previous code:

SELECT
ST_Y(ST_Transform(point.way, 4326)) AS lat,
ST_X(ST_Transform(point.way, 4326)) AS long,
poly.name as place,
point.name as location
FROM planet_osm_point point left join planet_osm_polygon poly
ON st_intersects(point.way, poly.way)
WHERE poly.admin_level='8'
AND point.sport='tennis'
AND ST_DWITHIN(point.way, ST_TRANSFORM(ST_SETSRID(ST_MAKEPOINT(11.4279876860241,48.1122856691474),4326),900913), 500)
AND point.name IS NOT NULL LIMIT 20;


texts: Official source for QGIS icons or banners?

texts:
I am going to use maps and analysis results generated by QGIS in a presentation. Besides citing the official reference i want to show a small QGIS icon or Banner on the first slide.
I found some nice examples on the net (1)(2), but i wonder if there is an official source providing a QGIS (and maybe GRASS) icon for such purposes ?


texts:
I have a tab delimited text file with a column "Status". It can have values OFF-AIR or ON-AIR. How do I control the marker properties like colour or out-line using the expression building. I am not looking for categorized symbol, I need the expression.
This is my syntax, but it is not working:

CASE WHEN  "STATUS"  = 'ON AIR' THEN  color_rgba(255,255,0,80)
  WHEN  "STATUS"  = 'OFF AIR' THEN  color_rgba(255,0,0,80)
END


texts: Create a grid of NDVI maps from different dates

texts:
I want to create a grid of 4 NDVI maps from different dates with the same parameters so I can check variance over a period

I tried to change an example named 'Linked Maps' in User Interface folder in order to do so, but I only get 4 empty maps without any layers
I still need to use inspector to see the NDVI value in each map (once I click a point in a map) in order to adjust the visualisation paramters for all the maps

var img1 = ee.Image('COPERNICUS/S2/20180302T081839_20180302T083044_T36SYB');
var img2 = ee.Image('COPERNICUS/S2/20180307T081801_20180307T081804_T36SYB');
var img3 = ee.Image('COPERNICUS/S2/20180312T081729_20180312T082903_T36SYB');
var img4 = ee.Image('COPERNICUS/S2/20180317T081651_20180317T082118_T36SYB');

var img1_ndvi = img1.normalizedDifference(['B8','B4']);
var img2_ndvi = img2.normalizedDifference(['B8','B4']);
var img3_ndvi = img3.normalizedDifference(['B8','B4']);
var img4_ndvi = img4.normalizedDifference(['B8','B4']);

var IMAGES = [ 'img1_ndvi',
              'img2_ndvi',
              'img3_ndvi',
              'img4_ndvi'];

var NAMES = ['Last NDVI Image',
             'Second last',
             'Third last',
             'Four last'];

var vis_params = {min: 0.1, max: 0.6, palette: 'FFFFFF, CE7E45, DF923D,     F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};

var maps = [];
IMAGES.forEach(function(image, name, index) {
  var map = ui.Map();
  map.add(ui.Label(name));
  map.addLayer(IMAGES[index], vis_params, name);
  map.setControlVisibility(false);
  maps.push(map);
});

var linker = ui.Map.Linker(maps);
// Enable zooming on the top-left map.
maps[0].setControlVisibility({zoomControl: true});
// Show the scale (e.g. '500m') on the bottom-right map.
maps[3].setControlVisibility({scaleControl: true});
// Create a title.
var title = ui.Label('NDVI maps for the Last 4 Images', {
  stretch: 'horizontal',
  textAlign: 'center',
  fontWeight: 'bold',
  fontSize: '16px'
});
// Create a grid of maps.
var mapGrid = ui.Panel(
  [
    ui.Panel([maps[0], maps[1]], null, {stretch: 'both'}),
    ui.Panel([maps[2], maps[3]], null, {stretch: 'both'})
  ],
  ui.Panel.Layout.Flow('horizontal'), {stretch: 'both'}
);
// Add the maps and title to the ui.root.
ui.root.widgets().reset([title, mapGrid]);
ui.root.setLayout(ui.Panel.Layout.Flow('vertical'));
maps[0].setCenter(35.07, 31.8, 6);


texts: Merging polygons into one single feature

texts:
I have 73 polygons of Protected Areas, and I want to merge this into one single feature of protected areas and calculate area total of all of them, I tried using the merge tool, but I keep getting the output with all 73 polygons and attribute table fields.
I also tried doing the merge with the Edit tool, but it tells me to add the merged feature into one of the protected area polygons and doesn't let me calculate geometry.
These are my protected areas

And I want to have just one field that contains them all

Any direction in how to do this?


texts: Draw Polygon using other lines in different shapefiles?

texts:
I have a scenario where i need to update a database. I have several shapefiles from which I get the data from.

I need to draw polygons inside certain borders, and the borders don't always come from the same shape file.

How can I achieve this?

We do have the Construct Polygons tool in the advanced Editing Toolbar. It does serve the purpose but it creates unnecessary polygons too.

Is there a better way to do this?

Need a polygon in the area between the three features. All these come from different shapefiles.

If i do this with Construct Polygons, I get the following:

But I need only the one inside not all the possible polygons.

Is there any tool already for this? Or If there is some way we can custom make a code for this.

Edit:
There seem to be number of ways to do this. But i am looking for something that will speed up the process as I am working on a huge database and every second saved would be helpful.

Something like this in ArcMap Would be great. Creation of polygon in the gap between two existing polygons


texts: Initialize modify interaction so it can only delete vertices, OpenLayers

texts: I want to initialize modify interaction in the way that the user can only delete vertices. I want to disable the other two functionalities which interaction offers like moving and adding new ones. Does anyone have some suggestion?


texts: How to install Python library (Openpyxl) in QGIS 2.18.16?

texts: I need to use GniPlannerFTTH, but this plugin does not work without the Openpyxl python library installed. How can I install Openpyxl in QGIS?


texts: C# implementation of concave hull / alpha shapes

texts:
I'm looking for a C# implementation of the Concave hull (alpha shapes) algorithm.
I already found several Java, Python and C++ implementation but no C# version.
The other implementation are quite complicated and will take a lot of time to convert so I'm looking for a C# version which is already working.

I'm already using NetTopologySuite and can do a delaunay triangulation. This will produce a convex hull, which I don't need.
Something like this is what I need: http://www.geosensor.net/phpws/index.php?module=pagemaster&PAGE_user_op=view_page&PAGE_id=13

Any suggestions?


texts: How can I denote the values of reclassification using variables?

texts:
I used the following lines of code to classify a raster. However, my problem is how can I denote the old and new classification values in the last code line using variables? Then number of values hence the number of variables may vary from raster to raster. I note that I know how to read these values from the raster attribute table. Thanks

' Create the Geoprocessor object
Set GP = WScript.CreateObject("esriGeoprocessing.GPDispatch.1")

' Check out any necessary licenses
GP.CheckOutExtension "spatial"

' Load required toolboxes...
GP.AddToolbox "C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx"

' Local variables...
Reclass_rast1 = "C:\LACONISS\LandSpaCES\Reclass_rast1"
Substract5_12 = "Substract5_12"

' Process: Reclassify...
GP.Reclassify_sa Substract5_12, "VALUE", "-3 0;-2 0;-1 0;0 1;1 0;2 0;3 0",  Reclass_rast1, "DATA"


texts: Stacked with "expected a raster or layer name" error

texts:
I am using the following python script snippet but I get a TypeError: expected a raster or layer name on the commented line below. I could not find a solution about how to convert the A1 raster (that is correctly created and stored on disk) to be added as addRas.

for raster in rasters:
    inputPath = os.path.join(inputFolder, raster)
    print "Processing raster: " + inputPath
    N1 = IsNull( raster )
    A1 = Con( N1, 0, raster )
    if i == 0 :
        # Got an error at the line below
        addRas = arcpy.Raster(A1)
        addNull = arcpy.Raster(N1)
        i = i + 1
        print i
    else :
        addRas = addRas + arcpy.Raster(A1)
        addNull = addNull + arcpy.Raster(N1)
        i = i + 1
        print i
avgRas = addRas / (i - addNull)


texts: How to select a number of rows based on current sort order

texts:
I would like to be able to select a certain number of rows from my attribute table based on their ordering and not their attributes. For example, if I sort by ObjectID, I then want to select the top 100 rows. If I then alphabetise the Name field, I want to be able to select the top 150 rows.

According to the help files on building a query expression:

  Because you are selecting columns as a whole, you cannot restrict the SELECT to return only some of the columns in the corresponding table because the SELECT * syntax is hard-coded. For this reason, keywords, such as DISTINCT, ORDER BY, and GROUP BY, cannot be used in an SQL query in ArcGIS except when using subqueries.

In this SO question it shows how to achieve this using a subquery, but it seems to depend on the WHERE clause being part of the subquery.


texts: No output in PyQGIS standalone script running Distance Matrix?

texts:
I am trying to create a standalone PyQGIS script which runs the Distance Matrix in QGIS 3, see below:

######### INITIALISE QGIS STANDALONE ################

import sys
from qgis.core import (
     QgsApplication,
     QgsProcessingFeedback,
     QgsVectorLayer,
     QgsField,
     QgsFields,
     QgsProperty,
     QgsProcessingFeatureSourceDefinition,
     QgsProcessingOutputLayerDefinition
)

#start QGIS instance without GUI
QgsApplication.setPrefixPath('/usr/share/qgis', True)
myqgs = QgsApplication([], False)
myqgs.initQgis()

######### INITIALISE THE PROCESSING FRAMEWORK ################

# Append the path where processing plugin can be found (assumes Debian)
sys.path.append('/usr/share/qgis/python/plugins')

#import modules needed
import processing
from processing.core.Processing import Processing

#start the processing module
processing.core.Processing.Processing.initialize()

######### Set Distance Matrix plugin parameters ############################

# I used this command in the QGIS3 python console
# >>> processing.algorithmHelp("qgis:distancematrix")
# to get the object types it accepts for each parameter

inputlayer = QgsVectorLayer('/path/to/myinputfile.shp', 'layer 1', 'ogr')
targetlayer = QgsVectorLayer('/path/to/mytargetfile.shp', 'layer 2', 'ogr')

inputfield = QgsProperty.fromField("InputPointID")
targetfield = QgsProperty.fromField("TargetPointID")

matrixtype = QgsProperty.fromValue(0)
nearestpoints = QgsProperty.fromValue(0)

outputlayer = QgsProcessingOutputLayerDefinition('/path/to/myoutputfile.csv')

params = {
    'INPUT' : inputlayer,
    'INPUT_FIELD' : inputfield,
    'TARGET' : targetlayer,
    'TARGET_FIELD' : targetfield,
    'MATRIX_TYPE' : matrixtype,
    'NEAREST_POINTS' : nearestpoints,
    'OUTPUT' : outputlayer
}

feedback = QgsProcessingFeedback()

res = processing.run('qgis:distancematrix', params, feedback=feedback)
res['OUTPUT'] # Access your output layer

The script runs without errors but I don't get an output file. I have looked at the PYQGIS 3 API reference and searched more generally but there are so few worked examples for standalone PYQGIS3.

I have also run the command below in QGIS3 and used the field types it recommends but it still doesn't work:

processing.algorithmHelp("qgis:distancematrix")

Any ideas what could be wrong with my script?

From Failed to create memory layers in QGIS application on Linux, I added the following code to check if the input layers were loading properly.

#Input Layers Validity Check
print(inputlayer.isvalid())
print(targetlayer.isvalid())

I got a False result on both layers.

Then using PyQGIS QgsVectorLayer() Loading Invalid Layer in Standalone Python Script? I found that the input layers were not loading because I had the wrong prefix. I replaced the below.

QgsApplication.setPrefixPath('/usr/share/qgis', True)

With this

QgsApplication.setPrefixPath('/usr', True)

It is now creating a blank csv file as the output (which is progress) but its still throwing an error when it runs, see below.  Its a step forward from where it was but still not fixed yet.

Traceback (most recent call last):
  File "/path/to/distmatrixqgis3standalonenogui.py", line 99, in <module>
    res = processing.run('qgis:distancematrix', params, feedback=feedback)
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 96, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 183, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.


texts:
I'm trying to switch to IIS on my MS4W MapServer install.  It works fine in Apache but on IIS I get:

  msProcessProjection(): Projection library error. proj error "no arguments in initialization list" for "init=epsg:3857" msProcessProjection(): Projection library error. proj error "no arguments in initialization list" for "init=epsg:3857"


texts: Exporting list of tables contained in *.mdb with ArcPy gives Indentation Error?

texts:
I'd like to export to excel each of the tables in a .mdb database as separate .xml files.

The code works until the last line which gives an unexpected indentation error, and I'm also struggling with os.path.join, and how to save the excel files to an output folder.

import arcpy, os

from arcpy import env
arcpy.env.workspace = r"C:\Users\afullhar\Desktop\Lab6\HydrogetWorkCopy\HydroGet.mdb"

outFolder = r"C:\Users\afullhar\Desktop\Lab6\ExportXML"

tableList = arcpy.ListTables()

for table in tableList:

    print (table)
    in_data = table
    out_file = os.path.join(outFolder, in_data + ".xml")
    export_option = "DATA"
    storage_type = "NORMALIZED"
    export_metadata = "#"
    arcpy.ExportXMLWorkspaceDocument_management(in_data, out_file, export_option, storage_type, export_metadata)


texts:
I am very new to PostGIS (started using it this morning).

My goal is to determine if a lat/long point is contained inside a polygon that is described in a shapefile (The shapefile has a single polygon).

My plan for attack is to import the shapefile into PostGIS (this is the tool I have to use) and run a spatial query. These are my steps:

1) Infer the SRID from the shapefile:

cat CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.prj
PROJCS["Calgary_3TM_WGS_1984_W114",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["World_Geodetic_System_of_1984_GEM_10C",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-114],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

and http://spatialreference.org/ref/epsg/3776/

My best guess is that the SRID I should use is 3776.  I have no background in GIS so I am just making a best guess.

2) Convert the shapefile into .sql file:

shp2pgsql  -g the_geom_3776 -s "3776" -W "latin1"  CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.shp > CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.sql

3) import the shapefile in PostGIS.

psql -d testdb -f CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.sql

I want to determine if a lat/long point is inside the perimeter of the City of Calgary.  This are the coordinates: 51.0644894,-114.081923.  Google tells me they are within the city limits so it must be true.

However, the following query returns false:

SELECT ST_Contains((select the_geom_3776 from calgis_citybound_limit), ST_SetSRID(ST_MakePoint(51.0644894,-114.081923),   3776));
st_contains
-------------
 f
(1 row)

I have a gut feeling that my projection is wrong but I can't figure out what the problem might be.


texts: (c++) Advanced Geospatial Analysis Library

texts:
I've been searching for a lightweight library for geospatial analysis, with low overhead. Yesterday I found TurfJS, looks great, simple and fast, yet, its made in JS.

Does anyone know similar libraries, also lightweight, but for C++? I tried searching, but I'm fairly new to GIS so I have no clue on how to refer to such libs.


texts: Problem with OpenLayers tile stitching printing - All black image

texts:
I am trying to get an image of my OpenLayers application.  I am using this http://www.mapsportal.org/olprint/ and http://trac.osgeo.org/openlayers/wiki/TileStitchingPrinting as a guide.

When I click my print button I get a print failed message and when I go to the TMP folder (c:\wamp\www\tmp) there are IMG files and a JPEG but the JPEG is all black.


texts: How Can I Clip a Basemap with an ascii file?

texts:
I'm new here but with many doubts.I work with ArcGIS 10 and I add a basemap layer from the ArcGIS online but I just want  the continent coverage to my map,and the area that I wanted I have an ascii file,like a grid.So how can I cut of the ocean part?Is it the "Clip" function the better option?Is there a script for python for this?
If someone had a answer please I'll be very grateful.

Thanks.

Nara


texts: Convert from Mercator to transverse Mercator projection in C++?

texts:
I'm writing code to import images from a WMS server into some software. However, the software expects the images to be in transverse Mercator, but the WMS server provides it in Mercator. This ends up with images that look like this:

http://i.stack.imgur.com/faPzE.png

http://i.stack.imgur.com/8bl4n.png

Is there some way to transform the image from Mercator to transverse Mercator, assuming I know the EPSG and bounding box of the image?


texts: Where are the page setting for the composer in QGIS 3?

texts: When I open the composer I can't see any obvious means of changing the page setup (size or orientation).  QGIS 2 had a tab for this.


texts: Moving shapefiles to folder using ArcMap ModelBuilder?

texts: How do I move the created shapefiles (black circles) into the folder (red circle) within ModelBuilder?


texts: Stationing Callouts using VBScript Labeling in ArcMap?

texts:
I have a shapefile that I'm trying to produce a specific label for using VBScript. I'm using the values in the [MEASURE] field.

"PI: " & INT([MEASURE]/100) & "+" & (ROUND(([MEASURE] - INT([MEASURE]))*100,2))

MEASURE = 33598.2794

Current Output = PI: 335+27.94

Desired Output = PI: 335+98.28

How do I change the VBScript I currently have In ArcMap to produce the desired output?


texts:
I am trying to load a geojson file. I get an error (404) if the file doesn't exist but nothing happens if the file exists. I have tried different versions of D3 ... v3, v4, v5

console.log("hello"); // displays
d3.json("/sites/default/d3_files/json/conposcovidloc.geojson")
.then(function(data) {
   console.log("data"); // nothing
})
.catch(function(error) {
   console.log("error", error); // displays 404 error if file not found
});

What could be the problem?


texts:
I run a project in ecognition but always get a reminder of error when saving it. The error message is:

Error while saving project file 'L:\projects\New Project.dpr'


texts: Alternative overpass api for OpenStreetmap locations query?

texts:
Currently i use the overpass api to query for "leisure=playground" ways and nodes to display them on a map in a webapp.

The problem is, first of all its very slow (even with querys limited to a high zoomlevel) and i query dynamically based on the currently displayed map. Which leads sometimes to multiple querys in a short  time, when the user drags the map around. This querys often get blocked, because of "too many querys".

The only solution i found so far is: Download the OSM data, filter them locally and upload the filtered data in a DB or to the maptile provider. But to keep the data up to date i would have to do this on a regular basis. This seems to be very inefficient.

Is there any (free or paid) services, which gives fast results on OSM data querys? I wouldnt mind if the data is a few days old, but it should be updated regularly.

Whats the best practice to build some kind of an OSM-Interface for a special type of location?

Update

Example query:

http://overpass-api.de/api/interpreter?data=[out%3Ajson]%3B%0A%28%0A%20%20way%0A%20%20%20%20[%22leisure%22%3D%22playground%22]%0A%20%20%20%20%2840.72781164387366%2C-74.19745445251465%2C40.76032570468407%2C-74.09960746765137%29%3B%0A%29%3B%0Aout%20center%3B

which is

[out:json];
(
  way
    ["leisure"="playground"]
    (40.72781164387366,-74.19745445251465,40.76032570468407,-74.09960746765137);
);
out center;

~3 seconds for 5 ways (directly executed in the browser

It is used in a MeteorJS project. But i don't see any major differences with different frameworks.


texts: I would like to know if there is QGIS version for a 32bit Mac (Intel Duo core processor). The newer versions require a 64bit processor.


texts: Splitting 750km street into 10m intervals?

texts: I am using ArcGIS 10.3. The length of my street network is 745127.157m and would like to split this into 10m equal intervals so that I can calculate the slope of each segment. The split tool does not allow beyond 999


texts: Python: gdal.Warp how to crate mosaic and don't save?

texts:
I would like create a mosaic but I don't want to save it, just have in memory or as a variable and then close. I want to do that because I would like to cut out the roof of the buildings and sometimes buildings stand on 2 raster files.
from osgeo import gdal

a = 'C:/Users/ja/Inne/Desktop/73776_931613_M-34-46-C-d-3-3.tif'
b = 'C:/Users/ja/Inne/Desktop/73776_933044_M-34-58-A-b-1-1.tif'

files_to_mosaic = [b, a]

g = gdal.Warp('C:/Users/ja/Inne/Desktop/output.tif', files_to_mosaic, format="GTiff",
              options=["COMPRESS=LZW", "TILED=YES"])

With the code I can make mosaic and then save it but I don't want to save that. How can I solve that?


texts: Openlayers 3 display in Firefox and Chrome

texts:
I'm experiencing some issues with displaying OSM and GeoJSON.

With Firefox v30 I see my GeoJSON data, but OSM is invisible.
With Chrome v35 I see OSM, not my GeoJSON...

What is the explanation? (except creating "Chrofox"..)

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM()
        }),
        new ol.layer.Vector({
            title: 'operation',
            source: new ol.source.GeoJSON({
                //projection: 'EPSG:3857',
                url: 'emprise.geojson'}),
            style: new ol.style.Style({
                image: new ol.style.Circle({
                    radius:2,
                    fill: new ol.style.Fill({color: 'red'})
                })
            })
        })
    ],
    view: new ol.View({ <!--2D marche plus...2014 06 29 -->
        center:ol.proj.transform([1.22, 47.435], 'EPSG:4326', 'EPSG:3857'),
        zoom:9
    })
});


texts: Exporting a Whitebox GAT Raster to GeoTIFF format?

texts:
I have Whitebox 3.3 and I am trying to export from the Whitebox GAT format to Geotiff.  However, I do not see any tools to do this.  I have the ExportGeoTIFF.groovy in the scripts folder, but when I run that script, I get "javax.script.ScriptException: groovy.lang.MissingPropertyException: No such property: name for class:".

Could somebody please help me understand how to convert my .dep files back to geotiff?

The data is Landsat imagery.


texts: How do I only use the top layer or feature clicked on a multilayer OpenLayers3 map?

texts:
I have a map where some states are clickable, and some cities individually within those states are clickable.  Each layer needs to open a different pop-up or overlay.  Right now I'm using forEachFeatureAtPixel (I've tried forEachLayerAtPixel also).  These obviously return each layer, so then my colorbox is duplicating in the instances where a pin or city feature is on top of a state feature.

When a pin is sitting on top of another layer, I'd like to only use the top layer instead of iterating through each feature.  I'm having a hard time identifying the clicked layer, or top layer, in my logic.  I'm not sure which method would be correct for this use case.


texts:
I am trying to learn and work with PyQGIS 3.

It is a very challenging task cause documentation is really limited and out of date.

What I am trying to do is to:

create a project
add a layer
save the project
export it as a png

I have done the first 3 steps (see code below).

import sys
import qgis
from qgis.core import QgsVectorLayer
from qgis.core import QgsApplication
from qgis.core import QgsProject
from qgis.core import QgsMapRendererJob, QgsMapSettings
from qgis.gui import QgsMapCanvas
from PyQt5.QtGui import QImage, QColor, QPainter
from PyQt5.QtCore import QSize

# PYTHON APPLICATIONS

# prefix path: the location where qgis is installed in the system
# the easiest way to find the prefix path to run the following commaand from the python console of qgis: QgsApplication.prefixPath()

#1. CREATE THE APPLICATION
QgsApplication.setPrefixPath('/usr', True)
# second parameter into false disables the gui
qgs = QgsApplication([], False)
#load providers
qgs.initQgis()

#2. CREATE THE PROJECT
# create project
project_name = '/home/dkar/workspaces/qgis/data/test_project_2.qgz'
project = QgsProject.instance()
new_project = project.write(project_name)

#3. READ/OPEN THE PROJECT
# modify the project e.g. adding more layers and save it.
new_project = project.read(project.fileName())

#4. ADD A VECTORLAYER
# add a layer in the QgsProviderRegistry
layer=QgsVectorLayer("/home/dkar/workspaces/qgis/data/pakistan_1.shp", "pakistan", "ogr")

if not layer.isValid():
    print("Layer failed")

#5. ADD LAYER TO THE REGISTRY
project.addMapLayer(layer, True)
project.write('/home/dkar/workspaces/qgis/data/test_project_2.qgz')

#6. SAVE THE PROJECT
all_layers = QgsProject.instance().mapLayers()

img = QImage(QSize (800, 600), QImage.Format_ARGB32_Premultiplied)
# set image's backgroud color
color = QColor(255, 255, 255)
img.fill(color.rgb())

#create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)

# THIS IS WHERE THE ISSUE APPEARS
# THIS IS WHERE THE ISSUE APPEARS
render = QgsMapRendererJob() # THIS IS WHERE THE ISSUE APPEARS

# set layer set (am I using here the prohect?)
lst = [layer.id()] # cause "layer" is the name of the variable
render.setLayerSet(lst)

# to remove the provider and layer registries from memory
qgs.exitQgis()

But when I try to do the last part, following the instructions from the existing documentation. I get issues because the specific function (QgsMapRenderer) doesn't exist anymore in QGIS 3.

I am not sure how to solve the issue.

I see in the documentation that this function:

QgsMapRenderer

was replaced by

QgsMapRendererJob

But then when I try to import and use it in the code as:

from qgis.core import QgsMapRendererJob, QgsMapSettings
    render = QgsMapRendererJob()

I get an error message:

TypeError: qgis._core.QgsMapRendererJob represents a C++ abstract class and cannot be instantiated

Any idea how to proceed with this?


texts: Why does QGIS use negative Y spacing in the default raster geotransform?

texts:
When loading a raster dataset that does not have a geotransform, QGIS overwrites GDAL default geotransform to use a negative Y spacing instead of positive.

This is done in qgisgdalprovider, here:

if ( GDALGetGeoTransform( mGdalDataset, mGeoTransform ) != CE_None )
{
    // if the raster does not have a valid transform we need to use
    // a pixel size of (1,-1), but GDAL returns (1,1)
    mGeoTransform[5] = -1;
}

Why? It can be a problem for example if you process such raster data with GDAL, and produce vector data to display later in QGIS, they will not overlap (regardless of CRS).

I realize that the identity geotransform will not happen in real data, but the topic of the default geotransform is important enough that I'm curious why QGIS made that decision. See here or here for example.


texts:
I'm trying to do a one-time query of all points on Instagram with a given tag and visualize them in CARTO. It looks like CARTO allows you to connect a dataset via a URL. Does this mean I can simply insert the URL from Instagram's Tag Endpoints API and specify the tag name? Or is there another step I'm missing?

It looks like this would also require me to have an access token. I found this GitHub page on how to create an access token and it seemed quite involved. Is there a simpler way to get an access token? I'm not really a developer and probably won't be creating any advanced applications, but I will likely want to pull data from time to time.

I just tested this out in CARTO using a token generated from Pixel Union but it didn't work and generated the error message below. It looks like an issue having to do with the access token. Could this be that the Instagram account that the access token is associated with is private? Or would it be for some other reason?

ERROR: {"meta": {"code": 400, "error_type": "OAuthPermissionsException", "error_message": "This client has not been approved to access this resource."}}


texts: How to define a Topocentric CRS in OGC WKT form?

texts:
I'm trying to construct a DerivedCRS in OGC WKT form that describes a 3D local rectangular coordinate system that is sometimes referred to as an LSR or ENU or Topocentric.  It has an origin described by a (lat,long,ht).

My best guess so far is this:

ENGCRS[
    "Topocentric example A",
    BASEGEODCRS[
        "WGS 84",
        DATUM[
            "WGS 84",
            ELLIPSOID["WGS 84",6378137,298.2572236,
                      LENGTHUNIT["metre",1.0]]
        ]
    ],
    DERIVINGCONVERSION[
        "Topocentric example A",
        METHOD[
            "Geographic/topocentric conversions",
            ID["EPSG",4326]
        ],
        PARAMETER["Latitude of topocentric origin",
                  55.0,ANGLEUNIT["degree",0.0174532925199433]],
        PARAMETER["Longitude of topocentric origin",
                  5.0,ANGLEUNIT["degree",0.0174532925199433]],
        PARAMETER["Ellipsoidal height of topocentric origin",
                  0.0,LENGTHUNIT["metre",1.0]]
    ],
    CS[Cartesian,3],
    AXIS["(X)",east,ORDER[1]],
    AXIS["(Y)",north,ORDER[2]],
    AXIS["(Zaq)",up,ORDER[3]],
    LENGTHUNIT["metre",1.0]
]

However, gdalsrsinfo -V says "Validate Fails - SRS is not well formed".  I'm using gdalsrsinfo version 2.1.2 (2016-10-24).

I have a large number of points to convert, and will need to do this from time to time, so I would like to be able to construct the transform I need and feed it to a tool (like cs2cs) to batch convert a bunch of points at once.

How do I construct a CRS definition for the Topocentric coordinate system that I can use as input to a conversion tool, given the location for the origin of that coordinate system as a WGS-84 lat,lon,ht?


texts: ArcGIS 10.1+: Allow a model tool to run on any computer

texts:
I've created a tool that runs from a model built in ModelBuilder in ArcGIS 10.1. Right now, I can run the tool with no problem from my computer, but no one else can run it successfully.

I've changed the tool's Environment Settings to run on a Scratch database that others have access to (Fig.1), and the remaining output data is defined by the user within the tool itself.

Unfortunately, the tool will still not run for anyone else because it's still trying to access my own default geodatabase.

How might I override this - telling the tool to only try accessing the Scratch gdb defined in the tool's Environment Settings for the intermediate data, and NOT my default geodatabase that Arc normally tries to use?

Fig. 1:


texts: GRASS GIS: How do I run a batch (bash?) script?

texts:
I am in a GRASS database (in command shell), where I want to run a loop function. For example, I want to do an operation on many raster files that are in the database. I know this must be done through a txt file. What i am confused about is where I should put this file, and how exactly I should call it.

I am also a bit confused if I should be writing a batch or a bash file. GRASS is UNIX based which makes me think I should be looking to bash, but I am in windows command shell, where batch scripting applies.

Thanks for your help!


texts: How to keep leaflet tooltips within visible map view

texts:
I have a map with panning and zooming disabled, creating a bounding box. This map has a geojson layer with features and tooltips displaying underlying data.

The issue is when the user mouses over features on the outer edges of the visible map area, the tooltips are outside of the div view. Because there's no panning enabled, the user is SOL to see the tip.

How can I keep the tips within the visible map view bound by extents?


texts: Calculating Angles at which polylines meet point features using ArcGIS Desktop or FME?

texts:
I am attempting to update old sewer network data to a new format.
Manholes features (points) are collected and pipes(line features) are given attributes in the manhole if they intersect.

For example, if 4 pipes intersect a manhole, attributes for depth, diameter, material and clock face are attributed for each of the 4 pipes in the manhole attribution.
I have no problem populating the depth, diameter and material fields, but I am having a tough time populating the clock attribute.

Think of the manhole as a clock face:  North is 12, East is 3, South is 6, and West is 9.  I would like to calculate the angle at which each line intersects the manhole (azimuth, angle, or some other if it exists).
I could write a script populating the clock attribute based on these angles, but I am having a hard time finding the tool to make this calculation.
I have access to ArcGIS 10.1 and FME workbench.  Maybe I have to write a script to make this calculation.


texts: Prevent losing WMS layers when loading project

texts:
I am building a couple data packages in QGIS - that is, QGIS projects full of useful data for my colleagues to use. The problem is, if I open the project without any internet, all the WMS layers will fail to connect and they will be lost from the project. This is not a big problem for me because I can work around it, but this is certain to be annoying when for my colleagues when I pass it on to them to use.

Here is an example. There are 12 layers that are currently failing to connect when I open the project (this is currently due to something server-side).

If I click cancel, the window disappears and I can work on the project. If click "OK", I'll be warned about losing these layers first:

In both cases, I'm taken to the project, and these layers are no longer in the layers list. If I save the project, this change will be permanent, and I'll need to manually add the layers again. That's not a big issue for me, but I'm passing this on to people who are not as comfortable with QGIS.

Is there any way to stop this happening?

I see a similar question here, unsolved: Prevent QGIS from loading WMS services when opening project


texts: Clustered points legend prints as separate RGB bands

texts:
I'm using one the ArcGIS clustering samples in order to generate clustered points. This achieves the desired effect on the map and legend.
However, when I include the Print Widget in this sample, the resulting print-out has the clustered points layer split into separate bands for red, green and blue. My updated sample is at https://codepen.io/slead/pen/bGwPgeJ and prints like this:

(Hit the Export button on the print widget at the left to generate the print, then click on the PDF link in the Exported Files section of the print widget, in order to see the print-out.)
How can I print the legend as shown on the interactive map, rather than as RGB bands?


texts:
Does anyone know of a chart or table that would explain the benefits / shortcomings of the different vector file formats, eg. .shp, .gml, .kml, .gdb.  It seems that ESRI shape is still quite prolific, but since we are constantly "trading" data with other places, pieces of the shape data tend to be forgotten (.shx, etc.).

We are thinking of switching to .kmz files but what would we lose by doing so?

I don't believe all of the attribute data comes over with the .kml/kmz files.

Is this comparison possible, or am I trying to compare apples and oranges?


texts:
I want to create a form for editing attributes in QGIS. I have a lot of question and hope somebody can suggest how they work or where I can read about the next elements of form:

Classification
List
Relation References (I have found some information about this one. http://blog.vitu.ch/10112013-1201/qgis-relations, but it doesn't activate)


texts: Analyzing the landuse in overlapping buffer area

texts: I want to do landuse regression. The reference points are 200m apart. I used buffer tool in Arcgis to draw areas of 500m to know the landuse area of the 500m radius buffer area. There were overlapping. Now, if I want to know, the area of 500m radius around every reference point. What should I do? There are lots of reference points.


texts: How to pass GeoJSON object to ST_GeomFromGeoJSON

texts:
As shown in the below posted code, I am trying to print and pass the value of feature['geometry'] to the method postgreSQLHelper.queryPostgreSQLForGeoJSON
the output of the log-statement:
DEBUG:root:feature['geometry']: {'type': 'Polygon', 'coordinates': [[[6.69223420722764, 51.1329165124371], [6.69222819326484, 51.1329081193748], [6.69221180776184, 51.1329013203209], [6.69219149698336, 51.1329073049609], [6.69217949340067, 51.1329244260894], [6.69217956546811, 51.1329430447612], [6.69218526127669, 51.1329578171713], [6.692200130301, 51.1329606289999], [6.69220511560768, 51.1329581192089], [6.69221372090576, 51.1329518135339], [6.69222707128915, 51.132939215911], [6.69223420722764, 51.1329165124371]]]}

which is a valid GeoJSON format,
but when I pass the object features['geometry'] to the method I receive an error:
self.cur.execute(sql,bindings)

psycopg2.errors.SyntaxFehler: FEHLER: Syntaxfehler bei "type"
ZEILE 1: SELECT ST_AsText(ST_GeomFromGeoJSON('{'type': 'Polygon', 'co...
Why, despite the fact that the GeoJSON object is well formatted, does this error appear, and how do I fix it?
I pass the code:
for feature in featuresArray[0:1]:
    logger.debug("feature['geometry']: {0}".format(feature['geometry']))
    postgreSQLHelper.queryPostgreSQLForGeoJSON(feature['geometry'])


texts:
I have a dataset of buildings as polygons in Qgis 2.14 (image below).
I'd like to create a raster of distance from these features (ideally with arbitrary cell dimension), for which each cell expresses the minimum distance from the closest polygon.
I've searched through all processing tools (GDAL, GRASS, SAGA), but didn't find anything that could make this kind of operation.
Is there any solution?


texts: Questions regarding the processing of mosaic Landsat 8

texts:
Approximate Translation of Original Question:

A question of theoretical and practical about the mosaic of satellite images. I have 3 different months of Landsat images, which need to accomplish the end MaxVer classification. At what point should I make the mosaic of 3 scenes? before or after sorting? Follow the steps ...

Acquisition >> Georeferencing >> Atmospheric Correction >> Classification >> Mosaic

or

Acquisition >> Georeferencing >> Atmospheric Correction >> Mosaic >> Classification

Original Question:

Uma dúvida de ordem teórica e prática a respeito do mosaico de imagens orbitais. Possuo 3 imagens Landsat de meses distintos, a qual necessito realizar ao final uma classificação MaxVer. Em que momento devo realizar o mosaico das 3 cenas? antes ou depois da classificação? Segue as etapas...

Aquisição>>Georeferenciamento>>Correção atmosférica>>Classificação MaxVer>>Mosaico

ou

Aquisição>>Georeferenciamento>>Correção atmosférica>>Mosaico>>Classificação MaxVer


texts: Cropping background map by shapefile in QGIS

texts:
I would like to crop OpenStreetMap by shapefile of the country in QGIS v3.12.1.
I knew some Plugin for that in the older version of QGIS however, it is not available anymore. So what are the options?
I am a user of Python GDAL and Rasterio so if there is some option, it would be also great.


texts: identification of Forest cover area by using eCognition Software

texts:
How to find forest cover area by using eCognition software?

I'm trying to find vegetation by using NDVI in the ruleset but it takes vegetation too with the forest.

Can I use any other filter along with the NDVI to identify only the forest cover area?


texts: I am trying to find a general raster widget to use in a custom Web AppBuilder app, but I am not seeing anything out of the box. I basically just need a widget to plug into the app that would display just raster layer. My goal is to then select a raster layer from the dialog window and be able to click or mouse over that layer and have info displayed as a popup about that raster layer. My question are 1) Are there any out of the box raster widgets that fit the bill and/or 2) Any suggestions for a custom widget that would work? I am checking custom widgets from Web AppBuilder - The custom widgets list (4/5/2018) , but not seeing anything popout.


texts: Image collection to Single Image in Google Earth Engine

texts:
I'm stuck with an issue in google earth engine.

I'm working with MODIS Burned Area Monthly Global 500m, for a certain polygon area:

var MOD = ee.ImageCollection('MODIS/006/MCD64A1');
var MOD1 = MOD.filterDate('2005-06-21', '2005-09-31');
var MODlim1 = MOD.map(function(img) {return img.clip(lim)})

The problem is the image collection has three images (one for each month), and I need to convert this into a single image, which contains the burned area for three images in only one.
Anyone knows what function can be helpful in this situation?
I need to repeat this for every year (until 2018), so with a single image, it would be easier to analyze the data.


texts: Creating a Map with Spatial filtered layer on the fly

texts:
I created a web application with User accounts.
based on the User type I need to filter certain point layer with a selected cities region.
Ie , when user log in, his accessible cities (city ids are retrieved from a db) will get. With these city IDs ,
select corresponding city geometries regions from the map. After doing an Union operation on all these polygons (city layer) ,
 give resulting final polygon to a spatial  Query filter , so that certain point layers need to be shown filtered by this polygon.

So that user can only see the his accessible points features filtered by cities regions.

I am using ArcGIS JS API V.3.3 with ArcGIS Server 10
How can I achieve this goal?
Is a solution possible via REST API/ JS API?
Geoprocessing solution also accepted.


texts:
I have a pretty immense amount of Zonal Statistics as Table computations to do, and would love to find a way to grab mean and variance values from the little .dbf tables ArcGIS creates and put them neatly in spreadsheets.

I'm trying to get basic descriptive statistics on ~40 species distributions across about 6 climate scenarios, so it adds up pretty quickly.

I am a Python neophyte at best, so I don't know that I would be able to come up with satisfactory script to do this, though I have messed around with it some.

Can any of you think of another way to do this that would hopefully save a huge amount of time?  At this point I'm pretty much cutting and pasting values into a spreadsheet, which could well take me until my 80th birthday.

Other ideas: Access + pivot tables?  Trying to run the stats in R instead?

Anyway, thanks in advance for your help.


texts: QGIS - select or export point attributes by polygon feature

texts:
I have two layers, one polygons of boundaries (similar to say local government areas) and a point layer of hospitals.  I want to select all of the hospitals within a particular area (ie a particular polygon attribute).  This could be to a new shapefile but the next step is I want to export a list.

I know how to export attribute features to CSV but how do I select only the points within one boundary area?


texts: How to erase multiple polygons from another polygon?

texts: I am trying to erase around 200 county road shapefiles from a state boundary.  The files become too large to work with when the roads are merged into one state roads file.  How can I iterate the erase process so that all the roads are removed from my one state shapefile, instead of manually erasing each road one after another? Note, I have buffered the roads so that they are no longer polylines. I am using Arcmap 10.2.2.


texts: Can I find the pixels of tiff file for given lat long if I have this information?

texts:
Driver: GTiff/GeoTIFF
Files: correctly_merged.tif
Size is 44835, 26402
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (13.159040050633186,11.913294314323917)
Pixel Size = (0.000003210050390,-0.000003210074723)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  13.1590401,  11.9132943) ( 13d 9'32.54"E, 11d54'47.86"N)
Lower Left  (  13.1590401,  11.8285419) ( 13d 9'32.54"E, 11d49'42.75"N)
Upper Right (  13.3029627,  11.9132943) ( 13d18'10.67"E, 11d54'47.86"N)
Lower Right (  13.3029627,  11.8285419) ( 13d18'10.67"E, 11d49'42.75"N)
Center      (  13.2310014,  11.8709181) ( 13d13'51.60"E, 11d52'15.31"N)
Band 1 Block=44835x1 Type=Byte, ColorInterp=Red
Band 2 Block=44835x1 Type=Byte, ColorInterp=Green
Band 3 Block=44835x1 Type=Byte, ColorInterp=Blue

I found this answer in JAVA but I cannot understand it.
Anyone who can tell me in python?

Convert latitude longitude pair to pixels in geotiff


texts:
I have a script I want to use to update a value for a group of selected features in ArcMap in a string field.  I am using the UpdateCursor method but I cannot get the output generated that I want.  I know my code is close but something is off.  Here is my code:

import arcpy

fc = r'D:\_data\sidewalk.gdb\repairs'

cursor = arcpy.UpdateCursor(fc)

for row in cursor:
    if row.FixedBy == '':
        row.FixedBy = 'Contractor'
        cursor.updateRow(row)

print "Processing complete"

del row

I am trying to update the string "FixedBy" field where values are blank, not 'Null', with the value 'Contractor'.  What is my code missing?


texts: Using memory layer as input and/or output with runalg?

texts:
Can a memory layer be used as an input and/or output with runalg?

With runandload I get this to work:

 processing.runandload("qgis:mergevectorlayers","memory:bufferOne", "memory:bufferTwo", r"memory:merged")
 processing.runandload("qgis:dissolve", r"memory:merged", True, '', r"D:\PythonTesting\dissolved.shp")

Parameter 2 and 3 are the input layers, which are in memory already, and parameter 3 is the output layer, which is written into memory and used in the following tool (Dissolve). This works fine, but when I try to use the same logic with runalg, then my merged layer in memory does not seem to be created, as the Dissolve tool never runs. The following, for example, would not work:

 processing.runalg("qgis:mergevectorlayers","memory:bufferOne", "memory:bufferTwo", r"memory:merged")
 processing.runandload("qgis:dissolve", r"memory:merged", True, '', r"D:\PythonTesting\datenschrott\dissolved.shp")

Can I assume that runalg either does not accept layers that are in memory or that it cannot output any?


texts: Processing gdalogr:cliprasterbymasklayer function runs with no results

texts:
Processing gdalogr:cliprasterbymasklayer function is used in PyQGIS and it runs, but no results were obtained. Code is below:

import processing

mask = "C:/Users/Siddarth/Desktop/Kamalesh/kamalesh/New folder/Mask/Mask1.shp"
input = "C:/Users/Siddarth/Desktop/Kamalesh/kamalesh/New folder/input/dist2faults2.tif"
output = "C:/Users/Siddarth/Desktop/Kamalesh/kamalesh/New folder/output/result.tiff"

processing.runalg('gdalogr:cliprasterbymasklayer',
                   input,
                   mask,
                   "",
                   False,
                   False,
                   False,
                   0,
                   0,
                   1,
                   1,
                   1,
                   False,
                   0,
                   False,
                   "",
                   output)

Are used parameters correct?


texts: Feature Collection of random points need to be converted in Multipoint

texts:
I know that with this function I can create n random points in the area geometry.
var fg_points = ee.FeatureCollection.randomPoints(geometry,n);

But it return a Feature Collection. I need to have a Multipoint.
Is there a code that gives me a Multipoint instead of a Feature Collection, or is there a command to transform the Feature Collection into a Multipoint?


texts: How to group my plugins in QGIS menu and toolbar?

texts:
I work on and test more QGIS plugins which falls into the same category (lets call it MyProcessing). Though they fall into the same category I want to keep them as separate plugins. I would like to have these plugins in the Plugins menu sharing same menu item (My Processing) and especially sharing the same toolbar (separate from standard plugin toolbar). Is there a way to do that?

In one plugin I create toolbar in interface:

self.toolbar = self.iface.addToolBar(u'MyProcessing')
self.toolbar.setObjectName(u'MyProcessing')

How do I check in other plugin if toolbar MyProcessing exists and add icon into it?

Same with the menu item, is there a way to check if menu item already exists and add other plugin action into it?


texts:
I am having problem showing a marker on page, I have this code which is implemented on my page, in console it is not showing any errors , the map is shown also with this code, but I cant see marker.

Map


texts: proj4's omerc projection (cs2cs) broken?

texts:
When I project longitude -60, latitude 0 to mercator, I get the
expected answer:

> echo "-60 0" | cs2cs +proj=lonlat +to +proj=merc
-6679169.45     0.00 0.00

However, when I project it to oblique mercator with very similar
parameters, I get a very different result:

> echo "-60 0" | cs2cs +proj=lonlat +to +proj=omerc +k=1.0 +lat_0=0 +alpha=0.0001 +lonc=0.00001
-8416502.68     -4.68 0.00

Shouldn't the results be approximately equal?

I couldn't use exactly 0 for parameters, since omerc won't initialize there.

I sense I'm doing omerc wrong: what's the right way?

EDIT: I am aware that omerc might have a minor bug for alpha<0: http://permalink.gmane.org/gmane.comp.gis.proj-4.devel/5773 but this shouldn't affect what I'm doing.


texts:
I have a script which reads a geodatabase table and calculate two fields based on existing fields.
This is my code:

import arcpy
from arcpy import env

#overwrite output
arcpy.env.overwriteOutput = True

#Linear reference
def locateFeatures(fc,rts,tbl):
    arcpy.LocateFeaturesAlongRoutes_lr(fc, rts, 'RID', '0.1 Meters', tbl,
    'RID LINE FMEASN TMEASN')

#add field Diff MEAS
def addFieldDif(tbl):
    arcpy.AddField_management(tbl,'DifFMEAS', 'DOUBLE')
    arcpy.AddField_management(tbl,'DifTMEAS', 'DOUBLE')

#calculate Diff MEAS
def calFieldDifMeas(fc):
    fdList = ['FMEASN','TMEASN','FMEAS','TMEAS','DifFMEAS','DifTMEAS']
    with arcpy.da.UpdateCursor(fc,fdList) as cursor:
        for row in cursor:
            row[4] = abs(row[0] - row[2])
            row[5] = abs(row[1] - row[3])
            cursor.UpdateRow(row)

fc = r'C:\Users\Gary\Documents\ArcGIS\Default.gdb\PplProt'
rts = r'C:\Users\Gary\Documents\ArcGIS\Default.gdb\PplRoute'
tbl = r'C:\Users\Gary\Documents\ArcGIS\Default.gdb\Locate_MEAS'
locateFeatures(fc,rts,tbl)
addFieldDif(tbl)
calFieldDifMeas(fc)
print'done!!'

When I run the script I get the error:

  Traceback (most recent call last):
      File ", line 31, in
      File "", line 21, in calFieldDifMeas   RuntimeError: A column was specified that does not exist.

Maybe this is a very silly error but I've started learning Python recently and I cannot find a fix anywhere. Can anybody give me a hint of what might be causing the error?


texts:
I'm trying to use Ogr2Ogr.exe to import some shapefiles. It came installed with QGIS 1.8. When I try and run it, i get the following error:

C:\Program Files (x86)\Quantum GIS Lisboa\bin>ogr2ogr -overwrite -f MSSQLSpatial "MSSQL:serve...<snip>

ERROR 4: Unable to open EPSG support file gcs.csv.
Try setting the GDAL_DATA environment variable to point to the
directory containing EPSG csv files.
ERROR 6: No translation an empty SRS to PROJ.4 format is known.

C:\Program Files (x86)\Quantum GIS Lisboa\bin>

Can anyone help?

It looks like it's an issue with some PATHS.

I just opened up a normal command prompt, when to that directory (listed above) and ran the command.


texts:
I have a set of x,y data used to plot transit pickups. I am using the multiple attributes function in the symbology tab to show how many different pickup locations the agency is traveling to. I'm using latitude and longitude coordinates for my multiple attributes.

I then want to be able to make the size of the symbol dependent on the count. Is there a way to do this within ArcGIS 10.3.1? Right now pickup points that are the same just show up on top of each other.

When I go to the symbol size button, I'm only able to use latitude or longitude for symbol size which is not what I am wanting.

I'd like to be able to show the more frequented locations.


texts: Loading print composer template and nothing turns up

texts:
I am with Qgis 2.14.1 Essen, I had stored several templates on print composer, but today trying to load again those templates and the working area turns up totally blank, and the computer seems struggling to do anything in to there (check the link)

My screen when I am trying to load a save print composer template

My question is, can I recover my info and start working again with this print composer?

To clarify: I try to use the print composer with and without my projects, just in a blank new file, but still the same problem.


texts:
What I want to do is cut polygons into new euqal size polygons, so far the only solution I found was using the Parcel Editor. However, I cannot use any of the functions because everything is disabled and I can't figure out what it is that I'm doing wrong.

I first converted my polygons to lines and created a topology with the 6 required rules. I validated the topology and checked for errors (no errors found). I created a new Parcel Fabric for the feature dataset where I store my polygons and lines (enabled information model). Then I used the tool "Load a topology to a parcel fabric" to import the polygons into the parcel fabric. I am able to select a parcel but I can't do anything with it. I loaded the lines into the parcel fabric as well but it doesn't change anything. Neither does using the workflow function.
I also have some red exlamation marks in the table of contents but I don't know how to change that. (see image)

Does anyone know how to solve this or another way to split polygons into eual size smaller polygons?


texts:
I have imported a csv into ArcGIS, then displayed xy data. For some unknown reason to me, the point are plotting incorrectly. I am using an ESRI base map in WGS 84 web mercator aux sphere and the display is set to decimal degrees.  When I mouse over the correct position in ArcGIS the bottom right coordinates of the mouse match the imported csv properties.

If I search using the coordinates in Find locations in ArcGIS, it plots indicates correctly, but the same point imported in via a csv does not as shown in the below image:

Anyone have any suggestions?


texts: Building a tool to both geocode and do topology check on the geocoded point (ArcGIS)

texts:
Never built a tool in ArcGIS. Would someone tell me if it's feasible to create a single tool (preferably using Python) which:

Geocodes an address
Checks to see if geocoded point falls within a boundary
Notifies user if point does or does not fall within boundary

Is there any advice on resources to consult, or ways to go about a solution as this is outside my realm of experience?


texts: Layer Join and Styling Trouble because of Field Types

texts:
I am trying to join a .csv-file with als rates and zip3 values to a shapefile with zip3 values and then color the map based on the als rates.
This seems very simple and I have done similar maps before but I am having trouble with the field types. Whenever I import the CSV it makes all of the field types strings and I am not able to change them even using the refactoring toolbox. Then when I join them they are still strings and I am not able to style them using graduated symbols because they are not numerical.
Here is a link to the shapefiles and the csv that I am using: https://drive.google.com/drive/folders/1X26BXzrCtjSGuAPweH3HPiGshcfGaU6U?usp=sharing
I am using QGIS version 3.10


texts: Keep only one geometry type when clipping polygons to sqlite with ogr2ogr

texts:
I have a regular grid of polygons and want to clip them using ogr2ogr. While using sqlite files as destination I stumbled over some wierd behavior:

When clipping polygons exactly at their borders, a mix of geometry types is created so that I end up with points, lines and polygons while I only would expect polygons. When using an ESRI Shapefile, I get errors, thus I add -skipfailures. However, eventually, the Shapefile will only contain one polygon. The figure at the bottom shows the example data.

So my Question is: Is there any possibility to define or select the output geometry type from the clipping? I stumbled over a similar post. Adding to this, would there be any solution, how to do this in one ogr2ogr statement?

I prepared a reproducible example:

#! /bin/bash
echo "ID, GM" > poly.csv
echo '1, "POLYGON((4249000 2863000, 4250000 2863000, 4250000 2862000, 4249000 2862000, 4249000 2863000))"' >> poly.csv
echo '2, "POLYGON((4250000 2863000, 4251000 2863000, 4251000 2862000, 4250000 2862000, 4250000 2863000))"' >> poly.csv
echo '3, "POLYGON((4250000 2862000, 4251000 2862000, 4251000 2861000, 4250000 2861000, 4250000 2862000))"' >> poly.csv
echo '4, "POLYGON((4249000 2862000, 4250000 2862000, 4250000 2861000, 4249000 2861000, 4249000 2862000))"' >> poly.csv

# vrt in order for reproducibility
echo '<OGRVRTDataSource>
    <OGRVRTLayer name="poly">
       <SrcDataSource>poly.csv</SrcDataSource>
      <GeometryType>wkbPolygon</GeometryType>
 <LayerSRS>PROJCS["ETRS89 / ETRS-LAEA",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],AUTHORITY["EPSG","3035"],AXIS["X",EAST],AXIS["Y",NORTH]]</LayerSRS>
 <GeometryField encoding="WKT" field="GM"> </GeometryField>
     </OGRVRTLayer>
</OGRVRTDataSource>' > poly.vrt

# convert vrt to sqlite
ogr2ogr -f "SQLite" poly.sqlite poly.vrt

# clip twice
ogr2ogr -f "SQLite" -clipsrc 4250000 2860000 4252000 2862000 out.sqlite poly.sqlite
ogr2ogr -f "ESRI Shapefile" -clipsrc 4250000 2860000 4252000 2862000 out.shp poly.sqlite -skipfailures
# check the result
ogrinfo out.sqlite -al
ogrinfo out.shp -al

Output:

# the sqlite
....
OGRFeature(poly):1
  id (String) = 1
  POINT (4250000 2862000)

OGRFeature(poly):2
  id (String) = 2
  LINESTRING (4251000 2862000,4250000 2862000)

OGRFeature(poly):3
  id (String) = 3
  POLYGON ((4250000 2862000,4251000 2862000,4251000 2861000,4250000 2861000,4250000 2862000))

OGRFeature(poly):4
  id (String) = 4
  LINESTRING (4250000 2862000,4250000 2861000)
....

# the shp
....
  POLYGON ((4250000 2862000,4251000 2862000,4251000 2861000,4250000 2861000,4250000 2862000))
....

Figure: visualization of the problem.
blue dashed line: exemplary clipping extent
red: location where the point is created
green: location where the lines are created


texts: Running external program on "identify" on WMS layer?

texts:
I have a QGIS project with several WMS layers. When I select one of these, the identify-function shows a windows with more information on the feature I clicked. It gets this data from the GetFeatureInfo-request to the WMS. Now I would like to pass a part of that to an external program. All I have found so far is how to do this with actions, but these seem to work only on vector layers.

Is there a way to accomplish this with a WMS layer and its GetFeatureInfo response?


texts: Exporting many tables from PostgreSQL to GDB using ogr2ogr

texts:
So basically I imported .gdb package into PostgreSQL with the following command:
ogr2ogr  -f "PostgreSQL"   PG:"dbname='my_database'" "ntest.gdb" "layername" -nln table_layername

I have about 50+ layers in one GDB, so I run this script for all 50+ layers. After all the tables have been edited/processed, I have to dump all my 50+ tables into one .gdb again. I have tried as suggested here
ogr2ogr -f "FileGDB" errors.gdb  PG:"dbname='my_database' tables='table_layername, table_layername2, etc'"

But this gives me errors:

ERROR 1: PQconnectdb failed. missing "=" after "table_layername2," in
connection info string

If I run the same command without tables part in PG, but specifying the table in the end it gives me:

ERROR 1: Unable to find driver `FileGDB'.

If I try the existing driver "OpenFileGDB", the error is:

ERROR 1: OpenFileGDB driver does not support data source creation.

I guess it is because OpenFileGDB is only read-only, whereas FileGDB is read and write access (source).
All the information I am finding is old (2012-2015)6 so my question is:

Where to find and how to install the `FileGDB' driver? Does it still exist?
If I could install it, how can dump all my tables into one .gdb file? Will my command work if the driver is installed?

Edit: My Operating System: Debian GNU/Linux 9 (stretch);
Kernel: Linux 4.9.0-8-amd64;
Architecture: x86-64;


texts:
How can I fit the map extent into the QgsComposerMap size?

printable_layers = QgsMapLayerRegistry.instance().mapLayers().values()
printable_layers_ids = []
bigger_printable_layer = None

for layer in printable_layers:
    printable_layers_ids.append(layer.id())

    if type(bigger_printable_layer) != type(None):
        if layer.extent().width() + layer.extent().height() > bigger_printable_layer.extent().width() + bigger_printable_layer.extent().height():
            bigger_printable_layer = layer
    else:
        bigger_printable_layer = layer

printable_layers_ids.reverse()

template_document = QtXml.QDomDocument()
template_document.setContent(QtCore.QFile(Printer.printing_template_file_path), False)

map_renderer = QgsMapRenderer()
map_renderer.setLayerSet(printable_layers_ids)

rect = QgsRectangle(map_renderer.fullExtent())
rect.scale(1.1)
map_renderer.setExtent(rect)

map_renderer.setOutputSize(QSize(800, 600), 300)

composition = QgsComposition(map_renderer)
composition.loadFromTemplate(template_document)

composer_map = QgsComposerMap(composition, 0, 0, composition.paperWidth(), composition.paperHeight())
composer_map.setNewExtent(rect)
composer_map.updateItem()
composition.addComposerMap(composer_map)

# add created map to composition layers list
composition.refreshZList()

# not working
# composition.moveItemToBottom(composer_map)
# must do that instead
composition.setAllUnselected()
composition.setSelectedItem(composer_map)
composition.moveSelectedItemsToBottom()

# must refresh
composition.refreshZList()
composition.refreshItems()

composition.exportAsPDF(printed_file_path + ".pdf")


texts: How to Join Intersected Points by Location

texts:
I extracted nodes from polygon parcel layer by Geometry Tools > Extract Nodes. However there are many crossed unnecessary points with same location after extracted nodes. There is a topic like this:  Join attributes by location

Is there any way to simplify unnecessary points?


texts: Deferred functions, Leaflet with GeoJSON

texts:
I'm trying to chain functions, and I make a mistake somewhere, but I can't find where...
First, I create a new GeoJSON layer. It works.
var dataLayer = L.geoJSON(null).addTo(myMap);
Next I add some JSON data, it works too.
function getData() {
        var defer = $.Deferred();
        returnData = $.ajax({
                type: "POST",
                dataType: "JSON",
                url: myURL,                          // a normal URL
                error: function (x, e, s) {},        // handle error
                success: ,                           // empty
                complete: function (jqXHR, status) {
                                                     // some other stuff
                        dataLayer.addData(returnData.resopnseJSON);
                                                     // some other stuff
                        defer.resolve();
                        }
        });
        return defer.promise();
}

The AJAX query is quite long, about 50-100 ms. I'd like to chain the functions, because the next function should use the data added by getData(). In an other part of the code, when there is a lot of time (eg. some interaction, click from the user), nextFunc() can use the data added by getData(), so I'm quite sure, that the code of the functions are OK.
But if I'd like to call an other function immediately after getData() completed
$.when(getData()).done(nextFunc());
the nextFunc() function doesn't wait, it starts immediately after getData() started.
How can I force nextFunc() to wait until getData() is done?


texts: Not able to register spatial view with sde

texts:
I am using 10.2 ArcGIS and SOL Server 2008 R2.
I have created a cross database spatial view in sql server using following syntax:

USE [pipelines]
GO

SET ANSI_NULLS ON
GO

SET QUOTED_IDENTIFIER ON
GO

CREATE VIEW [gis].[parcels_test]
AS
SELECT     a.Shape, a.asmt as ASMT, meg_assor_rep.gis.megaview_ext.Owner, meg_assor_rep.gis.megaview_ext.StreetNum,
        meg_assor_rep.gis.megaview_ext.StreetDir, meg_assor_rep.gis.megaview_ext.Street,  meg_assor_rep.gis.megaview_ext.StreetType,
        meg_assor_rep.gis.megaview_ext.SpaceApt, meg_assor_rep.gis.megaview_ext.Community, meg_assor_rep.gis.megaview_ext.Zip,
        meg_assor_rep.gis.megaview_ext.Mail1, meg_assor_rep.gis.megaview_ext.Mail2, meg_assor_rep.gis.megaview_ext.Mail3,
        meg_assor_rep.gis.megaview_ext.Mail4, meg_assor_rep.gis.megaview_ext.Acres, meg_assor_rep.gis.megaview_ext.LandUse1,
        meg_assor_rep.gis.megaview_ext.TRA, a.floor as Floor, a.OBJECTID, meg_assor_rep.gis.megaview_ext.Notes,
        meg_assor_rep.gis.megaview_ext.StreetAddr,a.Shape.STArea() as Prcl_Size_SqFt,meg_assor_rep.gis.landuse.Descr
FROM       pipelines.gis.prcl a
INNER JOIN meg_assor_rep.gis.megaview_ext ON a.asmt = meg_assor_rep.gis.megaview_ext.Asmt
LEFT OUTER JOIN meg_assor_rep.gis.landuse ON meg_assor_rep.gis.megaview_ext.LandUse1 = meg_assor_rep.gis.LandUse.LandUseCode
GO

view created successfully however, when i am trying to register it with sde using following command it is giving me invalid entity(29) error.

sdelayer -o register -l pipeline_test,SHAPE -e a -k Geometry -t Geometry -i sde:SQLSERVER:GISDATA -D pipelines -p *** -u gis

This view includes a feature class of a parcel fabric dataset.
The interesting fact is that when i am trying to register view created on different feature classes (other than parcel fabric) using above command its getting registered successfully.

So, Its most probably the issue with parcel fabric.

Is there any specific prerequisite to follow before creating view on a parcel fabric feature class?


texts: Receiving planet basemaps in earth-engine

texts:
I am working on satellite imagery from the NICFI Program which is available on planet.com I am trying to download those basemaps in earth-engine.
I followed this setup guide and have connected my gee account with planet.com, but I don't understand how to use their order API (documentation). I tried a few combinations of item Ids but can't figure out what the item id for a basemap or a quad will be.
Can someone help me out
endpoint -
https://api.planet.com/compute/ops/orders/v2

body -
{
"name":"test order",
"products":[
    {"item_ids":["364-1203"],
    "item_type":"PSScene4Band",
    "product_bundle":"analytic"}
],
"delivery":{
    "google_earth_engine":{
        "project":"planet-webapp",
        "collection":"test_quad"
    }
 }
}

response -
{
    "field": {
        "Details": [
            {
                "message": "Item ID 364-1203 / Item Type PSScene4Band doesn't exist"
            }
        ]
    },
    "general": [
        {
            "message": "Unable to accept order"
        }
    ]
}


texts: In qgis2threejs html output select deselect the layers not working

texts:
I used qgis2threejs to output 3d model with a dem and another viewshed output.
I want to tick untick the layers for displaying and not displaying the layer.But the select unselect option in right window of qgis2threejs in not displaying and hiding the layers.

I used additional DEM option to integrate the additional layer in the output with display type Map canvas image.

[![on selecting viewShed the viewShed layer in not appearing][1]][1]

I want to tick/untick that option to show/hide the layers in the 3d model.
I googled but could find the results.

I have created a model including that viewShed but the tick/untick option to hide/show the layer is the feature i want in the output.


texts: Two different kind of stops for network analysis

texts:
I have a network dataset all set up and I'm able to add various stops.

However, I can only add one kind of 'stop'. I'd like to have two separate stops, in this case, ambulance stations and care homes.

This would then allow me to find the total drive time from my ambulance station to the 30+ care homes on my network dataset.

Is there a way to have two different kind of stops?


texts:
I have very little experience in dealing with ogr2ogr, and currently I am trying to find out if there is a way I can automate batch conversion of .shp files to KML using windows services schedule.

So far I've got the below script In a .bat file based on what I have found on this site:

for /R %%f in (*.shp) do ogr2ogr -f "KML" %%~dpnf.kml %%f

I am able to run the batch file and convert successfully in FWTools shell, but windows cmd.exe tells me that 'ogr2ogr' is not recognised as a command. Is there any way I can set it up so that cmd.exe recognises it automatically and just runs the .bat file every time?


texts: In R, set NA cells in one raster where another raster has values

texts:
I have two rasters in R.  I want to set values to NA in the first raster where the second raster has values.  I think this should be simple, using the raster package, with two RasterLayer objects raster1 and raster2, both the same extent and snapped to each other. They are 29775x29930.

I'm doing:

newraster <- raster1[is.na(raster2)]

But this seems to take an unnecessary amount of memory, and keeps crashing R.  My computer has 8GB of memory.  Is there a less memory-intensive way to do this?


texts: Query pixel value at specific location using arcgis image service?

texts:
I am interested in USGS elevation and slope data and I found arcgis service that serves this data. I want to query elevation or slope data at specific location using arcgis image server, so that data is returned in json object.
I tried this request: https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/identify?geometry={%22x%22:-13681264,%22y%22:4858686}&geometryType=esriGeometryPoint&mosaicRule=&pixelSize=0.5,0.5&f=pjson&returnGeometry=false
But it returned multiple values, I assume those are values of overlaping rasters...
Does anyone knows how can I query just one layer (elevation) inside this service?


texts: How can I change the language of the QGIS Interface?

texts:
I want the QGIS interface in English, in order to share information with the GIS community. Using another language, when chatting with the GIS community, makes things more difficult. When I installed QGIS, it automatically set the language to French. It's very frustrating not having the opportunity to select my preferred language.

I am using Windows XP in French, but I want QGIS in English.


texts: Change order of added FeatureGraphicsLayer in GraphicsDataSet?

texts:
I am adding two FeatureGraphicsLayer to a GraphicDataSet like so:

 graphicsFunctionality.GraphicsDataSet.Tables.Add(featGraphicsLayer);

It looks like it puts the last one in on the very top of all the other FeatureGraphicsLayers.

I need to add one to the bottom of the GraphicsDataSet.

Is this possible?

Is there a ways to change the order of FeatureGraphicsLayers once they are in the GraphicsDataSet?


texts:
Styling a marker is possible with  pointToLayer I don't know how to get a marker change style for mouse over and for click etc....

    var map = L.map('map', {
            zoom: 14,
            center: new L.latLng(35.948782,8.139912),
            layers: L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png')
        });

    var geojsonMarkerOptions = {
    radius: 4,
    fillColor: "black",
    weight: 1,
    opacity: 1,
    fillOpacity: 0.8
};
    var piky1=new L.geoJson(pick1,{
    pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng, geojsonMarkerOptions);
    },
        onEachFeature: function (feature, layer) {
        layer.on({
        'mousemove': function (e) {
        //I want to make the point a red syle ?
                color: 'red';
                piky1.bindPopup(feature.properties.Name);
       }});}}).addTo(map);

JsFiddle


texts: Pretty simple question though I am struggling with how to operationalize this, I want to rename a shapefile based on the number of features it contains. My initial thoughts are that I would need to create a list of the shapefiles contained in a folder, for those shapefiles use arcpy.GetCount and where count >= X rename the shapefile. I used an if statement followed by arcpy.rename and it managed to rename all the shapefiles. Thanks for your help with this one.


texts: Mapping Multiple Variables - MapInfo Professional

texts:
I am a fairly new user of MapInfo and was trying to determine if there is a simple way to do the following:

I have these first two parts down:

Plot a list of retail stores within a given geographical area.
Using a thematic map to "value" a particular zip code based on sales coming in from a particular zip code.

This is the part I am having trouble with:

Add an additional layer that will "shade" or "fill" a particular zip code that I have designated for "purchase." Basically I want to be able to shade a zip code and the shading will appear over the thematic shading.. Is this possible?

I have had trouble trying to find anything in the guide/online.

Any help would be greatly appreciated.

Thanks!


texts:
I install FWTools in Windows7 and run FWTools Shell,after type gdal2tiles,but give error 'no module named osgeo' .

also i installed OSGEO4W ,but another error to happen


texts: Convert all files in a Linux directory from NITF to PNG using gdal_translate

texts:
I have a directory called /home/john/WAMI/WPAFB_TRAIN/ with .txt files that have names like:
20091021202517-01000100-VIS.ntf.r0
20091021202517-01000101-VIS.ntf.r0

I know that the code needs to say something like:
for something in something:
   use gdal_translate to convert the files in the directory to .png format

For the life of me, I can't get any of the answers to previous stackexchange questions to work.


texts: Given two lat/longs, how can I tell if they are within 1 mile of each other?

texts:
I'm trying to implement a very efficient check to see whether two points are within a mile of each other, or not.

My current approach is to compute the Haversine distance, and then check to see if it's less than a mile.

Efficiency matters in this case because I have to compute this yes/no flag for large record sets.

I only care whether they are within a mile - nothing else about the distance matters to me.

So, what is the most efficient way to tell whether two lat/long points are within a mile of each other?

In response to the comments, I'm doing this in SQL Server.  My code is below.

CREATE FUNCTION dbo.USR_UFN_HAVERSINE_DISTANCE
(
  @LAT1 FLOAT(18)
 ,@LONG1 FLOAT(18)
 ,@LAT2 FLOAT(18)
 ,@LONG2 FLOAT(18)
 ,@UnitOfMeasure NVARCHAR(10) = 'KILOMETERS'
)
RETURNS FLOAT(18)
AS
BEGIN
  DECLARE
    @R FLOAT(8)
   ,@DLAT FLOAT(18)
   ,@DLON FLOAT(18)
   ,@A FLOAT(18)
   ,@C FLOAT(18)
   ,@D FLOAT(18)
   ;
  SET @R =
    CASE @UnitOfMeasure
      WHEN 'MILES'      THEN 3956.55
      WHEN 'KILOMETERS' THEN 6367.45
      WHEN 'FEET'       THEN 20890584
      WHEN 'METERS'     THEN 6367450
      ELSE 6367.45  --km
    END
  SET @DLAT = RADIANS(@LAT2 - @LAT1);
  SET @DLON = RADIANS(@LONG2 - @LONG1);
  SET @A = SIN(@DLAT / 2)
         * SIN(@DLAT / 2)
         + COS(RADIANS(@LAT1))
         * COS(RADIANS(@LAT2))
         * SIN(@DLON / 2)
         * SIN(@DLON / 2);
  SET @C = 2 * ASIN(MIN(SQRT(@A)));
  SET @D = @R * @C;
  RETURN @D;
END;


texts: Why is MOD17A3 GPP is showing huge values in R?

texts:
I have downloaded Gpp_1km (MOD17A3.055)data from MOD17A3 thorugh https://lpdaacsvc.cr.usgs.gov/appeears/ from the years 2000-2010. However when I opened all the tif files in R, each tif file showed the same maximum and minimum value (-2147483648, 2147483647). I am very confused about this.

How come there is no change in GPP from the period of 2000-2010?

Secondly the GPP is supposed to be in the units of kg C m−2 yr −1.

How come these values of tif are so large?

My study area shapefile is available at http://rds.icimod.org/Home/DataDetail?metadataId=3924

The following is GPP of 2010:

class      : RasterLayer
dimensions : 2804, 5303, 14869612  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.325  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : C:/Users/Halima Usman/Downloads/GPPMODIS/GPP2010.tif
names      : GPP2010
values     : -2147483648, 2147483647  (min, max)

and all the rest of the tif files are identical to this one.

I would like to know how come the GPP never changes from 2000-2010 and why are the values coming so large?


texts: Clipping OSM Standard Basemap using QGIS

texts:
I am trying to clip the state of MN out of the OSM Standard Basemap, omitting everything outside of the state boundaries.  Everything I've read about clipping talks about clipping between two shapefiles, but the OSM Standard Basemap isn't a shapefile.
How do I clip this to get rid of the basemap outside of the state boundaries while retaining it within the state?


texts: Can the histogram size in the legend of a tmap be adjusted?

texts:
Really enjoying tmap, and very happy to see the inclusion of histograms. Can that hist object be modified in any way? See below how the natural histogram object is somewhat squished.

this_name = "Cumberland"
this_county = nc_cr[nc_cr$NAME == this_name,]
tm_shape(this_county) +
  tm_fill("grey")+
  tm_shape(this_county) +
  tm_borders("black") +
  tm_shape(all_exposures_cr[this_county,],
           size=2)+ #spectral or set1 RdYlGn
  tm_dots(col="type", palette="Spectral", legend.hist=T, title="Environmental Exposure Type", scale=1.2,size=.1)+
  #tm_markers(col="type", legend.hist=T, title="Environmental Exposure Type", scale=1.2,size=.1)+
  tm_layout(frame = F, inner.margins = 0, outer.margins = 0, asp = 0, legend.text.size = 1)+
  tm_legend(legend.outside=T, legend.outside.position="right")


texts: Export map with single point feature over many features in feature class

texts: I have a feature class with > 7,000 point features and ~ 2,500 unique ID's (integers). I need to create simple map print outs for each of the unique ID's all with the same zoom extents. I've been looking into Data Driven pages, but it seems like that wont do because I can't figure out how to have each page contain only points associated with a particular unique ID. Is this something I can do in DDP or should I write something with ArcPy?


texts:
Is it possible to add vector tile sources (using MapBox's protocol buffers format, or GeoJSON, or otherwise) to QGIS? If so, how?

(My question is a bit similar to Can the new Mapbox Vector Tiles be used as a QGIS basemap? but isn't restricted to MapBox sources, and isn't just about basemaps or OpenLayers.)


texts: How can a QGIS user plugin reflect the events of UI?

texts: I am developing a plugin for QGIS 2.8 which uses the data of the active layer. It has a dockable window. When the user changes the active layer, the contents of the plugin window's widgets must reflect to the event, without any additional user interaction. Please provide some ideas for a solution?


texts: Line is not displaying where lat and long define

texts:
I have created 3 lines by importing them into QGIS 3 with the Add Layer/ Delimiter Text option and it does display on the map (I have added an OSM layer to understand where the line is located) but the line is supposed to be close to Las Vegas, NV (long: -114.7092, lat: 36.2565) but it is off the coast of Africa - see image below. The extents are:
Extent
-116.0046389999999974,36.7258830000000032 : -115.9775159999999943,36.7416350000000023
I am using the WGS 84 Pseudo-Mercator project for both the lines and the OSM map.
Here is the text file I import:
oid;Line
TOPIC;LINESTRING(-115.977516 36.725883, -115.977516 36.741635)
YES;LINESTRING(-115.977516 36.741635, -116.004639 36.741635)
NO;LINESTRING(-115.977516 36.741635, -115.977516 36.741635)


texts: Using ModelBuilder to iterate through fields in shapefile and create raster from each

texts:
I have the exact same problem as Iterate Through Fields in Model Builder [sic] at GeoNet but all the solutions proposed there involved ArcPy and I'm not familiar with that program. I would like to use ModelBuilder to do the exact same thing: iterate through fields in a polygon shapefile and, for each one, create a raster based on the value in the field.
Any suggestions?
I'm using ArcMap and here's what I have so far.

The error message reads:

ERROR 000728: Field 1 does not exist within table
Failed to execute (Polygon to Raster).

I am trying to iterate through fields, not values of a single field. The shapefile does have some ID fields that are string format though and which I don't want to convert. However, there are 26 Double fields that I do want to convert.


texts: Arcgis SDK C#: Addin running code on arcmap startup

texts:
I know that the addins are 'lazy loading', but I want to initialize some variables on startup of arcmap. This would be loading user preference data for the workflow we have in place, updating the FC with initial data etc.

Any direction on how to do this would be great.


texts: How to check if WMS layer loaded in ArcGIS with arcpy

texts:
I'm extracting image of .mxd map containing 3 layers, one of which is WMS layer, using standalone script.  For some reason (I guess wms is unavailable at the moment) in about 10% of extractions WMS layer is not presented on image.
Is there any way to check if layer is loaded (drawn)? Something like prolonging the response time or refreshing layers right before extraction (but with assurance that layer is loaded).
I use ArcGIS Desktop 10.0

mxd = arcpy.mapping.MapDocument(mappath)
df = arcpy.mapping.ListDataFrames(mxd, '*')[0]
df.extent = newExtent
mxd.activeView = 'PAGE_LAYOUT'
arcpy.mapping.ExportToPNG(mxd, pngname, resolution=100)


texts:
Disclaimer: I don't know anything about maps, please excuse me.

In our mobile app, we use Cloudmade as our map tiles provider.
For certain places, we download multiple tiles for each zoom level (level 13 through 19) to be used offline.

The thing is that sometimes you need to download tiles for thousands of places at once, making it a very lengthy download.

We are trying to optimize this process and there are a few things that I've come up with, which I don't know if they even make sense:

Cloudmade serves the tiles as png files. We could try to cache those tiles in our server, and reduce their size.
Skip some zoom levels, let's say we only download one every two (13,15,17,19).
Download a higher quality version on a higher zoom level, let's say 10 and split it up in the client.

Do any of these sound reasonable? What else could we do?

Thanks a lot!


texts: Installing the latest release of QGIS on Windows 10

texts:
I have a bit of experience using QGIS and this software has given me a lot of trouble. I previously downloaded the stable 3.10 release on my Mac and there was a bug with GRASS tools, so I installed the latest release 3.16 on my Mac using the Kyngchaos package, and the GRASS and SAGA tools both did not work. I do not know how to do a clean installation on my mac, I feel like I messed up the installation process and I don't know how to fix it.
I just bought a Windows computer (because I feel it is more compatible with a lot of software) but I am still so hesitant to install QGIS because I do not want to install it improperly.
The installation instructions on QGIS for windows are vague, and I don't know what the difference is between the OSGeo4W network installer and the standalone installers.
Which one is better for me?
I want to be able to access all GRASS tools and SAGA tools, and I want to install plugins like the semi automatic classification plugin.
Do I need to install any software manually (Pthyon, GDAL)?
Perhaps somebody could give me some tips on how to get the most out of QGIS without needing to install a bunch of different software manually. I want to get this installation right the first time around.


texts: Someone told me about a R package with the definition of objects some regions and nations boundaries. I'm trying to find the border of Veneto, an italian region. What's the name of the package? Do you know an alternative way to find that border (I'm not sure that Veneto is included in the package)??


texts: Is it necessary to use cachewrite control to use cacheread in openlayers?

texts:
the question is very simple. i want to use the OpenLayers.Control.CacheRead control to load and store tiles when the connection is fast e.g. with wifi network. then use the same device in other places where there is slow or no connection e.g. gprs network.

the examples i have found are not clear on implementation?

For my mobile application, i would like to locally cache the tiles. This i think can be done by CacheRead,

my question is do i need to use CacheWrite to use CacheRead Control or it will be handled by the browser's cache?

and in case i use cachewrite, how can i clear the cache so that new changes in server are visible?

i am trying to cache OpenStreetMap tiles.


texts: Terraserver replacment - source for US Topo maps?

texts:
Now that Terraserver has been shut down (or at least that's how it appears), is there an alternative WMS server for US topo maps? Probably too much to hope to find one that provides vector layers rather than rasters.

About a year ago this similar question appeared (actually bulk download of topo maps), but back then Terraserver was still in operation and appeared in several of the answers.


texts: Allowing view parameters separated by comma in GeoServer?

texts:
I want to allow view parameters to my view separated by comma.

        My parameters are like **Z060401010014,Z060401010095,Z060401010096**

            Sql Query `**select * from plotboundarygeo where plcd in('%plcd%')**`

            And i am Using Reg Exp like
            **^[\w\d\s,]+$**

           callback parseResponse
           maxFeatures 5000
           outputFormat text/javascript
           request  GetFeature
           service  WFS
           typeName plotgeo:vwarea
           version  1.0.0
           viewparams plcd:Z070400030091,Z070400030041

When i use this condition and run my application in my i am getting error like

  1 feature types requested, but found 25 view params specified.

Can any one help how to resolve this issue?


texts: How do I map points on a 2D map to a 3D illustration of the same area?

texts:
I'm looking for advice on how to build an algorithm that maps any point on a map to a point on another map, interpolating where needed.

Let me better clarify with an example of what I concretely need...

I have a set of elements and paths on a real top down map (think google maps or any gps-like system) and I have another non top down representation of the same area. Specifically I'm working on a typical ski resort map of piste and slopes like this one.

This second map is more like an illustration and it's typically a three dimensional representation of a section of mountains. The viewing angle is often selected so to offer the clearest view of the whole skiing complex and there's no guarantee that size or proportions are actually correct. It's more an artistic work than a precise map.

I want portable devices to be able to map themselves on this map, using the "real" coordinates they receive from their GPS chips.

I obviously assume that the "mapping" process is based on providing myself a certain number of reference points between the two maps. Corners of the map, start and finish of every pista, aerial ways, specific buildings, mountain tops... the more I feed to the algorithm the more precise it would probably be.

Obviously the tricky part is for it to be able to map everything in between hard coded reference points, by looking for the closest reference and averaging everything in between.

I'm sure this is done often in many similar contexts (think of games that constantly have to map things showing player locations or such) but I have no idea where to look for a concrete example. Any link to tutorials, sample code and the mathematical theory behind it would be great.


texts:
I have a table with 13,000 points of geographic coordinates ( degrees, minutes , seconds)

How I can create a layer of points ?

QGIS or Arcgis . it is possible to create a layer of points with coordinated geography ?

or only with decimal coordinates ?

There are many points , it will be very time consuming to convert 1-1 .


texts:
I want to download an appropriate aerial image. Since they are huge, disk-space consuming, and time-consuming to download, I want to zero in on the right one. The directory lists images for my area (UTM zone 10 -- San Francisco and Oakland and beyond) named like this, in this order:

10SEG355985.tif
10SEG370880.tif
10SEG370895.tif
10SEG370910.tif

Can't find helpful metadata. Is the "37" degrees of north latitude? That would make sense. But what is the "0985" part? And the "SEG" part -- that is a UTM convention, correct?

The website is "Cal-Atlas Geospatial Clearinghouse" -- http://www.atlas.ca.gov/download.html#/casil/imageryBaseMapsLandCover/imagery/urban_imagery/0.3meter_urban_area/overflights/san_francisco_oakland


texts: Converting raster to TIN in ArcMap 10.1 is causing crashes?

texts:
I have several tiles (~200) of Lidar 1m resolution data from Geomatics group and need to convert it to TIN format (currently in .asc) to develop a flood model.

Is there a quicker way of doing this than doing each individually?
I have reached a point where now I go to convert a file, Arc crashes and does not even produce an error report.

Currently using tiles in the format of .asc. Converting them using the 3D analyst tool - from raster - raster to TIN. Using the max number of points as the set 1500000 and z factor 1. Have reached a point where when I am trying to convert tiles, the system crashes without being able to generate a report. Also, very slow in showing the TIN drawn data.


texts:
I'm using ArcGIS and my datasets cover most of the world, but I'm wanting to focus on just Iceland and I'm wondering how you crop the layers to a specific area?

I've read about the crop tool, but that uses a polygon as the 'cutting line' but how do I actually draw a graphic shape that can be used in such way?


texts: How do you calculate drive times around a point?

texts:
For example, if I have an Airport point and I want to create a 30-minute Drive Times buffer polygon around it, how would I do that?

I'm currently using QGIS 3.4.4.


texts:
I work in ArcMap and occasionally have to export all of my data to KML for our office employees who only have access to Google Earth Pro. I just updated the parcel information for the entire district and exported it to KML. I switched the link from the older parcel info to the new info.

Now, we can't click on a parcel in Google Earth to see its attributes unless you click on the precise center of the polygon. My version of parcels is hollow with an outline, while the previous version of parcels was shaded throughout.

What am I doing wrong?


texts: Plotting contours for geodetic data.(in Python)

texts: I need to plot contour plot for a geodetic dataset, but I have data for 15 locations of New Zealand. I have been able to generate a map a New Zealand using cartopy library of Python, but I am not sure how to get the contour plot, since the 'countourf' function of matplotlib takes in a grid of longitudes and latitudes of the whole concerned area, along with the data at each grid-point. but I only have data for few(15) specific locations.


texts:
I'm working on a map that has nine unique categories and Carto adds an "Other" category to the legend. There is no other data outside of my nine unique categories and I can't get rid of Other label on the legend.

Do you have any advice for how to fix this?