texts: Setting DATE as default column type for esriFieldTypeDate in ArcGIS Desktop?

texts:
If I use esriFieldTypeDate with ArcGIS 10.2 a DATE column is created in Oracle.

If I use it in ArcGIS 10.6 a TIMESTAMP(6) is used

Is there a way to set DATE as default column type for esriFieldTypeDate?


texts: Populating label fields in print composer

texts:
Is it possible to fill text fields in a print composer template using data saved in a csv file?

For example, every template I create has the same fields which I always have to manually enter:

Client name, site name, title

If I were to enter this information into a csv file could I then point the label to the appropriate data using an expression?


texts:
I have to edit and maintain multiple SLD's. Is there a tool to do bulk edits/ updates?

I have used tools like Atlas SLD to one SLD at a time, and could use XML editing tools to potentially edit multiple but not sure if I can edit multiple SLD's using SLD editors.

Thanks,
Ramesh


texts: MapServer - what are some options for a GUI for map file creation?

texts:
What are some possible options for a GUI to manage style and map file for MapServer.

I am interested in all solutions, QGIS plugin, ArcGIS mxd or MapInfo workspace converter or else.


texts:
Goal: Use python script to find overlap between various shapefiles in folders

Inputs: folder 1 with multiple shapefiles and folder two with multiple shapefiles

Method:

#First intersection
folder1_shp1 = "C:\..."
folder2_shp1 = "C:\..."
arcpy.Intersect_analysis(folder1_shp1 + ";" + folder2_shp1,OutputFolder + "/Intersection 1.shp","ALL","#","INPUT")

#Second intersection
folder1_shp2 = "C:\..."
folder2_shp2 = "C:\..."
arcpy.Intersect_analysis(folder1_shp2 + ";" + folder2_shp2,OutputFolder + "/Intersection 2.shp","ALL","#","INPUT")

#Third intersection
folder1_shp3 = "C:\..."
folder2_shp3 = "C:\..."
arcpy.Intersect_analysis(folder1_shp3 + ";" + folder2_shp3,OutputFolder + "/Intersection 3.shp","ALL","#","INPUT")

#Forth intersection
folder1_shp4 = "C:\..."
folder2_shp4 = "C:\..."
arcpy.Intersect_analysis(folder1_shp4 + ";" + folder2_shp4,OutputFolder + "/Intersection 4.shp","ALL","#","INPUT")

Problem: So you can see it's a pretty simple code, and it works. However, when I run it and I check the output shape files, some of the output files look like they haven't been fully processed. .

So I have to go back in and manually use the intersect tool on that particular intersection.

I use this code on many folders and there is no pattern as to which intersect won't fully process. It could be the first intersection, the second and the forth, all of them. The code doesn't throw an error and it proceeds to the next intersection where it may or may not fully process that one. I have tried rerunning the code (get different output scenario where some intersections processed more, some less) and adding in time breaks between intersections (no change). I have to do these intersections for hundreds of folders so checking each output and manually fixing them is pretty awful. Any thoughts would be helpful.


texts:
I expected to get 'false' from the following statement:

SELECT ST_IsValid(ST_GeomFromText('POLYGON((0 0,10 0,10 10,10 10,0 10,0 0))'));

Instead it returns true! Even if the vertex in 10,10 is followed by another one with the same coordinates. Maybe I misunderstood something and/or I am making some mistake. Where can I find a complete explain of that makes a geometry not valid?


texts:
(I asked this question on stackoverflow, but since there were no solutions, and it has a GIS component, I thought I'd ask here).

I have an excel file with minimum and maximum projected x and y coordinates. I used the function rect to make rectangles

The code looks something like this, (the plot boundaries are different in my real script, but you get it hopefully):

op<- par(bg= "white")
plot(c(0,1000), c(0,1000), type = "n", xlab="", ylab= "", main= "poly test")
xmin=0;ymin=0;xmax=1000;ymax=1000
poly_test<- rect(xleft= xmin, ybottom=ymin, xright= xmax, ytop= ymax, lwd=3)

So far so good. But now I want to export it to QGIS or some other GIS application, or at the least have R recognize that these are spatial objects. When I try this,

poly_test2 <- SpatialPolygons(poly_test)

it gives an error. First that error was is.list(Srl) is not true. Then I turned poly-test2 into a list, and now I get "" Cannot get a slot('area") from an object of type "null"

What should I do? I've also tried setting the xmin, ymin etc as coordinates first, but it doesn't work either.

I also want to label my polygons based on an id listed in the excel file, how to do that in R? And will that label show up after exporting, in QGIS?

EDIT2:
There was a projection in the projection file,.. I forgot to set something in QGIS , I think..
Anyway, I adapted the script by including a "for loop", since I have a spreadsheet with a bunch of minimum and maximun x and y coordinates. The for loop looks something like this:

for(i in 1: length(Excel_File$column))
{

e<- extent(Excel_File$XMIN[i],Excel_File$XMAX[i], Excel_File$YMIN[i],Excel_File$YMAX[i])
}

And then the rest of the script you suggested. When I apply plot(e), I see only one polygon. How do I get to see a plot of all of the polygons I created? Is my for loop not good? Do I need to crete e first with either e<-list(), or e<- data.frame etc?


texts: Gdal: How iterate through all pixels of a raster image?

texts:
Given a topographic raster image.

How to iterate through all pixels of this raster image ?

Solution preferably with gdal, in console. Some other libraries/script could be installed.

Related: Gdal: How to get a pixel's value from raster image?, How to conditionnally assign a new value to pixels of a raster image?.


texts: In my ArcMap project I’ve added text to the map in callout boxes to identify certain districts. Is it possible switch these on and off like the data layers? I don't want them to show in everything I export.


texts:
I'm using OpenLayers v6.1.1 and I would like to ask you to confirm that my understanding is correct. Take a look at below image:

Am I right that there might be a situation that if I would take all tile images and put them next to each other then such map might be bigger to what I actually see on the screen?

Or the tile images should always match what is on the screen (nothing more)?

Why I'm asking...

I wanted to take picture of what I currently see on the screen and my idea was to simply get all tile images (get the grid using getTileGrid() and then forEachTileCoord() for my current view extent) and just concatenate them to have the map image, but the image is not what I see on the screen. Sometimes there is much more on the image edges or the image starts from the different place on the map then I see it on the screen.

Note:

I'm not using Node.js - any kind of server side (Node.js side) image manipulation will not work for me
I'm just including ol.js and the page and play with vanilla JavaScript


texts:
I'm a newbie to Spatialite so forgive me if the answer to my question is obvious...(I tried looking online but couldn't find the answer to this question).
I am trying to use the function ST_touches in Spatialite Gui but it says it does not recognize that function. Is there anything I need to install in order to use this function? And how?
Thanks so much,
Ana


texts:
I'm working with the arcgis javascript API on a webmap.

I have a feature layer on some server which I'm consuming in my webmap and now I want to get the different symbols/icons of that layer one by one in order to display it in a text next to the map (but NOT using the legend widget).
Is that possible? I'm basically looking for something like a featureLayer.getSymbols() function but that doesn't seem to exist.
Any help is greatly appreaciated!


texts: Calculating travel time along rail network in ArcGIS Network Analyst?

texts:
I have a custom railway line network with stations in ArcGIS 10.3. The numbers between each station represent travel time between the two stations.

I would like to be able to select a station and have returned a selection of all stations within a set journey time away.

I am using ArcGIS 10.3 for Desktop and I am very new to Network Analyst.


texts: Getting AddIn to appear in customize mode categories?

texts:
I was hoping that someone here might have tried this, it seems like it could be a simple error somewhere. I was running the sample AddIn that was provided by ESRI and it worked fine until I tried to add it to ArcMap.

Once I added the Simple_Add_In esriaddin file to ArcMap it would compile and I would be able to add the toolbar component of the sample but all the buttons on the toolbar would appear as "missing" instead of buttons as well as the commands would not appear in the customize menu, and I cannot add them as buttons.

I have installed ArcObjects into Eclipse IDE.


texts: QGIS startup error (out of range pointer)

texts:
I downloaded QGIS (for Windows) and installed the software. However, when I start QGIS, I get an error message:

"The Instruction at 0x04c20d55 referenced memory at 0xffffffff. The memory could not be read."

Any ideas how to solve this problem??


texts: How can I join values of Point features to raster dataset in ArcGIS?

texts: I have to join point values to rasters. Every points are located at the geometrical centre of the raster cells. Every single cells have to include the point values that are in their centre. Unfortunately, spatial join works only with vector features. I used point raster tool and after that I used Mosaic to new raster tool because the rasters are not in the same layer. After the mosaic the new raster, I am not able to join because there are multiple rasters with the same ID. I need exactly the opposite of what extract values to points tool does. Namely, I want to extract point values to raster. Does somebody know how can I do that?


texts: Neighborhood analysis based on occurrence of value

texts: How can I do, either in ArcGIS, grass or another tool, the following: for each cell (pixel), if there is at least one cell with value '4' in a neighborhood of 5 cells around it, reclassify that cell to 0, otherwise don't change the value? Thanks.


texts: QGis : superimpose 2 rasters georeferenced (adjustment)

texts:
I am pretty new in QGis and I encounter the following problem.

I have 2 rasters, let's say, "cadaster" and "tranchot". The raster "tranchot" was georeferenced from the raster "cadaster" with QGis. After the georeferencing, QGis shows the 2 rasters perfectly superimposed and built on the same projection.

Today, I want again to superimpose this 2 rasters but "tranchot" looks like a thumbnail in the corner of "cadaster". From the metadata, the pixel makes 0.1143 meter for "cadaster" and 0.02831 meter for "tranchot".

Probably there is an handling to be done in QGIS to superimpose these two rasters, probably simple, yes, but which one ?

In the future, I'll need to compare vectors from "cadaster" with the raster "tranchot". Is there a mean to make rasters and vectors superimposable permanently ?

Many thanks in advance for your help


texts: arcpy.SetParameter() - Problem passing escaped string as geoprocessing service output

texts:
I have a geoservice that get some text stored from a table in order to return a json string to a javascript client.

However, the arcpy.SetParameter() seems to have problem to returns the escaped new lines character (\\n)

my_dict = {"key1" : "test is a string with\nsome escaped characters\nsuch newlines"}

my_json = json.dumps(my_dict)
print(my_json)
>>> '{"key1": "test is a string with\\nsome escaped characters\\nsuch newlines"}'

arcpy.SetParameter(0, my_json)

in the javascript console the response is the following which is not a valid json:

{"key1" : "test is a string with
some escaped characters
such newlines"}

Is there a solution? Someone can confirm that the problem is with arcpy.SetParameter() function?


texts: How to handle spatial join 10 chars field name limit?

texts:
When I try to perform a spatial join between a point and polygon file, I get a notice saying spatial join will not be performed because following field names exceed 10 Characters. Kindly help to rectify this.

Regards.


texts: Setting up MapServer with existing XAMPP

texts:
I am trying to set up MapServer with my existing XAMPP on Windows.
I downloaded ms4w_4.0.5.zip and extracted it into my C://-folder next to my xampp-folder. Now I would like to use XAMPP to run my MapServer WFS, but all documentation I found so far is always for setting up MapServer with the Apache server within the zipped folder...
I found this question: How to install mapserver on Windows/Apache? from 5 years ago, but it did not really help me. Does someone have something like a step-by-step guide to also help others with this in the future?


texts: how to get the boundaries of a city in KML from google maps

texts:
im trying to write a google maps application. and i wanted to only highlight one city in the world ( im using the masking method explained: https://github.com/vasile/geomask)
but now im bugged with the question of how to only show the boundaries of the city as displayed in google maps.
can any one suggest how to extract the boundaries of the city from google maps in KML file so that i can use the same. in my project.


texts:
In a python script, I can't access the geometry of my tables :

SQL_req = 'SELECT id_parcel, id_melba FROM parcels, network WHERE st_intersects(parcels.geometry, network.geometry)'
grass.write_command("db.select", input='-', stdin=SQL_req, stdout = DEVNULL, stderr = DEVNULL)

gives :

DBMI-SQLite driver error:
Error in sqlite3_prepare():
SELECT id_parcel, id_melba FROM parcels, network WHERE st_intersects(parcels.geometry, network.geometry)
no such column: parcels.geometry

Shapefiles were imported with v.in.ogr and other sql requests are correctly ending. That's when I try to use spatial functions that it don't work.

Any suggest or example?


texts:
Basically I have a permalink which contains the co-ordinates of a temporary layer when I enter this permalink into the browser the layer shoes and then goes. I used use a vector layer the permalink used to work fine but then I rasterized the layer after doing that the layer does not work. The permalink which contains the layer  looks something like this:

r=POLYGON((15523179.169885%20-4076448.6241544%2C15482208.92273%20-4058715.2335947%2C15490464.121784%20-4097545.2439582%2C15523179.169885%20-4076448.6241544)).

After rasterizing the layer I store the layer in the database before it was simple a  vector layer.

I am using geoserver, openlayers and javascript. The permalink is created using Jquery.

I am not sure who to make the layer show permanently when the user enters the permalink into the browser


texts:
Before I start, this is going to be a general question but with specific cases. I've been trying to understand GIS tools and applications to achieve the following:

Build my own basemap
Publish the basemap online
Add my own data layers to the online basemap
Perform analysis (on-the-fly)
Add interactivity to the map

So, In my first attempt, I built a spatially enabled database on PostGIS, queried this using PHP, overlayed it on google maps using their API and added interactivity using their api itself. I carried out analysis like area, denisty etc. on the fly using postGIS functions.

For my next attempt, I downloaded TileMill and rendered some really beautiful maps and added zoom level functions to it. I also exported my MBtile files. When I came across trying to publish this online, I hit quite a few names, but got lost and confused:

TileStache
TileCache
TileStream
Mapserver
Mapnik
GeoServer
Mapbox

So, I understand these are mapservers. But what exactly does a mapserver do? Why is it used? Why can't i server my tiles from a normal web Apache server? And for MBtiles, how many servers are out there which lets you publish it directly?

I can also export my basemap as a png/vsvg format. Can these formats be used to make a basemap through a mapserver?

When researching interactivity I came across:

Leaflet
ModestMaps
OpenLayers
Mapstraction

I understand these are Javascript Mapping API's. But on what base map can I use these? Can it be used on Google maps? Can I use them on my own base maps too?

I basically want to understand the differences between these tools and which would be the best software stack to build a fully customized, interactive map with my own basemap.

Sorry if my question sounds to general. I will definitely be here to clarify my questions


texts: How to remove narrow 'notches' in polygons?

texts:
I am currently working on a project where i am doing a union on several feature classes. The union normally creates several polygons which share edges/borders. however after running dissolve on the union output, i have noticed it also generates these small "alleys" or notches into the polygons:

(zoomed in)

(zoomed out)

I was wondering if there is a way to remove these narrow notches in the dissolved feature?

So far i have tried Simplify Building (ArcGIS Pro) with a simplification tolerance of 100 and 1000 meters and a minimum area of 0. This didnt resolve the issue. I would like to avoid doing this manually in an editing session, as there are many of these.


texts: GeoServer not infilling Polygons in WMS?

texts:
I am trying to make a WMS service for complicated Polygons using geoserver pointed at postGIS.

When some of the polygons are being rendered by the WMS the polygon fill crosses the geometry of the polygon, see the image for a better idea of whats going on.

Why is the polygon not being styled correctly?

My geoserver version is 2.7.6, the polygon has a srid of 3857

SLD for the Polygons in question :

<ogc:PropertyIsEqualTo>
              <ogc:PropertyName>feature</ogc:PropertyName>
              <ogc:Literal>Depths</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <se:MinScaleDenominator>300000</se:MinScaleDenominator>
          <se:MaxScaleDenominator>1000000</se:MaxScaleDenominator>
          <se:PolygonSymbolizer>
            <se:Fill>
              <se:SvgParameter name="fill">#a5bfdd</se:SvgParameter>
              <se:SvgParameter name="fill-opacity">1.0</se:SvgParameter>
            </se:Fill>
            <se:Stroke>
              <se:SvgParameter name="stroke">#000001</se:SvgParameter>
              <se:SvgParameter name="stroke-width">1</se:SvgParameter>
              <se:SvgParameter name="stroke-linejoin">bevel</se:SvgParameter>
            </se:Stroke>
          </se:PolygonSymbolizer>


texts: How to change the extent of a raster after converting its coordinates?

texts:
I changed the coordinate system from Conical coordinates to lat/long UTM-39N. It's OK, but the new layer doesn't have appropriate extent, so it can't be plot on other layers.
I exported the raster file with "dataframe extent and spatial references", but it doesn't work. Also, I used define project and raster project tools. still doesn't work.
what should I do?


texts: Labels not working on remote Mapserver installation

texts:
I have Ubuntu 14.04/Mapserver7/apache2 installed on both my dev machine and on a remote server. As far as I can tell both installations are identical. For some reason labels are not showing up on the remote version. I built the mapserv executable on each installation with identical configurations so all the required libraries should be there. Any suggestions for troubleshooting?

MapServer version 7.0.0 OUTPUT=PNG OUTPUT=JPEG SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=WFS_SERVER SUPPORTS=WCS_SERVER SUPPORTS=FASTCGI SUPPORTS=GEOS INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE

I have identical FONTS directory on each installation using truetype fonts. There are no errors showing in the mapserv error log.

It is servicing a WMS request.


texts: Landsat 8 - Surface reflectance panchromatic band ?

texts:
I ordered landsat 8 surface reflectance high level data products via Espa: http://espa.cr.usgs.gov/ordering/new.

Is it correct that there is no surface reflectance panchromatic band available after processing?


texts: How can I center Openlayers div on screen?

texts:
I have a div centered on screen (using width:75%;height:50%;margin: 0 auto;)
When I init the map div with OL , it auto aligned to the left.
How can I keep having the div in the center of the window?


texts: How to use cartodb.js to create layer - not just tilelayer - from named map

texts:
This gis.stackexchange question asks how a layer can be displayed after initializing a named map.

I don't have enough reputation to ask a follow-up comment to javisantana's response.  It reads in part:

Why didn't you use cartodb.js to get that layergroup? if you are going to use
cartodb.js it's better to use it to get that layergroup

My question is:  how do you use cartodb.js to do that?  I can do it just fine using the instructions on the Maps API page and can also successfully create a leaflet tileLayer in a fashion similar to that which javisantana's suggests.  However, I would love to be able to use cartodb.js and then be able to do a bit more than just display a static tileLayer (e.g. replace the named map's template variables at client run time).

It seems like it might be possible to do this calling cartodb.createLayer, but my attempts to provide layermetadata have not been successful.  The documentation on the cartodb.js page is a little unclear:

cartodb.createLayer(map, { layermetadata })
Layer metadata is always in the form: { type: 'LAYER_TYPE_NAME', options: {...} }


texts: Using Postgis spatial index in st_distance operation based upon explicit SRID column

texts:
I have two large tables with data that spans over multiple regions of the U.S and needs a defined SRID column in atleast one the tables to denote which coordinate system should be used for area/distance calculations:

CREATE TABLE table_a(
  id uuid,
  geom geometry(MultiPolygon,4326)
);
CREATE INDEX ON table_a USING GIST (geom);
CREATE INDEX ON table_a USING GIST (st_transform(geom,2876));
CREATE INDEX ON table_a USING GIST (st_transform(geom,2877));
CREATE INDEX ON table_a USING GIST (st_transform(geom,2926));
CREATE INDEX ON table_a USING GIST (st_transform(geom,...));

CREATE TABLE table_b(
  id uuid,
  geom geometry(MultiPolygon,4326)
);
CREATE INDEX ON table_b USING GIST (geom);
CREATE INDEX ON table_b USING GIST (st_transform(geom,2876));
CREATE INDEX ON table_b USING GIST (st_transform(geom,2877));
CREATE INDEX ON table_b USING GIST (st_transform(geom,2926));
CREATE INDEX ON table_b USING GIST (st_transform(geom,...));

--note there can be multiple applicable SRIDs for a given geometry
CREATE TABLE table_a_srid(
  id uuid,
  table_a_id uuid REFERENCES table_a,
  srid integer
);
CREATE INDEX ON table_a_srid(table_a_id)

The following query is designed to apprehend the closest record of table_b for each record of table_a within a distance of 500ft using the SRID from table_a_srid:

SELECT DISTINCT ON (a.id)

  a.id,
  b.geom

FROM

  table_a a
  JOIN table_a_srid c ON (a.id = c.table_a_id)
  JOIN table_b b ON (st_dwithin(st_transform(a.geom,c.srid),st_transform(b.geom,c.srid),500))

ORDER BY

  a.id, st_distance(st_transform(a.geom.c.srid),st_transform(b.geom,c.srid))

Query plan shows that a sequential scan is used instead of the relevant spatial transformation index for the given SRID -resulting in completely unacceptable performance.

How can I perform this operation efficiently without resorting to UNION ALL for every single SRID defined?

postgis_full_version = "POSTGIS="2.2.0dev r12887" GEOS="3.5.0dev-CAPI-1.9.0 r3995" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.0, released 2014/04/16" LIBXML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER"


texts:
A client has a very old Certificate of Ownership document for a plot of land.
We need to show the outline of the plot on Google maps to potential developers.

Sadly the description of the plot is in archaic coordinate system. How do I get this represented on Google maps?

Here is an additional piece of information I could scrape from the document:

Those are some form of coordinates for point A.4589 in the diagram.

Seems I've made some good progress

Map Location


texts: Syncing Google Earth My Places data in OS X?

texts:
Since Google hasn't implemented a sync service for Google Earth's My Places data, I'm trying to make this work manually between my home and work computers, one OS X the other Windows.  Windows side works fine - I've just redesignated the KML path in registry (HKEY_CURRENT_USER\Software\Google\Google Earth Plus\KMLPath) to a Dropbox.  The problem's on the Mac side.  I've tried replacing the myplaces.kml with a symlink, but this failed as the new object is somehow automatically overwritten by a backup (even when the backup in the same folder (/Users/username/Library/Application Support/Google Earth/) has been manually deleted).  So I'm trying by adding a 'Network Link' folder to My Places, linking to the dropbox kml.  This half works - edits from Windows side load up in Mac side, but edits/additions in Mac side placed in the networked folder are not saved upon exit.

Does anyone know the reason for this failure to sync, or of an alternative method that works on OS X?


texts: Install QGIS 3.4 with GRASS support in Ubuntu 18.04

texts:
I can install QGIS3.X but cannot run a GRASS process. I get the popup:

  This algorithm cannot be run :-(
  GRASS 7 can't be found on this system from a shell. Please install it or configure your PATH environment variable.

I have been putting repos from the following into my source.list https://qgis.org/en/site/forusers/alldownloads.html#debian-ubuntu
and after apt-get update and apt-get upgrade, using the following command to install

sudo apt-get install qgis python3-qgis qgis-plugin-grass

NB: I did have a previous (working) QGIS installation (2.18 in Ubuntu 16.04) that had I upgraded to 3.6 (after upgrading to 18.04) using the https://qgis.org/debian repo which I tried removing in order to reinstall using

sudo apt-get remove -y qgis python-qgis qgis-plugin-grass grass-gui grass grass-dev
    sudo apt-get purge qgis python-qgis qgis-plugin-grass grass-gui grass grass-dev
sudo apt-get autoremove

setting the repo in my sources.list file (to one from the link above) and updating.

I still get the popup.

After every attempt I also tried to run GRASS (or grass74) from the console and each time receive:

File "/usr/bin/grass", line 1673
sys.stdout.write("%s\n" % grass_version)
                                       ^
TabError: inconsistent use of tabs and spaces in indentation

(or similar for grass74)
I don't normally run GRASS from the command line  so I don't know if this is what I should expect. (I guess not, but can't tell, and can't find an example online of what I should find, or means to solve this).

At this point I'm stumped as to the problem or pathway to solution barring reinstalling 18.04 from scratch (not preferred).

Here are some of the things I've thought of but can't answer:

Is my purge process incomplete? What else could I check or do here?

Do none of the repos currently provide a version of qgis and grass that work together? Is there an alternative working repo?

Do people need to manually locate grass and input the location info into their system path or some parameter setting in qgis as suggested in some outdated stack-exchange queries? If so how?/where?


texts: Error while rendering Leaflet map on Shiny

texts:
I have been trying to reproduce a Leaflet map on a Shiny webpage. I have made this code so far:

library(shiny)
library(shinydashboard)
library(leaflet)
library(rgdal)
library(dplyr)
library(DT)

Map <- readOGR("C:\\Users\\XXXXX\\Documents\\NationalElectorateMap\\NationalElectorateMap\\national-esri-fe2019\\COM_ELB_region.shp")

webdata_visitors <- read.csv("C://Users//XXXXX//Desktop//webdata_visitors.csv")

webdata_visitors <- webdata_visitors[order(match(webdata_visitors$Electorate, Map$Elect_div)),]

bins <- c(0, 100, 500, 1000, 10000, 350000, 50000, 75000, 200000, 500000)
pal <- colorBin("RdYlBu", domain = webdata_visitors$user_id, bins = bins)

ui <- dashboardPage(
    skin = "blue",
    dashboardHeader(title = "Users"),
    dashboardSidebar(),
    dashboardBody(
    fluidRow(width = 12, box(leafletOutput(outputId = "mymap")))))

server <- function(input, output) {

    data_input <- reactive({input$webdata_visitors})

    labels <- reactive({
        paste("
", data_input()$Electorate,"
",
                    "
", "Number of Visitors: ", round(data_input()$user_id, digit = 0), "
",
                    sep = "")
    })

    output$mymap <- renderLeaflet({
        ElectMap <- leaflet()%>%
            addTiles('http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png') %>%
            setView(-25.2744, 133.7751, 4) %>%
            addPolygons(data=Map,
                        weight=1,
                        smoothFactor = 0.5,
                        color = "white",
                        fillOpacity = 0.8,
                        fillColor = pal(data_input()$user_id
                        ),
                        label = lapply(labels(), HTML)) %>%
            addLegend(pal = pal,
                      values = data_input()$user_id,
                      opacity = 0.7,
                      title = "Number or Visitors",
                      position = "topright")
       ElectMap })
    }

shinyApp(ui=ui,server=server)

However, when I run this code I see this error:

  Warning: Error in as.vector: cannot coerce type 'closure' to vector of
  type 'list' [No stack trace available]

I am having trouble interpreting what this error means and I am not sure how to fix my code. Could anyone help me out?


texts: Characters displaying correctly in File Geodatabase with ArcMap 10.1 but incorrectly with 10.2?

texts:
I work with data from a variety of European countries, and am having trouble getting the non-standard characters (e.g. Turkish alphabet) to display correctly in my geodatabase tables.  Sure, i expected this, but when I was running ArcGIS 10.1, they all displayed just fine!  I tried the registry edit workaround suggested on the ESRI website, to no avail.  Anybody have any better suggestions for reverting all of my '?' characters to their correct letters?

And no, i have no idea which character set they came from originally, nor which coding system they should be using.. no idea at all.


texts: What would be a good symbol for the map legend itself?

texts:
I am reformatting a mobile web map. Many of our users are equipped with mobile phones rather than large tablets, so screen real estate is at a premium. I've decided to remove the legend from default view for these devices. Users will have to click on or hover over a small control button to view the legend.

So naturally, I now need a symbol that conveys "legend" to place on the control button. Leading candidates:

A key: intuitive, but it says "key" more than "legend". A difference many users might not care about, but I do...
A question mark: it's engaging and implies it can answer questions about the map.
Will Smith
The word "legend". This is a last resort. The other control buttons have icons, so this would kind of throw a wrench in the works, but it definitely conveys the concept of a legend.

I'm not completely satisfied with any of these, which is why I've come to you. What says legend without saying legend?


texts: Creating label with number circle and field "1v" in ArcMap?

texts:
How to create a label which combines a circle number and field "1v" in ArcMap?
For example:

I need, a line the labels. My code label expression is:

Function FindLabel ( [TOTAL_VIV] )
FindLabel = "<FNT name='Arial' size='6'>"&[TOTAL_VIV]&"v"&"</FNT>"
End Function


texts: qgis2leaf exporting transparent polygons as all black

texts:
I am creating a qgis2leaf webmap from QGIS 2.8.1.  It is comprised of several polygon and point layers.  I am able to make my polygon layers semi transparent, but not 100% transparent.  I want to be able to show only the border of my polygon and have a transparent fill.  Howevere, when I set this my webmap shows all black polygons that should have a transparent fill.  I would just convert it to a line feature but I want to be able to click the polygon to view the features.  Is this a bug in qgis2leaf or is there a workaround?

QGIS Map with transparent polygons.

qgis2leaf webmap with all black polygons.


texts: How to calculate area of the output of r.water.outlet?

texts:
I wrote the following code to use r.water.outlet module in GRASS GIS:

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
gsetup.init(gisbase, gisdb, location, mapset)

#csv file reading and importation
import csv
rows = list(open('tc_sta.csv'))
totalrows = len(rows) - 1

#loop through the csv(coordinates) file in r.water.outlet module
f = open('tc_sta.csv', 'r')
element = list(csv.reader(f))
i = 0
j = 0
while True:
    if i <= totalrows:
        g.run_command('r.water.outlet', drainage = 'dra', basin = 'b' + str(i + 1), east = element[i][j], north = element[i][j + 1])
        i = i + 1
    else:
        break

#print all files in grass database to check results of the module
print g.read_command('g.list', _type='rast')

I stored the coordinates in tc_sta.csv, and input them to r.water.outlet to calculate drainage area of each point.

Although the output raster files were fine, I had to export them as ESRI shapefile and use QGIS to calculate area of polygons, and finally store the area value of each point in a csv file containing coordinates and its related drainage area.

I felt that the procedure was not efficient, and is there a way I could calculate and store the output area immediately after I run r.water.outlet?


texts: Calculating chopper rescue time with network analysis

texts:
I have a layer with locations for 3 different fire stations and one hospital with a helipad. I have created a network dataset to calculate how long it takes for the firemen to arrive to various locations around the map layer.

Now I want to calculate the rescue time for the chopper from the helipad assuming it flies at a speed of 100 knots. However I'm unable to find a way to do this. I have googled and found ways to calculate flight times but none that involves the network analysis tool.

Any ideas?

EDIT: Example for what I want to do, but with chopper times instead.


texts:
I am trying to run a spatial join on ArcPy. The files are quite large, and I am running the code in Idle. The code seems to "stop" without producing an error.
My code is as follows:
a = shapefile1.shp
b = Default.gdb\\shapefile2
y = shapefile2.shp
arcpy.SpatialJoin_analysis(a, b, y, "JOIN_ONE_TO_MANY", "KEEP_ALL")
print("finished join 1")

The code stops without printing "finished join 1". However it does not produce an error. I am not sure if it finishes the join either, it produces a large file, but I have no way of checking if it is complete.


texts: Using Windows Batch File (.bat) for executing QGIS Python program without opening QGIS Software

texts:
I have a Python Program for doing some geoprocessing in QGIS. The starting of the Python program is given below.

import sys
from qgis.core import *
import qgis.utils
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import Qt
from qgis.core import QgsApplication
from PyQt4.QtGui import QApplication
QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
QgsApplication.initQgis()
sys.path.append('C:/OSGeo4W/apps/qgis/plugins')
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *
from PyQt4.QtCore import QFileInfo
import processing
import os
import ogr
from osgeo import ogr, gdal
import math
import xlrd

..............................many other lines.................

I want to execute this python program from a Windows Batch File, that is .bat file; even without opening QGIS software.

For that, I created a .bat file as given below.

SET QGIS_PREFIX_PATH=C:\OSGeo4W\apps\qgis
call "C:/OSGeo4W/bin/o4w_env.bat"
SET PYTHONPATH=C:\OSGeo4W\apps\qgis\python;%PYTHONPATH%
SET PATH=%PATH%;C:\OSGeo4W\apps\qgis\bin

cmd /c python "C:/Users/Sreeraj/.qgis2/processing/scripts/Task3.py"

Unfortunately, when I am executing this .bat file, in C:\WINDOWS\system32\cmd.exe , I am getting an error:

python.exe has stopped working
A problem caused the program to stop working correctly.
Windows will close the program and notify you if a solution is available.

What modification do I have to make in .bat file or in my Python program, in order to solve this issue ?


texts:
I am trying to remove all the xyz connections in a QGIS3 project via PyQt5. I am hitting the following error for the code bellow:
for group in getGroupXYZ(self):  NameError: name 'getGroupXYZ' is not defined
 from qgis.PyQt.QtCore import QSettings

 s = QSettings()
 xyz = 'qgis/connections-xyz'

 def getGroupXYZ(self):
     s.beginGroup(xyz)
     groupsXYZ = s.childGroups()
     s.endGroup()
     return groupsXYZ

 def removeAllConnectionsXYZ(self):
     for group in getGroupXYZ():
     s.beginGroup(xyz)
     s.remove("")
     s.endGroup()

 self.dlg.wipeAllConnectionsXyz.clicked.connect(
     self.removeAllConnectionsXYZ)


texts: When merging tif-files - set crs to EPSG:25833 - ETRS89 / UTM zone 33N

texts:
When mering tif-files using rasterio merge, i need to set the crs to be EPSG:25833 - ETRS89 / UTM zone 33N. CRS are read using QGIS -> properties.
I use the script below, and it gives me EPSG:8687 - Slovenia 1996 / UTM zone 33N.
The code is based on https://automating-gis-processes.github.io/CSC/notebooks/L5/raster-mosaic.html
from tempfile import mkdtemp
import rasterio
from rasterio import Affine
from rasterio import windows
import math
import numpy as np
import os
import glob
#File paths
dirpath = "C:\Python\Raster\Mosaic"
out_fp = os.path.join(dirpath, "Mosaic.tif")

#Select the tif files
search_criteria = "d*.tif"
q = os.path.join(dirpath, search_criteria)
print(q)
C:\Python\Raster\Mosaic\d*.tif

dem_fps = glob.glob(q)

#List fileds:
dem_fps
['C:\\Python\\Raster\\Mosaic\\dtm1_33_103_117.tif',
'C:\\Python\\Raster\\Mosaic\\dtm1_33_103_118.tif']

sources = [rasterio.open(raster) for raster in dem_fps]

memmap_file = os.path.join(mkdtemp(), 'test.mymemmap')
first = sources[0]
first_res = first.res
dtype = first.dtypes[0]

# Determine output band count
output_count = first.count
# Extent of all inputs
# scan input files
xs = []
ys = []
for src in sources:
    left, bottom, right, top = src.bounds
    xs.extend([left, right])
    ys.extend([bottom, top])
dst_w, dst_s, dst_e, dst_n = min(xs), min(ys), max(xs), max(ys)

out_transform = Affine.translation(dst_w, dst_n)
# Resolution/pixel size
res = first_res
out_transform *= Affine.scale(res[0], -res[1])
# Compute output array shape. We guarantee it will cover the output
# bounds completely
output_width = int(math.ceil((dst_e - dst_w) / res[0]))
output_height = int(math.ceil((dst_n - dst_s) / res[1]))
# Adjust bounds to fit
dst_e, dst_s = out_transform * (output_width, output_height)
# create destination array
# destination array shape
shape = (output_height, output_width)
# dest = np.zeros((output_count, output_height, output_width), dtype=dtype)
# Using numpy.memmap to create arrays directly mapped into a file
dest_array = np.memmap(memmap_file, dtype=dtype,
                   mode='w+', shape=shape)

dest_profile = {
    "driver": 'GTiff',
    "height": dest_array.shape[0],
    "width": dest_array.shape[1],
    "count": output_count,
    "dtype": dest_array.dtype,
    "crs": "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs ",
    "transform": out_transform
}
# open output file in write/read mode and fill with destination mosaick array
with rasterio.open(
    os.path.join(r'C:\Python\Raster\Mosaic', 'test.tif'),
    'w+',
    **dest_profile
) as mosaic_raster:
    for src in sources:
        for ji, src_window in src.block_windows(1):

           print(ji)
           r = src.read(1, window=src_window)
            # store raster nodata value
            nodata = src.nodatavals[0]
            # replace zeros with nan
            r[r == nodata] = np.nan
            # convert relative input window location to relative output
window location
           # using real world coordinates (bounds)
           src_bounds = windows.bounds(
               src_window, transform=src.profile["transform"])
            dst_window = windows.from_bounds(
                *src_bounds, transform=mosaic_raster.profile["transform"])

           # round the values of dest_window as they can be float
           dst_window = windows.Window(round(dst_window.col_off), round(
               dst_window.row_off), round(dst_window.width),
round(dst_window.height))
            # before writing the window, replace source nodata with dest
nodata as it can already have been written (e.g. another adjacent country)
            # https://stackoverflow.com/a/43590909/1979665
            dest_pre = mosaic_raster.read(1, window=dst_window)
            mask = (np.isnan(r))
            r_mod = np.copy(r)
            r_mod[mask] = dest_pre[mask]
            mosaic_raster.write(r_mod, 1, window=dst_window)

os.remove(memmap_file)


texts: qgis 2.8 openlayers plugin : google & OSM

texts:
I'm working currently for a big project with QGIS 2.8 (windows). I would like load satellite images (google or OSM for instance) as baselayer. I know the usefull plugin called openlayers plugin. I would like ask 2 questions please.

1/ Am i allowed to load and to generate maps with google earth images satellite (via qgis/openlayersplugin) in a professional context ?

2/ If it not the case,i could  maybe use openstreet map satellite base layer (free licence). But it seems to be impossible to load OSM satellite image via your plugin openlayers layer?

Is there a way to load Openstreetmap satellite image baselayer into Qgis?


texts: Differences between ExtJS and Dojo Framework

texts:
I have been using OpenLayers and the ESRI's JavaScript API and found that both using different set of JavaScript libraries have some implications on the quality of application that they deliver. Dojo has three dominions, including dijit and dojox but there is no such aspect in ExtJS.

What I have come to believe is that, this is the most prominent feature that distinguishes Dojo from ExtJS. This provides the opportunity to change the look and feel of Dojo based apps more easily as compared to ExtJS. The look and feel of ExtJS can be very easily noticed.

I am still at a pretty early stages of learning and so might have a wrong notion, but I would definitely appreciate if you could suggest as to what makes ExtJS framework different from Dojo Framework?

Also, I keep stumbling on these terminologies of Framework and Libraries, could you also clarify, when are they Framework and when are they Library or do both mean the same?


texts: How do I convert WKT geometry from 3D to 2D

texts:
I'm working with kml files using Python gdal.ogr module. When I extract feature's geometry using feature.GetGeometryRef(), I get:

POLYGON ((44.184069865521963 41.971859752507207 0,44.355248797385123 41.976261187945859 0, ...

However, I do not need this Z value - I need just a plain 2D WKT. After hours of searching I found FlattenTo2D() method, but feature.GetGeometryRef().FlattenTo2D() returns None

So how do I convert 3D WKT to 2D?


texts: QGIS not including complete symbology when exporting from print composer to PDF

texts:
I have a vector layer that is displayed with both cross-hatching and a border.  When I export to PDF, only part of the symbology is exported - the cross-hatching shows up without the border.  It shows up fine in the print composer window, and if I export to JPEG.

Screenshot of portion of print composer:

and (roughly) same view from PDF:

This problem is only showing up in 1 category of a layer displayed as catagorized values.  The other categories in the layer are set up somewhat similarly, but aren't running into the same issues.  The symbology for this category is composed of three layers, a simple fill for the border, and two line pattern fills for the cross hatching.  I've tried changing the order, but it doesn't appear to make a difference in the output.

This issue has been consistent with this particular layer, so it doesn't appear to be a memory issue. It still occurs if I create a new blank page in print composer and add just a map with only the offending layer.


texts: Relative paths not working in model builder

texts:
I checked "Store relative path names" in my model, but ArcGIS seems to ignore the settings - the paths are always the same as on my computer.
Plus, the "Store relative path names" unchecks automatically sometimes (I didn't find what triggers it).

How to ensure the relative paths are used? Why does ArcGIS ignore the "Store relative path names" option?


texts: Convert Esri ASCII grid raster to WKT polygons efficiently

texts:
I need to transform a large number of Esri ASCII grid files to polygonized WKT format, e.g. as CSV. Currently, I use GeoTools with the ArcGrid reader and a polygon extraction process. For parallelism this function is executed via Spark. Still converting 30-40 GB of the raster data is very slow. Do you see any  possibilities to improve performance?

Would you recommend GDAL / OGR over GeoTools for this purpose? Is there another library which is faster?

edit

Also as these are many small files I want to point out that it is related to https://stackoverflow.com/questions/44322914/spark-read-zip-file-which-contain-multiple-files-with-cpu-intensive-parsing where I try to parse the zip files in parallel.


texts: Google Earth Engine - NDVI time series that excludes two middle years

texts:
Please assist me to better capture temporal changes in NDVI using 2000 to 2017 years, where I want to exclude year 2002, and 2005 as these were drought years. If I exclude these will this affect my Mann Kendall analysis to check for monotonic trends?

var collection = ee.ImageCollection('MODIS/006/MOD13Q1')
                   .filterDate(ee.Date('2007-11-01'), ee.Date('20117-04-30'));


texts: Updating Oracle Spatial geometry column using SQL?

texts:
I am trying to update a table B geom column with the coordinates from table A.
Table A having coordinates in two separate columns x and y.
Table B to have the coordinates from Table A into Table B's geom spatial column.

Table A and B have an ID in common for each record.

I usually use FME to do this but this time I'd like to do it using pure SQL.

Any suggestions on how to do it?


texts: How can I find the segments that precede and follow a segment in a multiline

texts:
I have a multiline. Along it I have nodes which form intersections between the different segments of my multiline.

My goal is to inform each segment of the segment that precedes it and the segment that follows it in two columns of my multiline table.

I tried this query

select cm_code,
  (select cm_code from thd.t_cheminement n where st_intersects((p.geom), st_buffer(n.geom, 1))) cm1,
  (select cm_code from thd.t_cheminement n where st_intersects((p.geom), st_buffer(n.geom, 1))) cm2
from thd.t_cheminement p;


texts: Qgis standalone script using eclipse IDE -No data provider available issue

texts:
I know, the question is asked few times also in this forum but
i apply already the presented solutions without any chance to fix the problem.
In General the source code looks like the following:

from qgis.core import *

# supply path to qgis install location
QgsApplication.setPrefixPath("C:/path/to the/apps/qgis", True)
print QgsApplication.showSettings()

# create a reference to the QgsApplication
qgs = QgsApplication([], False)

# load providers
qgs.initQgis()

if len(QgsProviderRegistry.instance().providerList()) == 0:
    raise RuntimeError('No data providers available.')

mylayer = QgsVectorLayer("C:/path to/the actual /layer.shp","DB_tracks", "ogr")
if not mylayer.isValid():
print "Layer failed to load!"
# When your script is complete, call exitQgis() to remove the provider and
# layer registries from memory
qgs.exitQgis()

Now the actual issue is the "no data provider available" exception
thrown by the QgsProviderRegistry.
By researching  the directory to the providers seems set correct.
The directories are set according
to the respond of QgsApplication.showSettings() as shown below:

Application state:
QGIS_PREFIX_PATH env var:
Prefix:     C:/NONBKU/QGIS Essen/apps/qgis
Plugin Path:        C:/NONBKU/QGIS Essen/apps/qgis/plugins
Package Data Path:  C:/NONBKU/QGIS Essen/apps/qgis/.
Active Theme Name:
Active Theme Path:  /themes\\icons/
Default Theme Path: :/images/themes/default/
SVG Search Paths:   D:/Users/Username
    C:/Path/QGIS Essen/apps/qgis/./svg/
User DB Path:   C:/Path/QGIS Essen/apps/qgis/./resources/qgis.db

By calling than QgsProviderRegistry.instance().providerList())
the following exception is thrown:

Traceback (most recent call last):
  File "C:\path\workspace_eclipse\DB_QGIS_Python\src\python_DB_test.py",
line 18, in <module>
raise RuntimeError('No data providers available.')
RuntimeError: No data providers available.

So any ideas which issue forces the "no data providers available"
exception?


texts:
How can I convert an OSM file to Shapefile? Is there a good tool or a good web page.

I tried this link https://code.google.com/archive/p/osm2shp/, but there is an error.


texts: Setting background RGB value to invisible for multiple rasters in ArcGIS Desktop?

texts:
I'm using ArcGIS Desktop 10.0.

I need to set my background rgb value of 255,255,255 to transparent. In the layer-Properties Windows > Symbology I can do it for one layer, but can I batch set this, or use arcpy? I have to do it for around 300 layers.


texts: How can i set layer polygon in center in QGIS using python?

texts:
I am trying to set a fixed zoom level whenever I layer is shifting on  different polygon. I have coordinates for all the polygons , I need each polygon to be on center with certain zoom level so that programme can take printout of the layout automatically.

PS : For print functionality , I already have code. I have made print layout program but it's make image centered polygon images with zoom level.

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.utils import *

mapRenderer = iface.mapCanvas().mapRenderer()
c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)

x, y = 0, 0
w, h = c.paperWidth(), c.paperHeight()
composerMap = QgsComposerMap(c, x,y,w,h)
print composerMap
c.addItem(composerMap)

item = QgsComposerScaleBar(c)
item.setStyle('Numeric') # optionally modify the style
item.setComposerMap(composerMap)
item.applyDefaultSize()
c.addItem(item)

dpi = c.printResolution()
print dpi
dpmm = dpi / 24.5
width = int(dpmm * c.paperWidth())
height = int(dpmm * c.paperHeight())

# create output image and initialize it
image = QImage(QSize(width, height), QImage.Format_ARGB32)
image.setDotsPerMeterX(dpmm * 100)
image.setDotsPerMeterY(dpmm * 100)
image.fill(0)

# render the composition
imagePainter = QPainter(image)
sourceArea = QRectF(0, 0, c.paperWidth(), c.paperHeight())
targetArea = QRectF(0, 0, width, height)
c.render(imagePainter, targetArea, sourceArea)
imagePainter.end()

image.save("C:/example/Ashwin/out2223.png", "png")


texts: Symbolizing Overlapping Oil & Gas Polygons in ArcGIS Desktop?

texts:
I am mapping an area where there are a lot of overlapping features because of certain Oil & Gas wells being drilled within the same geographic area, but to different depths and maintained by different companies.

There might be one large Oil & Gas unit (polygon) that is drilled to 5,000' and another smaller unit (polygon) within the same geographic boundary that is drilled to 10,000'.

Each polygon is symbolized based on company, so I have been messing with symbol levels to draw certain ones on top, but I don't feel like this is the best way to get ArcMap to automatically draw smaller ones on top and keep the larger ones on bottom. Is this possible?


texts: How to update only non-empty fields with a user input value in a model?

texts:
I have a model which iterates through a GDB, populating the same attribute field in every feature class with a parameter value.  In this way I can mark all FC's in the geodatabase as to what project they are a part of.

However, sometimes a GDB contains FC's which are migrated from other GDB projects, with the attribute already populated with a value designating it as part of the other project.  Currently my model stomps on this value and overwrites it with the new project code, whereas I want to retain the existing information.

I would like to have the model evaluate the field for the record it is working with, and if the field is not empty, to leave the field alone and move to the next record.  However if the field is  then it should write the new project code to the field and move on.

I've looked into a bit of the if/then questions and answers but have not found a simple way to detect if a field is empty or not, and if not, how to skip it?

I am now trying to write some Python to handle it and have the following, however I am not sure how to get the parameter value from the user entry in the model (ie. UserInputFromModelParameter), or what command to send to make the script skip to the next record if the script is true.

# set up variables
def BlockName:

def IgnoreField:

# set block name from user inpuit
BlockName = [UserInputFromModelParameter]

# set values to watch for
IgnoreField = [NULL, ""]

# evaluate the field

    if [BlockID] == IgnoreField
            MoveToNextRecord();

    else:
        [BlockID] = BlockName();


texts: Storing geographic coordinates for multivariate statistical analysis?

texts:
In terms of Principal Components Analysis or Multiple Discriminant Analysis, I have used dummy variables to represent a group of like-features that are located in the same geographic region. For example, in the analysis of two metropolitan cities, the census tracts that make up city one are classified as 0 and city two as 1. This binary classification groups the census tracts into two regions.

I am wondering how I could store geographic coordinates (centroid?) of these census tracts so that I could introduce a spatial component to the multivariate analysis. This is dissimilar to the approach mentioned above since I am not grouping any features; but rather, I am interested in determining whether or not spatiality is a statistically significant independent variable (in terms of regression), or if spatiality is a significant part of a component/factor (in terms of PCA/FA).

This would be simple enough if the variable was a single measured value; but alas, we have x and y coordinates to think about.


texts: Using Plugin in standalone PyQGIS script

texts:
Typically, when I want to use a QGIS processing tool for my PyQGIS standalone script I write something like:

processing.run('buffer:...'{…})

Now, I would like to use the functionality of the QGIS Plugin NNJoin in my script. Is that possible and if yes how is it done?

I'm using QGIS 3.4

ADDITION:
I found the following documentation for a processing tool in 3.10.
join Attributes by nearest

If I implement this in my standalone script, it won't work in my QGIS 3.4  right?


texts: Making overview with linked detailed maps in QGIS?

texts:
i want to create 4-5 detailed maps, and then an overview map that has rectangles that markes the locations of the detailed maps. Best case would be if the rectangles were links to the detailed map.

I have managed to make the overview and the rectangles, and think I know how to make the detailed maps one by one. The rectangles are drawn manually on the overview map.
All layers that I plan on using are in the same QGIS project.

The final product that I hope to create would be the overview map on the first page, followed by the detailed maps. The overview would have links to the detailed maps, click on the rectangles to open the right detailed map.

Is there a way to link the detailed maps to the overview in QGIS 2.8 print composer?

I'm still not sure if I'm doing something wrong, when I try to create the overview. I want it to show the extents of all of my detailed maps.

Is this possible?


texts: Generating DSM and DEM from processed LiDAR (TIFF) on ArcGIS Pro

texts:
I have a "processed" LiDAR that is available as TIFF, I do not know how it was processed although I've read about the method of converting a LAS dataset to Mosaic, terrain, etc.
I was basically just given the tiff below to generate a DSM and DEM, it is basically just flat.
I have searched on how to create dem and dsm surfaces from the lidar raster but can't find anything about it, although there's information of accomplishing that if it were a point clouds.
So how can I possibly create DSM and DEM from the Lidar tiff?


texts: Arc is changing decimals!! My scale is set to 2, why do my numbers in the select by attributes-unique values have more decimals?

texts:
My scale is set to 2.  Yet arc (using 10.2) is showing (in my select by attribute - get unique values) my number of 26.67  as 26.670000000000002.   My value of 25.25 is shown correctly as 25.25 in the same field.

Again my value of 32.60 is shown as 32.6000000000001.

I have tried replacing the values with the field calculator, but it goes back to these long decimal numbers all on its own.

Why is it doing this and how do I fix it?

These values are not more than 2 decimals, showing more is wrong. but it has to be a double because there are some negatives ..


texts: Attempted to load an infinite number of tiles while using proj4leaflet and EPSG:32614

texts:
I get error "Attempted to load an infinite number of tiles" while using proj4leaflet and EPSG:32614.

 crs = new L.Proj.CRS('EPSG:32614',
  '+proj=utm +zone=14N +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs', {
      resolutions: [
          8192, 4096, 2048, 1024, 512, 256, 128
      ],
      origin: [329370, 3970950],
      bounds: L.bounds([218128.7031, 6126002.9379], [1083427.2970, 7692850.9468])
  });

  options = {
     crs: this.crs,
    layers: [
      tileLayer('https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
        { attribution: '...' })
      ],
      zoom: 16,
      center: [329370, 3970950],
      attributionControl: false,
  };


texts: How can I find the percentage of representation of classes?

texts: I have done a classification using the semi-automatic classification plugin, and now, I would like to express each class in percentage of the study area (e.g. forest represent 20% of the study area and so one). Does anyone know how to get that ?


texts: Displaying directional bandwidths in QGIS

texts: I have a traffic model and want to display the results in QGIS. I have a shapefile for the road network and can easily annotate the flows onto this. However, I want to display the information as bandwidths, separate by direction. I have tried this but the bandwidths for the two directions overlap.


texts: Center text value of a field in polygon and display

texts:
I would like to display the value of a field over the polygon.

In this case, the field name is Subname,
and i would like to use a specific font.

Would like to perform this action for each polygon on the layer.


texts: Define color for label-buffer with data-defined override/expression

texts:
In QGIS 3.8 (Win 10) I have a Point-Layer with labels in a certain distance and leader lines, created with geometry-generator. The leader lines in some cases cross below the label and are masked by a text-buffer in white. However, some of the labels are above the sea - so in these cases, it would be better not to have a white buffer. I tried to make the buffer in blue by defining the color with data-defined override and using the following expression:

if (("NAME" = 'name 2a' or  "NAME" = 'name 3' or "NAME" = 'name 3a' ),blue, white)

However, the color does not change. I also tried to add the color-code in hex (#65baf2 instead o blue and #ffffff instead of white), but with the same result. How can I set different buffer-style for certain labels?

See screenshot:


texts: multipatch features are not visible in 3D ArcGlobe

texts:
I have created some 3D features in ROCKWORKS software, which represents lithology data. Then I converted it to 3D shapefile, and then I added it to ArcGlobe as well as ArcScene. It is clearly visible in 3D in ArcScene but in ArcGLobe it is not visible as 3D.  Instead it is just visible as a 2D Feature class.  Then later I have tried some extrusion options, but that graphic rendering is not good like ArcScene, and I have tried some other options like some sketch up files I have converted to multipatches but the same problem .  Is there any way to add 3D symbology from ArcScene to ArcGlobe?  Actually, I want to show lithology data like in the following link, and I want to create one globe service, how to do it.....

http://blogs.esri.com/Info/blogs/arcgisexplorerblog/archive/2011/04/15/3d-subsurface-profiles.aspx

Thanks and regards
Srinivas


texts: Converting GIS features into multiple CAD layers based on attribute values

texts: I am often requested to provide my water or sewer data in in a dwg or dxf format.  I can easily move data back and forth between the two platforms using either the native import or export functions ArcGIS (Info level 9.3.1 or 10) or AutoCAD (Map 3D 2010).  I am looking for a way to sort the features of a given featureclass into different CAD Layers based on their attribute values.  For example if my water pipe data has an attribute values for both Diameter and Material. On export I  would like to all of the 8 inch Ductile Iron pipes to be placed in the same CAD layer and all of the 10 PVC pipes in another CAD layer.  I have done this in the past by manually selecting the pipe similar pipes and exporting just the selections but this can be tedious and time consuming.  I have tried looking for a tool to do this, as I seem it is somethings others would be interested in, with no luck.


texts: Separating ("high enough") trees from "lower height" vegetation? In abstract (non-RGB) .xyz point cloud

texts:
Suppose I want to select trees that are "high enough", but the data contains other vegetation, like bushes and smaller trees.

My data is abstract (non-RGB) .xyz point cloud. I.e. there's only coordinates, no colors etc.

Consider as an example a "North European" forest, so e.g. like:

Forest Example Picture

How do I approach it?

Naive idea:

filter based on height/Z-value, i.e. discard points that are "below" some value.

Problem: If I filter and discover what points belong to "high trees", then how can I recover the "lower part" of those trees? To e.g. know the height of the tree?

Other ideas:

Classify in "3D k-nearest" sense to first group points that belong to a certain plant. This could allow me have point classes: bushes, trees, ... And I could remove point classes.

Problem: To know, which kind of classification works "consistently" for different inputs.

Some more advanced clustering algo, but ATM I'm not sure, which "view on clustering" would be most suitable for differentiating vegetation. E.g. density-based or something else.
A combination of different approaches. E.g.

Use naive "filter by height" to find some points belonging to higher trees.
Find some of these points from the cloud prior to the naive filter and do some clustering using these points as "anchor points". So that the clustering algo wouldn't have to care about other points than those related to the higher trees.

I'm currently using Python and pdal library. I've tried filters.cluster with varying parameters and it seems "decent", but it's not a very detailed implementation. I have hard time understanding what clusters to retain and what not.

This experimenting has suggested that I might also try doing several runs of the clustering, each getting rid of something, until my data is "clean" enough to cluster individual trees even.

I'm open to any tool, but this question is geared more towards identifying "theoretically suitable" approaches. Rather than, implementation-specific.


texts: Call ArcGIS Repair Geometry Function from Oracle 11.2

texts:
Could anyone please point me to the right direction on how to call an ArcGIS Function (i.e. Repair Geometry) from PL/SQL? All information I found so far concerns the other way round.

This is needed because Oracle's validation criterias differ from the ones used by ArcGIS causing some geometries produced within Oracle not drawn by ArcGIS services (see http://support.esri.com/en/knowledgebase/techarticles/detail/36047).

The geometries are stored directly in Oracle.

I thought about calling arcpy directly PL/SQL but couldn't find any examples for that.
Is there a way without using ArcSDE?

Any thoughts on this would be very much appreciated.


texts: How to save the edited features on a layer back to PostGIS table?

texts:
I have PostGIS + OpenLayers + GeoServer. I have developed an app which can show layers from PostGIS DB, allow user to edit, add, view features on a vector layer. Now, my need is to save the edited or newly created features back to PostGIS DB.

Is it possible to achieve it using only OpenLayers??


texts: How to solve this problem with proportional size data in qgis print composer

texts:
I would like to export a map with proportional circle representation for the number of inhabitant of each city.
I've used the assistant size of QGIS and created a proportional legend within qgis options.

Then when in the print composer everything is fine in the legend settings panel but within the frame it's crappy as you can see.

What should I do?

QGIS 3.10.1


texts: GeoServer sub-domain, where to place a HTML file

texts:
I would like to host a website, which is addressable as a sub-domain of my GeoServer. For example http://mygeoserver.com:8080/geoserver/web/my_website.html.

My question is: where do I have to place the html file in my folder structure so that I can access and run it with the http request above? How can I make the web page available in the web directory of GeoServer?

The GeoServer is installed under Ubuntu.

Edit (1):

Thanks for your comment BradHards.

Background:
The background to my question is that I want to have an html file with an OpenLayers WMS map, which shows my GeoServer data. In this map I need WFS as an overlay to the WMS map. Up to now, I can run the HTML file with WMS fine, but I can't get the WFS overlay to work. In this forum thread (WFS layer not showing?) someone mentioned a hint according to which the OpenLayers/WMS/WFS html file must be located within the GeoServer web directory in order to get it to work.

My steps:
First I installed Tomcat6 under Ubuntu 12.04 LTS (sudo apt-get install tomcat6).
After that, I installed the GeoServer version 2.2-SNAPSHOT (apt-get install opengeo-geoserver). It is running well.

My tries:
I tried a bit for myself and I am close to a solution.
In the folder "var/lib/tomcat6/webapps/geoserver/data/www" are two HTML demo files, ol-demo.html as demo for OpenLayers and wfs-t.html as demo for WFS-T. I can successfully call these files in the browser (localhost and from remote) by (http://myaddress.com/geoserver/www/ol-demo.html). But when I additionally place my HTML file in this folder I can't reach them from the browser in an analogous way. Instead: HTTP 404 the requested resource() is not available.

Hope this helps to be able to better answer my question.


texts: So basically, I've written a script that generates a lot of polygons over the UK with certain criteria. The output it a lot of shapes like the following - I wanted to know if, in QGIS, I'm able to join shapes like this together and fill in the holes by making a big polygon over the whole thing.


texts: Creating homogeneous groups of areas from raster using QGIS

texts:
I'm trying to create homogeneous groups of areas from a raster which contain elevations data (in the red area).

currently I am doing it manually creating new polygons for each groups.

Is there a way to do it automatically using QGIS?


texts: Does Google Maps API or OSM provide location/place boundaries

texts: I've been digging around today in Google Maps and Leaflet API's to find if there is an easy way to get these boundaries instead of needing our users to manually outline these things (like drawing polygons on the map)?


texts: Printing at defined scale and paper size from ArcMap?

texts: How to print a drawing to a defined scale (e.g. 1:2500) in ArcMap, and to a given paper size (say A3)?


texts: PyQGIS QgsVectorLayer() Loading Invalid Layer in Standalone Python Script?

texts:
In my standalone python 3 script, QgsVectorLayer() is loading an invalid layer. When I use the exact same function and inputs in the QGIS GUI python console, the layer loads fine. I am not sure what is missing in my standalone script. I have double checked the paths and made sure they are correct. I used

QgsApplication.prefixPath()

to check the correct path for input in the

QgsApplication.setPrefixPath()

function within my standalone script. The path and inputs I use in the

QgsVectorLayer()

function in my standalone script are the same as those used in the GUI python console. I'm not sure why loading a vector layer in my standalone script it failing. The vector object is created, but .isValid() returns False.

Here is my standalone script:

import sys, os, time

sys.path.extend([r'C:\OSGeo4W\apps\qgis\python',r'C:\OSGeo4W\apps\Python37\Lib\site-packages'])

#modify environment variables to find qgis and qt plugins during qgis.core import
os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = r'C:\OSGeo4W\apps\Qt5\plugins'
os.environ['QT_PLUGIN_PATH'] = r'%QT_PLUGIN_PATH%;C:\OSGeo4W\apps\Qt5\plugins;C:\OSGeo4W\apps\qgis\qtplugins;C:\OSGeo4W\apps\qgis\plugins'
os.environ['PATH'] += r';C:\OSGeo4W\apps\qgis\bin;C:\OSGeo4W\apps\Qt5\bin;C:\OSGeo4W\\bin'

from qgis.core import *
from qgis.gui import *

# supply path to qgis install location
QgsApplication.setPrefixPath(r'C:\OSGeo4W\apps\qgis', True)
#QgsApplication.setPluginPath('C:\\OSGeo4W\\apps\Qt5\\plugins\\platforms')
#print(QgsApplication.systemEnvVars())

# create a reference to the QgsApplication
# setting the second argument to True enables the GUI, which we need to do
# since this is a custom application
qgs = QgsApplication([], True)

# load providers
qgs.initQgis()

##########################
# Write your code here to load some layers, use processing algorithms, etc.
canvas = QgsMapCanvas()
canvas.show()
layer = QgsVectorLayer(r'C:\Users\Matt\OneDrive\FarmProject\Kankakee_Parcels\K3_TaxParcels.shp', 'Kankakee', 'ogr')
if not layer.isValid():
  print('Failed to open the layer')

# add layer to the registry
add_layers = QgsProject.instance().addMapLayer(layer)

# set extent to the extent of our layer
canvas.setExtent(layer.extent())

# set the map canvas layer set
canvas.setLayers([add_layers])
canvas.refresh()
time.sleep(30)
########################

# When your script is complete, call exitQgis() to remove the provider and
# layer registries from memory
qgs.exitQgis()


texts: Cause of expected an indented block error in Python script?

texts:
I have been trying to list all the field names of a shapefile but it is showing 'expected an indented block'.

What should I do?

import arcpy

fieldList = arcpy.ListFields("C:\Users\Mridul\Documents\GIS\Spring 2011\650\streets_dev.shp")

for field in fieldList:
print field.name + " is a type of " + field.type + " with a legth of " + field.length

rows = arcpy.UpdateCursor("C:\Users\Mridul\Documents\GIS\Spring 2011\650\streets_dev.shp", "SHAPE_len>400")

cnt = 0

for row in rows:
    row.STREETNAME = "OLD " + row.STREETNAME
    rows.updateRow(row)
    cnt +=  1
    print cnt

del row
del rows


texts: Creating Layouts for selected windows in MapInfo

texts:
I have a workspace from which I need to create a lot of different layouts, which need to be kept. Each time I create a layout which has five windows (2 of which are the same, copyright and key), I then minimize them all and open a new map window for the next site.

When I then open a new layout window for all currently open windows, every window appears in the layout and I have to delete all but the ones that I need.

Is there a way of specifying which ones I want to be included each time?


texts:
My UpdateCursor isn't working in arcpy.  Here's my code:

urows = arcpy.UpdateCursor(firstDatabase+ os.sep + "Temp")
for urow in urows:
    urow.Import_Date = time.strftime("%d/%m/%Y")
    urow.File_Name = str(table)
    urows.updateRow(urow)

It will always error out on the Import_Date field saying "error executing function".

The Import_Date is a 'date' field.

Any ideas?

Here is a screen shot of the loop when I add "arcpy.AddMessage(urow)" in the code and # out the other stuff.


texts:
There may be something I'm missing, but is it possible to use the arcpy module outside arcmap?

My situation is that we have one license on a computer which is shared over a network.  I need to do some basic processing and I'd like to do it via python without having to log on to the GIS computer.

Is this possible?

Is it as simple as copying the location of the arcpy module into another instance of python?

Or, can I just direct python to the directory arcpy is living in?


texts: GeoServer To Vue2-leaflet projection issue

texts:
I am using GeoServer (using Tomcat in localhost) and to send a WMS file to my Vue2-leaflet application.

My projection in GeoServer is ESPG 4326

This is my Leaflet code :

<template>

WMS Layers with Layer selector

        NOTE: Use layer symbol on the top right to select which layer you
        want to display. More information about WMS (Web Map Service) can be
        found on the
        leaflet.js WMS example page.

    <l-map
      :zoom="zoom"
      :center="center"
      style="height: 80%"
    >
      <l-control-layers />
      <l-wms-tile-layer
        v-for="layer in layers"
        :key="layer.name"
        :base-url="baseUrl"
        :layers="layer.layers"
        :visible="layer.visible"
        :name="layer.name"
        layer-type="base"
      />
    </l-map>

</template>

It seems that I have a projection problem when I am rendering into my browser:

How I can fix this?


texts:
I'm new to OpenLayers,currently in existing project they have created markers in openlayers map by using following code,

// vector layer creation
vLayer = new OpenLayers.Layer.Vector(layerName, {
        projection: map.displayProjection,
        extractAttributes: true,
        styleMap: new OpenLayers.StyleMap({'default' : MarkerStyle}),
        strategies: [new OpenLayers.Strategy.Fixed()],
        protocol: new OpenLayers.Protocol.HTTP({
            url: listurl,
            format: new OpenLayers.Format.GeoJSON()
        })
    });

map.addLayer(vLayer); // Adding the vector layer to the map

To display markers,coordinates are fetching from the url, but now we are having coordinates information in a hidden block of HTML like as follows

    <input type="hidden" id="latitude_1" value="17.43189">
    <input type="hidden" id="longitude_1" value="78.447167">

    <input type="hidden" id="latitude_2" value="15.43189">
    <input type="hidden" id="longitude_2" value="76.447167">

Now instead of fetching coordinate information from url,I want to use above information(which is in hidden block of HTML) to show markers,how can we modify existing code ?
I've chance to modify HTML code according our convenience.


texts: Using ModelBuilder to Summarize Multiple Feature Classes?

texts: I am trying to create a model that will iterate through all the files in a Database, return the counts of each feature class, then add them all together to get a final count and put the result in a table. I have created an iterator that goes through the database and summarizes each feature class, but now I want to get a total count of the features.


texts: Export ArcMap cartographic representation to a symbol?

texts:
Using Cartographic Representation (CR) I've created a nice feature symbol with a filled interior and sharp outline. It looks great on a map, on screen and printed. However I can't duplicate that symbol in the legend, as the regular Symbol Property Editor doesn't have the same level of control and edit options as the CR marker editor does.

Here's the CR marker editor, and the resultant symbol:

And here's the closest I can get with regular symbol editor, just not good or close enough. With stacked level character symbol the two don't line up at all, and with halo effect the effect is thicker on the right side than the left.

So, how can I export the CR symbol so I can use it in the legend (and elsewhere)?

We're not using live legend as it doesn't afford enough control over the result, for instance the size of the CR symbol can't be changed, even after converting to graphics.


texts: ProgressBar for a function that has to communicate with QObjects

texts:
My Plugin uploads Layers to a WebGIS and the QGIS GUI turns inactive during this process. So that the User knows that the Program hasn't crashed I want to implement a ProgressBar.
My first try gave me a black Window that only refreshes after the task was done. So I learned that I have to use threading. But I'm not sure if I should use threading provided by Python, PyQt or PyQGIS.
Also like the QGIS Documentations says I can't use data from the main task, so I'm not sure how to continue from here. All the tutorials I found, only build threads with time.sleep but I guess I need another depth of complexity.


texts: Is there a 3rd party ArcGIS viewer (builder) for Javascript?

texts:
Esri provides an ArcGIS Viewer for Flex as well as an ArcGIS Viewer for Silverlight.

Strangely though, there are no plans for a similar viewer for Javascript.

Are there any 3rd parties that have created a "viewer" for javascript?

Though they are referred to as "viewers", they are actually viewer-builders: they allow administrators to create web applications without writing any code. The javascript viewer blogged about here does not allow user to choose tools to include in the app, as can be done with the flex and silverlight versions.

Update

I have developed some capabilities that have been packaged as a combination of an SOE and an Add-in for the silverlight viewer.  This works great - when the admin builds a viewer application, he can pick my add-in tool from the list and configure it to use the URL of a server where my SOE is installed.  However, from what I can tell there are no add-ins for the javascript viewer.  Or am I missing something?


texts:
I want to warp a vector layer (set of polygons, but not a raster map) to fit a map. The problem is, the vector layer coordinates are from 0 to 360 and the basemap layer coordinates are from -180 to 180

So I tried subtracting -360 from polygon coordinates, but polygons are broken (maybe because of meridian)

Eventually what I want is this picture, (just for description, I duplicated the layer). Can somebody help please?

The polygon file is GeoJSON.


texts: QGIS Cross Section Profile Distance using Ellipsoid?

texts:
I am working with QGIS 1.8 and 2.6 and can't get any of the Cross Section Profile tools to provide the ellipsoid ("proper") distances. All the distances use geographic coordinates which are too long/far.

For example, the straight line distance from Vancouver BC, to Kelowna is about 270km. My cross sections all say it is 420km.

The QGIS distance measurement tool allows you to use the Ellipsoid (I set mine to WGS 84) and give me 271km.

Can anyone shed light on how to get any of the profile tools to give me the "proper" distance?


texts: Python program for "Comparing total number of pages of an mxd to total number of PDFs created from that mxd"!

texts:
I need some help on my program. I have the parameters set to ask user to pick an mxd and the folder where those data driven PDFs are located. I started with Search Cursor to create a list from the Sheet_ID field of an mxd. To get total number of pages of that mxd. Now I need to write in for the program, to compare the list created from the Sheet numbers to the number of PDFs created and tell me if there is a difference in the total count. For Example: If there are 1 out of 11 sheets for the mxd, then there should be 11 PDFs created and so on. MY code is below:

#Purpose: Compares SheeID field to the toal PDFs created from Data Driven Pages Program
Import arcpy, os, string, sys

#Make parameters for people to choose mxd and folder of PDFs to compare
mxdList = string.split(arcpy.GetParameterAsText(0), ";")
dir = arcpy.GetParameterAsText(1)

#Use Search Cursor to go through Attribuite Table to get Sheet number info
rows = arcpy.SearchCursor(mxdList, "Sheet_ID")

#Use list from Sheet_ID field to match the number of PDFs created to get difference


texts: 3D Interpolation using point cloud in GRASS GIS

texts:
I've imported a set of 3D points (z,y,z) into GRASS GIS. I'm having trouble: a) visualising my dataset in 3D using nviz or the standard data viewer (I can only get a 2D plan view of my data) and, b) interpolating them into a solid 3D voxel object. My dataset is pretty straightforward and there's only 1520 points.

The GRASS user interface is not that user friendly and I can't find a suitable tutorial to follow. Any ideas or pointers would be very much appreciated. I'm sure I should be using the v.vol.rst tool, but when I input my parameters, I get the following error message:
ERROR: Column type of wcolumn is not supported (must be integer or double).

Any ideas?


texts: Using a postcode sector map with a population attribute, create clusters of those sectors to have a total population within a defined range

texts:
I have a postcode sector map in QGIS and one of the attributes is 'population'. I am trying to create a map showing what postcode sectors would need to be clustered together in order to split the whole country up into groups of sectors so that each group had a total population within a specific minimum and maximum range. I don't really care about the postcodes that fall between the clusters, they can be treated as noise for the purpose of this part of my task. I would prefer the clusters to be blobs rather than worms of sectors joined together weaving across the map. I have tried all the different types of clustering I can find and various plugins but none do the job, or i don't understand them enough.  Can anyone suggest where I should be looking or what I should be trying to find? (not knowing the correct terminology isn't helping at this point either).
I've tried. K-Means, Density Based Clustering, Group-average agglomerative clustering and many others but none seem to be easily implementable in QGIS.


texts: Geotools JTS, find most distant points in polyline or polygon?

texts:
is there any convenient API in JTS that would help to

find two most distant points in line
find two most distant points in polygon


texts:
As the title states I'm trying to open a REST/JSON service in QGIS, but without luck so far. Now I'm wondering if this kind of service is even possible to open in QGIS with current, existing tools. I'm using the new QGIS 3.0, but I tried with the latest 2.x.x version too.

Here's an example of the URL:
https://tie.digitraffic.fi/api/v1/data/weather-data
(open in browser in order to view JSON tree)

As far as I know QGIS is capable of opening REST services with add layer --> add ArcGIS featureserver layer, but that hasn't worked in this case. Then on the other hand I know ArcGIS has a specific JSON-toolset which might imply that this wouldn't be an issue with ArcGIS.

EDIT: As you can see in the tree, there's no geometry or coordinate information. Presumably that will be joined later based on id.

Is there anyone with more knowledge about the matter?


texts: How to write a python code for Add Feature Class To Terrain using interface?

texts:
I am trying to write a code that creates a Terrain from lidar data (LAS). The problem is that the inputs are in Feature Dataset in GDB but I don't know how to write code so that using an interface. I won't write the whole code just the part that it's not working. Las data is saved in the new Feature Dataset created in the same script. I need to somehow replace D:\Google Drive\GDB2.gdb\data2\las with just las. I tried everything - by removing '' and it doesn't work. The error that I always get is: executing tool - Failed to execute (GDB). Can anyone help me? I wrote the parts of the code that are significant for creating terrain, the last line doesn't work:

import arcpy
GDB_Location = arcpy.GetParameterAsText(0)
GDB_name = arcpy.GetParameterAsText(1)
GDB_file = arcpy.CreateFileGDB_management(GDB_Location, GDB_name)
out_dataset_path = GDB_file
out_dataset_name = arcpy.GetParameterAsText(2)
feature_dataset = arcpy.CreateFeatureDataset_management(out_dataset_path,out_dataset_name)

arcpy.DefineProjection_management(feature_dataset, "PROJCS['WGS_1984_World_Mercator',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],UNIT['Meter',1.0]]")

....

l_g = 'las'
las_data = '{0}\{1}'.format(feature_dataset,l_g)
arcpy.LASToMultipoint_3d #in this step in created las

fd = feature_dataset
terrain = 'terrain'
terrain2 = '{0}\{1}'.format(feature_dataset,terrain)

    # Process: Create Terrain
arcpy.CreateTerrain_3d(fd,lidar_terrain, Average_Point_Spacing, "50000", "", "WINDOWSIZE", "ZMIN", "NONE", "1")

# Process: Add Terrain Pyramid Level
arcpy.AddTerrainPyramidLevel_3d(terrain2, "", "'1 5000'")

# Process: Add Feature Class to Terrain
arcpy.AddFeatureClassToTerrain_3d(terrain2, "'D:\Google Drive\GDB2.gdb\data2\las' Shape Mass_Points 1 0 1 true false las_ground_embed <None> false")


texts: Retrieving CSV from GitHub Repository for later use in ArcGIS Online?

texts:
I am trying to use data from GitHub at this location: https://github.com/factbook/factbook.csv/blob/master/data/c2147.csv

Is there a way to convert this to csv or utilize in some other more interesting way?

In the end, I would like to publish the resulting layer to ArcGIS Online.  I have the country boundary file to join already.


texts: Dynamically create leaflet popup via JavaScript Object

texts:
I want to dynamically populate a marker popup via a JavaScript object.

I have a JavaScript object like so:

{address: "92101 S Martin Ave.", address2: "", agent_name: "mr smith",
animal: "sheep", car: "honda", city: "Chicago", county: "Cook County", ID: "1",
latitude: "41.718986", longitude: "-87.550646", postalcode: "60617",
stateProvince: "IL", status: "GS"}

I want to be able to loop thru all of properties within the object and display those within the popup.

I have tried the below:

for (var j = 0; j < Object.keys(data).length; j++) {
                var marker = L.marker([data[j].latitude, data[j].longitude]).addTo(markers);
                marker.bindPopup(data[j]);

This does not work.


texts: Removing empty polygon from sf object in R?

texts:
I have such object:

> p_mb
Simple feature collection with 52 features and 0 fields (with 9 geometries empty)
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 585001.4 ymin: 211001.9 xmax: 609998.1 ymax: 236998.6
epsg (SRID):    NA
proj4string:    NA
First 10 features:
                         geometry
1  POLYGON ((591051.4 211855, ...
2                   POLYGON EMPTY
3  POLYGON ((588206.4 212570, ...
4  MULTIPOLYGON (((585004.6 21...
5  POLYGON ((600614.5 216112.9...
6  POLYGON ((603033.8 215125.2...
7  MULTIPOLYGON (((608168.4 21...
8  POLYGON ((600009.4 212984, ...
9  POLYGON ((603203.4 217336, ...
10 POLYGON ((585519.4 217546, ...

obtained from negative buffering other sf object. My question is simple, how to delete those empty slots to not get error and warnings:

> p_pol <- sf::st_cast(p_mb0, "POLYGON")
There were 23 warnings (use warnings() to see them)
> warnings()
1: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
2: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
3: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
4: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
5: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
6: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
7: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
8: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
9: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
10: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
11: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
12: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
13: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
14: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
15: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
16: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
17: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
18: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
19: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
20: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
21: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
22: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only
23: In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only

> spl_mbuf <- as(p_pol, "Spatial")
Error: length(srl) > 0 is not TRUE

I think I can skip those warnings.

Found workaround: just used p_uni <- sf::st_union(p_mb) after buffer, but this is only emergency exit.


texts: Does anyone know if there are any applications can support exporting Shapefile to GeoJSON and keeping the style?


texts: Transparency on black/white raster after clipping in QGIS?

texts:
I'm new to QGIS.  It seems as if this problem doesn't show up with ArcGIS.

I clipped a black/white raster. Unfortunately the surrounding background after clipping is also black. When I set the transparency to 0, black areas within my raster get white. I'ld like to avoid that. So here's what's happening in pictures:

So one workaround could be to just overlay the vector layer and raster, then it's perfect:

But I need this layer as only one raster or vector and not only as an overlay in QGIS.

So how can I save that?

Converting the raster into a vector or vector into a raster for some reason doesn't work as I would expect it to do.


texts: Extract number of vertices in each polygon in R

texts:
Simple question: how to extract the number of vertices in each polygon in R?

I've found one way that works, but this is very slow (when x is a SpatialPolygonsDataFrame):

y <- fortify(x)
n.verts <- aggregate(y, by=list(y$id), FUN=length)[,2]
x@data$nverts <- n.verts


texts: Stuck on how to list loaded layers in QGIS 3 via Python

texts:
I've stuck on a issue about getting a list of loaded layers. I've two layers loaded, one spatial (point, line or area) and another one non spatial (postgresql table).
I've tried this code but it lists only the spatial layers:

layername = []
for i in iface.mapCanvas().layers():
    layername.append(i.name())
print(layername)

I've tried QgsProject.instance().layerTreeRoot()

root = QgsProject.instance().layerTreeRoot()
ids = root.findLayerIds()
print(ids)

but it gives me a layer name plus an id.
How can I get a list of layer names from all loaded layers (spatial and non-spatial) in QGIS 3?


texts:
I want to use gdal_calc.py to reclassify a raster using numpy's conditional statements, but I am running into an error message. It appears that gdal_calc.py cannot access the numpy conditional functions. Does anyone know how to link numpy in the gdal_calc.py system call?

system('"C:\\Miniconda2\\python.exe" "C:\\Miniconda2\\Lib\\site-packages\\osgeo\\scripts\\gdal_calc.py" --co="COMPRESS=LZW" -A D:/xx.tif --outfile=D:/yy.tif --calc="A*logical_or(A<=0)" --NoDataValue=0')

Here is the error message:

0 .. evaluation of calculation A*logical_or(A>=1) failed
Traceback (most recent call last):
File "C:\Miniconda2\Lib\site-packages\osgeo\scripts\gdal_calc.py", line 326,
in <module>
main()
File "C:\Miniconda2\Lib\site-packages\osgeo\scripts\gdal_calc.py", line 323,
in main
doit(opts, args)
File "C:\Miniconda2\Lib\site-packages\osgeo\scripts\gdal_calc.py", line 272,
in doit
myResult = eval(opts.calc)
File "<string>", line 1, in <module>
NameError: name 'logical_or' is not defined


texts:
This question is very similar to How to define layer's CRS and avoid the CRS dialog in PyQGIS?, but it's doesn't work for me.

This is my try:

vlayer.setCrs(QgsCoordinateReferenceSystem(25830))
iface.mapCanvas().mapRenderer().setDestinationCrs(QgsCoordinateReferenceSystem(25830))

But the dialog pop up  asking for CRS to choose, and when I ask for the CRS of the layer, it show a diferent CRS:

print "CRS: " + vlayer.crs().geographicCRSAuthId()

It show CRS: EPSG:4258 and It must be 25830


texts:
This is a follow up question to my other question: How to speed up raster/map algebra in PostGIS?

I am much attracted to the way how PostGIS handle vectors and rasters in an interoperable manner. A practical issue though (as shown in the link above), is that pqraster is often much slower in many operations when compared to other mainstream packages such as ArcGIS.

My question is specific about the algorithms and methods implementing pgraster (e.g. map algebra):

Are there specific, technical reasons that make pgraster necessarily much slower than, e.g. ArcGIS, or is the slowness only because of a temporary issue or lack of interests/optimization?

For example, is there something inherently more costly (time consuming) in storing regular grid data in a relational database (as with PostGIS) than using a file system (as with ArcGIS)?

I remember reading somewhere that someone joked about pgraster was going to remain slow for the next 10 years, which is why I am curious.

Can someone with the expertise shed some light into this?


texts: Calculate Field Index Based on Feature Attribute Value in PyQGIS?

texts:
Is there a way of finding to which field/column a feature attribute belongs to?

For example, lets say i know feature[x] has an attribute of 0.17, I also know this must be between fields[12:17].

Can I ask QGIS to return the field name with that attribute for that particular feature?


texts: Updating z value of geometry using updateCursor?

texts:
I am trying to update the Z value of a polygon using update Cursor and it just doesn't work!

I have been exploring this site and trying different codes without success.

I have a shapefile with polygons that represents buildings, and some of the vertices have Z values, others have z=0.

What am I doing wrong?

import arcpy

sheet_name = "RPN_2609"
GDB_path = r"D:\Buildings_Height\Case_Study\Part_B\{}.gdb".format(sheet_name)
bldg = GDB_path + "//" + "RPN_2609_BLDG_No_Zero"

with arcpy.da.UpdateCursor(bldg, ["SHAPE@","OID@"]) as cursor:
  for row in cursor:
      for vertice in row[0]:
          for pnt in vertice:
            try:
                if pnt.Z == 0:
                    pnt.Z = 1
                    cursor.updateRow(row)
            except:
                pass


texts: Converting DMS to DD and projecting to NAD83 StatePlane PA South, point is in Ohio

texts:
I'm having issues plotting a point that I converted from DMS to DD on a map projected in State Plane PA South

I'm given a point...

  77° 30’ 10.3” W, 40° 8’ 12.0” N

And need to convert it to decimal degrees for ArcGIS so...

  77.502861,40.136667

Yet when I add this to my map and add a basemap, the point is thrown off to the west in the wrong direction...

My point has the correct projection system...

  NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702

Looking in my Layers -> Properties -> Coordinate System I have the same projection...

  NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702

Yet the point is somewhere off in Ohio not Pennsylvania.

Could someone help me in where I'm going wrong in my process?

I don't know which part is incorrect my conversion or how I'm projecting the data on the map.