texts: Mapproxy generating plain tiles using mapnik

texts:
I am using MapProxy 1.12 and Mapnik 3 for map creation.
Plain tiles are genererating, no map. Some tiles are showing 404 error(The requested tile is outside the bounding box of the tile map.).
sample.yaml
services:
  tms:
    use_grid_names: true
    origin: 'sw'
layers:
  - name: district
    title: District Layer
    sources: [dist_cache]

caches:
  dist_cache:
    grids: [distgrid]
    sources: [dist_tms]

sources:
  dist_tms:
    type: mapnik
    mapfile: /var/www/Mapfiles/district.xml
    layers: district
    transparent: true

grids:
    distgrid:
        base: GLOBAL_GEODETIC
        srs: 'EPSG:32643'
        bbox: [291627,904686,958569,1426831]
        origin: 'sw'
        bbox_srs: 'EPSG:32643'

globals:

district.xml
<?xml version="1.0"?>
<Map background-color="blue" srs="+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs" buffer-size="128" minimum-version="0.7.2">

<Style name="District_style_line">
<Rule>
     <MaxScaleDenominator>6600000000</MaxScaleDenominator>
     <PolygonSymbolizer fill="#ebb1af"/>
     <LineSymbolizer stroke="rgb(0,0,0)" stroke-width="0.98" stroke-linejoin="round"/>
</Rule>
</Style>
<Style name="District_style_text">
    <Rule>
    <MaxScaleDenominator>6600000</MaxScaleDenominator>
    <TextSymbolizer face-name="DejaVu Sans Bold" size="8" fill="black" wrap-width="10" avoid-edges="false" placement="interior">[name]</TextSymbolizer>
    <MinScaleDenominator>2000000</MinScaleDenominator>
    </Rule>
</Style>
<Layer name="District_Boundary" srs="+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs">
    <StyleName>District_style_line</StyleName>
    <StyleName>District_style_text</StyleName>
     <Datasource>
        <Parameter name="dbname">dummydb</Parameter>
        <Parameter name="extent">291627,904686,958569,1426831</Parameter>
        <Parameter name="geometry_field">the_geom</Parameter>
        <Parameter name="host">localhost</Parameter>
        <Parameter name="password">password</Parameter>
        <Parameter name="port">5432</Parameter>
        <Parameter name="srid">32643</Parameter>
        <Parameter name="table">select name,the_geom from districts</Parameter>
        <Parameter name="type">postgis</Parameter>
        <Parameter name="user">testuser</Parameter>
    </Datasource>
</Layer>

map.js (openlayer2)
var layer = new OpenLayers.Layer.TMS('TMS district', 'http://localhost:8080/tms/', {layername: 'district/distgrid', type: 'png', isBaseLayer:true,transitionEffect:"resize" });
map.addLayer(layer);

Tiles generated like


texts: Converting PNG map to polygons with colors and value using ArcGIS for Desktop?

texts:
I need to use this map (PNG file) in which each color is a population density:

I've already located it using Georeferencing. Now I need to convert it to shapefile with polygons, and narrow down its colormap to 5, using ArcMap 10.3.

How can I convert it to polygon, with the right tolerance (because conversion from raster tool just made too may polygons), and how to introduce line as boundaries, not another polygon.  And after that, how to classify about 200 colors in the image into 5, and add value to them?


texts: Unable to plot values for individual fields in my table using ArcGIS for Desktop?

texts:
I'm having problems adding tabular data to an ArcGIS project and need some help! Example of my table is below (please note that the fields have been edited since the capture was taken to remove any spaces or special characters and to shorten them).

ArcGIS displays the x and y data without a problem. However, I am unable to plot the values for the individual fields in my table (example below – I added an ObjectID field at one point after receiving error messages for the lack of one, but it still didn’t read this); when I attempt to do this in symbology it displays no data on the map whatsoever.

Via symbology I was hoping to add a point for each of my fields for the x and y co-ordinates in order to indicate the different values in the three fields, but it only plots one point for the co-ordinates and there doesn't seem to be an option to modify this.

I've had more joy plotting graphs based on this data, but I'd like to add data labels to the graphs - is it possible to do this? Also, on the layer menu on the left hand side and when I create a legend, there is an image of a bar chart with '400' next to it - what does this relate to?


texts: Calculate field using data from excel table

texts:
I need to interpolate data that I have in a excel file using data from an attribute table of a point feature.
I want to create a new field on the attribute table with the result of that interpolation. The problem is that my data is on excel and I don't want to join the excel table with the attribute table. I just want to use the data from excel table to obtain a result for my attribute table.

Is there any way to do that?


texts: OSM Nominatim Server Installation Projection Error

texts:
Environment:

Ubuntu 12.4
Postgresql 9.1.4
PostGIS 2.0.1
Latest Nominatim installation instructions posted at http://wiki.openstreetmap.org/wiki/Nominatim/Installation

The installation is failing with (command built by the Nominatim php build script):

/osm2pgsql/osm2pgsql -lsc -O gazetteer --hstore -C 1100 -d nominatim /home/user/osm_data/central-america.osm

The -l tag is causing the following error:

Projection code failed to initialise

osm2pgsql SVN version 0.80.0 (64bit id space)

ERROR: Error executing external command: /home/user/Nominatim/osm2pgsql/osm2pgsql -lsc -O gazetteer --hstore -C 1059 -d nominatim /home/user/osm_data/central-america.osm.pbf

In trying to troubleshoot, where is the -l tag pulling the projection code from?  I can get the import to run without the -l but the import later fails with an SRID error which I can only suppose results from this initial error.

Can anyone enlighten me on what I should be looking at to fix the issue?

Thanks in advance.


texts: Selecting the smallest number for each ID-group using ArcMap

texts:
I have attribute table with duplicate "ID" numbers, i.e. "1,1,2,2,3,3,....." and there is a unique number value in another field. I want to select the smallest number for each "ID". How can I perform that?
Below is an image of my Attribute Table


texts: Which methods are most suited to classifying urban areas from Landsat imagery?

texts:
I am trying to use Landsat data to map urban extents over time. I need to perform a classification into built up / non built up areas. I am new to image classification. I have access to the following software:

ArcGis
ENVI
Ecognition
GRASS / any other open source software you can recommend.

Which methods and applications are most suited to this analysis? I would appreciate if you could help me understand which methods or applications yield the highest classification accuracy.


texts:
In a recent article by Stewart and colleagues in IJHG I stumbled upon interesting technique of visualizing data using ringmaps. Some more information referenced in the article about this technique here and here.

[Source]

From the article:

  Ring
  maps were created in Adobe Illustrator through the application of a custom
  script that dynamically drew, distributed, and symbolized all graphic map
  elements.  The values for symbolization were read from a Comma Separated
  Value (CSV) file that contained all county attribute data.
  Three county-level ring maps were developed.

Is there a way to automate implementation of a map like that in ArcGIS Desktop 10?


texts: How is the first vertex in a polygon chosen?

texts:
I was wondering about the method/specification of choosing the first vertex in a polygon. Let's say that I have polygon that has 4 vertices, ([0 0] [10 0] [10 10] [0 10]). Which of them will be the first node when iterating through them (and why)?

My thought is that the vertex with minimal x-value is chosen (and if there are several vertices that share the same x-coord then the vertex with minimal y-value among them is chosen),also the nodes above would be represented this way in a GeoJSON:

[[[0 0]
  [10 0]
  [10 10]
  [0 10]
  [0 0]]]

I have not found any source that confirms or denies this.

I am using OSM data that is stored in PostGIS and the polygons will be examined by simple software that I wrote.


texts: GeoServer CQL Filter IN and strToLowerCase

texts:
I am trying to query GeoServer with a CQL Filter using the IN clause.

If I send a request like this I get data back.

http://server.com/geoserver/server/wfs?service=wfs&version=1.0.0&request=getfeature&typename=server:states&CQL_FILTER=state_name%20in%20(%27Texas%27)&outputFormat=application/json

but if I send a url with strToLowerCase around state_name it will not work.

http://server.com/geoserver/server/wfs?service=wfs&version=1.0.0&request=getfeature&typename=server:states&CQL_FILTER=strToLowerCase(state_name)%20in%20(%27Texas%27)&outputFormat=application/json

Is there anyway around this?


texts: Spatial join results different in ModelBuilder and right-click join?

texts:
I am trying to incorporate a simple spatial join in a model but the spatial join options in ModelBuilder are different from the right click-->join options. How can I do the join as in the image below in the data management-->

Above process works fine but I have to do it manually.

The process below fails/doesn't give the result I want and I need to incorporate this in the model as otherwise I have to do this manually and then complete the process in a second model.


texts: How to use PrintTask with preserveScale and WebTiledLayer and when to use the Print dijit?

texts:
I have some enquiries relatively to the use of PrintTask, so as to export the map in various formats:

Should I use the esri/dijit/Print dijit like in this example ?  or it is better to usee the PrintTask like in this ? which one is suggested ?
I want my map to have specific dimensions, and show the entire word at a specific scale? How can I tell the map not to show the current map extent? but the entire world ? I tried setting  like in the examples here the preserveScale to true and the exportOptions to template.exportOptions = {
width: 500,
height: 400,
dpi: 96
};  But the printed map was empty (white) with just the logo on it .
In case I need to use tiled layer as basemap should I only  use the ArcGISTiledMapServiceLayer? I tried with WebTiledLayer but the background map is white, despite the fact that the feature services items are exported/printed


texts: QGIS 3.10 export layer group into arbitrary GIS formats e.g. KML, SHP and etc

texts:
I would like to export the layer group created in QGIS into the GIS file extensions as: .kml, .shp, etc.

So far, the "Export" option available for this layer group gives a chance to export as a "Layer Definition File" in .qlr extension only.

I would like to have more options for saving this groupped layer.


texts: How to get overlapping polylines displayed correctly in Leaflet?

texts:
The situation:
I work on a web application where routes (represented by polylines) are displayed, that are partly overlapping, meaning some segments belong to more routes.

The data:
The polylines are included as GeoJSON files (Feature Collection) within Leaflet. Unless the routes are not overlapping, the display works fine and the routes are highlighted in a different color when they are selected by the user.

The problem:
The problem appears when routes have partly the same segments - where only the route that is later loaded with the GeoJSON file displays correctly, as it is put on the other route that has been loaded before. Thus, some routes (the ones that are beneath the later loaded ones) are not displayed correctly, but cut off at the position where the route on-top overlaps.

Possible solutions:
I have been looking for solutions for several days now, but I cannot figure how to deal with this problem. My question is similar to this question here (https://stackoverflow.com/questions/42251451/leaflet-tooltips-for-overlapping-polylines), but I could not get it to work.

Another way I thought it could work is, when the user clicks on a route, the ID is read out, and this route of the GeoJSON file could be loaded again and put on top of the map. What I don't know is how to access a specific polyline within a GeoJSON file and how to deal with it, when a segment is clicked, where multiple routes overlap.

This is the code I have so far:

var routes = L.geoJson(routes,
{
    style: style_route,
    onEachFeature: onEachFeature
});
routes.addTo(map);

var defaultRoutes = "";
function onEachFeature(feature, layer)
{
    layer.on('click', function (e)
    {
        routeID = feature.properties.RouteID;

        //style the selected route differently
        var clickedFeature = e.target;
        clickedFeature.setStyle(style_routeClicked);

        //reset the style for the route if a new route is selected
        if (defaultRoutes) defaultRoutes.setStyle(style_route);
        defaultRoutes = clickedFeature;
    });
}


texts:
I have the following style for an OpenLayers vector layer.

Color and Name come from Geoserver/PostGIS.

I red that setting labelSelect to true, user can also click the label, not only the point and fire an event such as featureclick. But no reaction. Clicking the point works, clicking the label does not. For my case, the following code is for points, lines and polys. What am I missing?

Thanks

var style_d = new OpenLayers.Style({
'fillColor': '${p_color}',
'fillOpacity': 0.8,
'fontColor': 'black',
'pointRadius': 12,
'fontWeight':'bold',
'strokeColor': 'black',
'strokeWidth': 1,
'fontFamily': 'Arial',
'fontSize': 15,
'label': '${p_name}',
'labelAlign': 'tc',
'labelOutlineColor':'white',
'labelOutlineWidth':3,
'labelOutlineOpacity':0,
'labelSelect': true,
'labelXOffset': 4,
'labelYOffset':4});


texts: Batch scripting r.in.lidar in pyQGIS for LiDAR DEM output?

texts:
I need to use pyQGIS to write a script that will run the r.in.lidar GRASS 7 on a set of LAS files.

I am terrible at Python, but I need this for my task. I would like to process as a batch. How can I set up my script in the PyQGIS console? And how do I set input in the arguments?


texts:
I currently have a layer file - in this case one of the United States.  This layer contains a shape for each state so if I were to open the TOC Control and expand the "United States"  I would see Alabama, Alaska, etc etc.

Currently, if I wanted to I could change the color of all the states through the following code:

IRgbColor rgb = new RgbColorClass();
ColorDialog myColorDlg = new ColorDialog(); ///This is the windows color picker

DialogResult res = myColorDlg.ShowDialog();

if(res == DialogResult.OK)
{
   IFeatureLayer featureLayer = layer as IFeatureLayer;
   IFeatureClass featureClass = featureLayer.FeatureClass;

  rgb.Red = myColorDlg.Red;
  rgb.Blue = myColorDlg.Blue;
  rgb.Green = myColorDlg.Green;

  IGeoFeatureLayer geoFeatureLayer = layer as IGeoFeatureLayer;

  ILineSymbol lineSymbol = new SimpleLineSymbolClass();
  lineSymbol.Width = 1.0;
  lineSymbol.Color = (IColor)rgb;

  ISimpleFillSymbol fillSymbol = new SimpleFillSymbolClass();
  fillSymbol.Outline = lineSymbol;
  fillSymbol.Color = (IColor)rgb;

  ISimpleRenderer simpleRenderer = new SimpleRendererClass();
  simpleRenderer.Symbol = (ISymbol)fillSymbol;
  geoFeatureRenderer.Renderer =  (IFeatureRenderer)simpleRenderer;

  myAxMapControl.Refresh();
}

The question now is, how can I change each state separately.  I'm aware that if I use the built in TOC control I should be able to do this with ILegendClass; however, I wish to use a different TOC Control.


texts:
I am following a tutorial developing a Python plugin in QGIS. After completing the script for the plugin I got the following message the next time I started QGIS. Warning I am a complete novice with Python programming if this is the issue

Couldn't load PyQGIS.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
  File "C:/OSGEO4~1/apps/qgis-ltr/./python\qgis__init__.py", line 26, in
    from qgis.PyQt import QtCore
ImportError: No module named PyQt


texts: What can I do in qgis to manage a WMS frequently unreachable?

texts:
In my project I have a layer linked to an external WMS server (provided by a local institution). The problem is that it is somewhat unreliable: it is down rather often.

Is there something I can do to correctly manage this situation?

At the moment it simply removes the layer, so if I were to save the project I should re-add the WMS later (when it is back on, unpredictable).

A far better behaviour would be just ignoring it and show some kind of warning when trying to get data out of it


texts: The minimum bounding circle of geometry that crosses the 180th meridian

texts:
Background: I've been given a task where I have to calculate the minimum bounding circle of Alaska (for the purpose of running a Roeck test).  I'm using the Census' state level geometries available here.

The problem: Alaska's minimum bounding circle gets calculated "incorrectly". There is a string of islands located between ~172 and 180 degrees longitude, while the rest of Alaska is from -180 to -129 degrees longitude.  This results in my minimum bounding circle covering the whole globe.  While it is technically a correct minimum bounding circle of the given points on a euclidean plane, it clearly doesn't represent the true minimum bounding circle of Alaska.  (I'm using PostGIS's ST_MinimumBoundingCircle.)

The solution: The way I made it work was to

convert Alaska to well known text
subtract 360 from all of the longitudes between 170 and 180 degrees
Put Alaska back into the database.

(The code will be in an answer.)

This is a dissatisfying solution because it's not generalized (it won't work in all cases), and it seems harder than necessary.

My question: What are the options for dealing with this problem?  Has someone written a more intelligent version of the MBC tool that will find the "true" MBC instead of the over-large one?  Is there a simpler way to transpose the geometries than transposing individual points?


texts: Using ArcPy outside ArcMap gives error no module named arcpy?

texts:
I want to use some arcpy function outside arc map
i used following code but it ends up with error no module named arcpy
can anyone help

import sys
sys.path.append("D:\Installed Softares\ArcGIS\Desktop10.1\arcpy");
import sys.arcpy
Radiancetemp = "C:\\Model_Uz\\radianMap"
# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("Float(\"%LE71970282002228EDC00_B6_clip.TIF %\" * 0.03705882) + 3.2", Radiancetemp)


texts: org.geotools.jdbc.JDBCDataStore stores shape files in the database but files stored in own tables

texts: Using geotools to store shape files in a PostgreSQL database. Each shape file is being stored in it's own table in the database. Is there a way to store the files in a single table? Also, the tables created for each file have a different structure. Is there a way to check the shape files to know the difference in tables that will be created and group the files accordingly?


texts: ST_Within for TIN or 3D objct in PostGIS

texts:
I would like to query if a point is located in the TIN object. My TIN object is closed (convex hull) object (i.e. surface of a round object). I tried to use ST_Within, but the error message returns as:

  ERROR: Unknown geometry type: 15 - Tin

My query is as follows:

SELECT ST_Within(mypoint,mytinobject) As mypointinmytinobject
FROM
(SELECT ST_GeomFromText('POINT(x,y,z)', SRID) As mypoint,
(SELECT geom FROM mytintable WHERE gid=1) As mytinobject
) As myobjects;

What are the equivalent 3D query functions like ST_Within for 3D TIN object in PostGIS?

It is a complex 3D anatomical region. It is not solid, it is a surface. If it is represented as Polyhedralsurface, how could we query if a point is located inside this closed surface (or a region). Do I need to make it solid body. Could it be done in PostGIS?


texts: Sort layers in table of contents in QGIS 3?

texts:
When I load a bunch of layers into QGIS, they seem to load in reverse alphabetical order. I'd like to sort them alphabetically.

In this thread a solution is proposed, but that does not work in QGIS 2.18 or 3.4:

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

mw = iface.mainWindow()
lgd = mw.findChild(QTreeWidget, "theMapLegend")  # get ref to object by type/objectName
lgd.sortItems(0, Qt.AscendingOrder)  # sort first column (Qt.DescendingOrder to reverse)

How can I do this in QGIS 3?
I tried:

from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui import *
from qgis.utils import iface

mw = iface.mainWindow()
lgd = mw.findChild(QTreeView, 'theLayerTreeView')
lgd.sortItems(0, Qt.AscendingOrder)

It "works" until the last line, then I get the error:

  Traceback (most recent call last):
    File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 90, in runcode
      exec(code, self.locals)
    File "", line 1, in
  AttributeError: 'QgsLayerTreeView' object has no attribute 'sortItems'

I don't know where to search for a method that will sort the layers as in the first code example.


texts: Individual Tree Segmentation segmenting each tree into many pieces

texts:
I am a first time user trying to help my community inventory trees over a threshold height using public Lidar data.  My goal is to extract hulls and/or approximate GPS coordinates of each tree exceeding the target height.
I followed along with the individual tree segmentation examples (1,2) but in my results trees are generally segmented into 4-8 parts based on the colorization in the image (included).
I first thought this had to do with the watershed parameters, but I didn't hit on the right combo in the workflow.  Later I tried the li2012 algorithm and my results were no better.
Can someone suggest either the key parameters I should be changing, or a different workflow I should be following?
print(Sys.time())
library(lidR)

las = readLAS("/Users/brianvanvoorst/Desktop/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.las")
las = classify_ground(las, csf())
print ("Classify done")
las = normalize_height(las, tin())
print ("Normalize done")
algo = pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
print ("Pitfree done")
chm  = grid_canopy(las, 0.5, algo)

plot(chm, col = height.colors(50))
# smoothing post-process (e.g. two pass, 3x3 median convolution)
ker = matrix(1,3,3)
chm = focal(chm, w = ker, fun = median)
chm = focal(chm, w = ker, fun = median)

plot(chm, col = height.colors(50)) # check the image
algo = watershed(chm, th = 4)
las  = segment_trees(las, algo)

# remove points that are not assigned to a tree
trees = filter_poi(las, !is.na(treeID))

plot(trees, color = "treeID", colorPalette = pastel.colors(100))
print(Sys.time())

Okay here is my latest source and result
require(lidR)
require(rlas)
require(rgdal)
require(tictoc)

las <- readLAS("/Users/brianvanvoorst/Desktop/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.las", filter="-keep_class 1L")

#dtm <- grid_terrain(las, algorithm = knnidw(k = 8, p = 2))
# Error: No ground points found. Impossible to compute a DTM.
#las_normalized <- normalize_height(las, dtm)

# Create a filter to remove points above 95th percentile of height
filter_noise = function(las, sensitivity)
{
  p95 <- grid_metrics(las, ~quantile(Z, probs = 0.95), 10)
  las <- merge_spatial(las, p95, "p95")
  las <- filter_poi(las, Z < p95*sensitivity)
  las$p95 <- NULL
  return(las)
}

las_denoised <- filter_noise(las, sensitivity = 1.2)

chm <- grid_canopy(las_denoised, 0.5, pitfree(c(0,2,5,10,15), c(3,1.5), subcircle = 0.2))

plot_dtm3d(chm)

ker <- matrix(1,5,5)
chm_s <- focal(chm, w = ker, fun = median)

algo <- watershed(chm_s, th = 4)
las_watershed  <- segment_trees(las_denoised, algo)

# remove points that are not assigned to a tree
trees <- filter_poi(las_watershed, !is.na(treeID))

# View the results
plot(trees, color = "treeID", colorPalette = pastel.colors(100))


texts: Position a polygon to maximise number of points within it

texts:
I have a data set with over 1 million points. I need to find location for polygons (circle) that will maximise the number of points within it.

Requirements:

Radius of the circle depends on where the centre of the circle lies (different areas will have different radii)
Some points may have over 5000 neighbours within the circle, so solution will need to manage the size of the data (and connections found)

Assumptions:

If this makes analysis easier, points in one circle may not lie in another circle (i.e. circles don't overlap)
There are no limits on how many circles can be made

I originally tried to centre the circles on the points themselves so I can run it through a loop, but ideally I'd like to maximise the points wherever it lies.

My previous attempt included looking at number of points within radius (relationships) from each point so that I can order them DESC of density, and for the highest density remove its child points from other relationships. Then using the next biggest density and so on.. This seems very inefficient (size of the interim data became too big and slow).

I have a feeling using PostGIS or raster type solutions will be best suited for this? Can anyone please throw some light?


texts: Can QGIS preserve layers when exporting a PDF?

texts: I saw this posted about a year ago (Can QGIS preserve layers when exporting to PDF?) but figured it was worth asking again -- is this possible? My usual workflow with ArcGIS is to export to AI then edit the layers that way, is there some way to do that from QGIS? Now I just get everything in one huge layer. Thanks.


texts:
Is it possible to get a list of street names (without house numbers) by bounding box with Overpass or Nominatim?
I got this query, but it is with house numbers so the responce JSON is big

[out:json][timeout:25];
// gather results
(

  node["addr:street"]({{bbox}});
  way["addr:street"]({{bbox}});
  relation["addr:street"]({{bbox}});
);
// print results
out body;
>;
out skel qt;


texts: Automatically move points so that they don't overlap (reposition into cluster pattern)

texts:
I have workorder points in a SDE.ST_GEOMETRY feature class (Oracle 18c/10.7.1. geodatabase).
I serve up the point feature class to an external system via a Map Service.

Problem:
The workorder points in the feature class often overlap. So, in the external system, users are having issues where they don't realize there are multiple points at a given position.
To solve this issue, I want to create a trigger in the feature class that would automatically move features as they are created. The trigger would reposition a point so that it does not overlap with other points (possibly into a cluster pattern).
What would it take to do this?


texts: CRS USA_Contiguous_Albers_Equal not working

texts: I have selected the USA_Contiguous_Albers_Equal_Area_Conic, yet the map is showing a different projection.  It worked when I started the project and something went wrong after I did a "print composer".  Now it doesn't work even if I start a new project and re-do all the steps.  Any idea how I can fix this?


texts: Why d3.js works only with GeoJSON violating right-hand rule?

texts:
Why d3.js works only with GeoJSON that is left-hand formatted?
{
    "type": "FeatureCollection",
    "features": [
      {
        "type": "Feature",
        "properties": {},
        "geometry": {
          "type": "Polygon",
          "coordinates": [
              [
                [0, 0], [10, 0], [0, 10], [0, 0]
              ]
            ]
          }
        }
      ]
    };

is supposedly well formed and not broken - https://geojsonlint.com/ shows it as follows
while d3.js shows it as a hole in the world (filled area across entire world except that area)

My, likely bungled, d3 attempt code is hosted at https://matkoniecz.github.io/wat/valid.html (with correct GeoJSON). Version with points in reverse order, raising warnings elsewhere works: https://matkoniecz.github.io/wat/map.html
What is going on here? I am likely missing something but I run out of ideas. Maybe d3 is simply not supporting GeoJSON? And all examples that I found are broken/confusing because people who made them were unaware that GeoJSON is not supported?
Code is also posted below:
<!DOCTYPE html>
<html>
<head>
    <meta charset="utf-8" />
</head>
<body>
    <div id="content">
      <svg width="800px" height="600px">
        <g class="generated_map" id="generated_map"></g>
      </svg>
    </div>

    <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/6.6.2/d3.min.js"></script>
    <script>
    const width = 800;
    const height = 600;
    var geojson = {
    "type": "FeatureCollection",
    "features": [
      {
        "type": "Feature",
        "properties": {},
        "geometry": {
          "type": "Polygon",
          "coordinates": [
              [
                //[-179, 80], [179, 0], [179, 80], [-179, 80] // works? why?
                //[0, 0], [10, 0], [0, 10], [0, 0] // hole in the world
                [0, 0], [0, 10], [10, 0], [0, 0]
              ]
            ]
          }
        }
      ]
    };

    var projection = d3.geoMercator().fitSize([width, height], geojson)

    var geoGenerator = d3.geoPath()
    .projection(projection);

    function update(geojson) {
      var u = d3.select('#content g.generated_map')
        .selectAll('path')
        .data(geojson.features);

      u.enter()
        .append('path')
        .attr('d', geoGenerator)
        .attr("stroke", 'black')
        .attr("stroke-width", "1")
        .attr("fill", "blue")
    }
    update(geojson);
    </script>
</body>
</html>


texts: Intersection algorithm that correctly handles 180° meridian and poles?

texts:
Is there a well known intersection algorithm that correctly handles 180 meridian and poles?

For example, let's say we have a list of latitude and longitude values that represent Antarctica.  Let's say we also have a simple polygon that represents an airplane.  We want to know if the airplane is over Antarctica.

Generic 2D polygon intersection does not work in this case if you simply use latitude for y and longitude for x, because the flat coordinate system will have edges at -180 and 180 longitude, and -90 and 90 latitude.  Antarctica will go off the page on three of the edges.


texts: How can I merge segments to one segment in QGIS (picture inside) and than get the center of each road?

texts:
I want to run a process that will make a table with lots of roads into distinct table of roads and than calculate the center of each road.

How can I auto merge segments with same name to one long segment ?

Attached example of 1 road built from 3 different segments , I want to make it one road and have it center point. (so I can zoom on it from OpenLayers)


texts: Iterrtuples + GeoPandas DataFrame + Plotting

texts:
I am not well versed in Python (or GeoPandas) and I am currently having a simple issue with plotting. To spot some inconsistencies, I need to individually plot all the geometries in the geodataframe. If you are wondering why I am using enumerate() and itertuples() is because I need to find where the inconsistencies are (but I haven't implemented yet in the code!).
Here's the code. Hopefully is something simple enough:
for i, row in enumerate(faci_shp.itertuples(index=False)):
  row.plot()

But I get this error:

--------------------------------------------------------------------------- AttributeError                            Traceback (most recent call
last)  in ()
1 choma_shp = faci_shp.iloc[facilities, :].copy()
2 for i, row in enumerate(faci_shp.itertuples(index=False)):
----> 3   row.plot()
AttributeError: 'Pandas' object has no attribute 'plot'

Is it clear?
Thanks


texts: Can you go back and view the equation you put into the field calculator?

texts: If you have an attribute table column that you populated with values using the field calculator, is there a way to go back and view the calculation you input? This is days later outside so I am no longer in an edit session but I would like to check it again since there are some weird values.


texts: Saving all the project layers in one folder in QGIS

texts:
There is a nice way to upload all the layers from one folder in QGIS.
Is there  a way to download all the layers from a project and save them in a new folder?


texts: Adding field to ArcGIS Map Service without feature access?

texts:
I have published a mapservice without feature access on a remote machine only with a file geodatabase (shapefile on the disk).  Because of limits of authority to create a spatial database(don't have permission to create sqlserver express 2008 for build arcsde to enable feature access for this map service).

However, I have to add a new field into the attribute table of this mapservice via the ArcServer Javascript API .Just like the operation on ArcGIS desktop ,as shown below:

Is this possible? Any hints to implement it?


texts: Is a PostGIS raster column more efficient than a geometry column in this scenario?

texts:
We have several different PostGIS tables that have a timestamp column and a geometry column. The geometry column contains polygons that are queried by a GeoServer instance to be served for display.
The same information is stored as vector points for each coordinate in "analysis" tables. The only other use case we have is querying these analysis tables for single point values, at specific timestamps.
The data is loaded with GDAL from raster files, and is always global. My question is, would there be any benefit to switch the pipeline to store the rasters in raster columns (either in-db or out-db) instead?
I've seen advice against using raster columns if not performing relevant computations, but the official PostGIS documentation doesn't offer any information on the drawbacks of using them.


texts: Viewing what properties exist in GeoJSON file loaded into MapShaper?

texts:
Suppose I have a large GeoJSON file loaded into Mapshaper. I want to convert a property to use as an ID field on export, but I don't know what the property naming schema is. The file is too large to conveniently check in a text editor.

Is there a way within Mapshaper to view available properties? I've tried using -inspect, but you have to select one record by a property -  which is hard when you don't know what the identifying property's name is! I can't find anything relevant searcing for properties in the command reference.


texts: Symbology levels are ignored in ArcGIS 10

texts:
I have a problem migrating from a ArgGIS 9.3 to a 10 geodatabase (both sde and  file gdb). The symbology levels stored in the 9.3 database are not transferred.

When migrated to ArcGIS 10.0, symbology levels are ordered by alphanumerical order, and not the order assigned by me.

I tried to migrate the data both by copying the featureclass and by upgrading the database -- both with the same dissapointing result: The symbology levels go back to the alphanumerical order of the categories. This is a BIG problem for me, since I have a massive amount of symbology levels in lots of featureclasses each having different representations.

Have any of you similar experiences? More importantly: Is there a fix for this?


texts: Display Labels in Google Maps from Spreadsheet Mapper

texts:
How would I go about displaying labels on placemarks on the map canvas where the source kml comes from the Spreadsheet Mapper?

This particular use case is where we want to share and control the source spreadsheet in Google Docs, and have a display interface that is fairly ubiquitous, i.e. on thin blade client browsers, android browsers. Not either of these display devices will allow the use of Google Earth (thin blades don't have it installed cause of those well meaning network admins.... and droid GE is not able to display limited share files AFAIK).

It is highly desired to be able to edit the source data in a GDocs SS, since the placemark locations will rarely change, but the labels and comments should pretty straightforward edits.


texts: Exporting polygon coordinates from R to Google Earth

texts:
I've created a list of lat long coordinates that represent the perimeter of a location. I wish to plot this polygon as a line in Google Earth rather than as push pins.

I have the XML/KML file format for plotting the polygon in Google Earth but is it possible to script this in R to populate this particular field in the KML file?


texts: Simplifying geometries and calculating non-spatial averages

texts:
I'm quite new to using proper spatial tools, and have a simple question that doesn't immediately show up on google.

I have a table that contains multiple rows of empiric data (say, speed of a vehicle) that also has spatial information connected to it. Every speed measure (real number) has a linestring geometry to it.

Many of those Linestrings are touching each other on the end points.

My task is to:

visualize all groups of touching lines as separate multilines
simplify each of the resulting multilines to reduce the amount of points and consequently data to be displayed (i.e. when we "zoom out", so I need configurable loss of precision)
aggregate average speed for every segment in each simplified multiline

Now, I understand that spatial simplification and scalar aggregation probably use different algorithms, but what's the common approach to this?

UPDATE:

Here's some data:

geom                                             |speed
LINESTRING(25.51369 65.01514, 25.51362 65.01521) |21.6324731983004
LINESTRING(25.51362 65.01521, 25.51369 65.01514) |20.7681850053341
LINESTRING(24.95117 60.1898, 24.95112 60.18981)  |15.5648781827192
LINESTRING(24.8263 60.20353, 24.82685 60.2035)   |19.7643783824494
LINESTRING(24.82685 60.2035, 24.8263 60.20353)   |15.2009950253671

I have understood that ST_Union(), ST_LineMerge() and ST_Simplify() can take care of building a joined multiline and reducing its complexity.

Part 3 of the task still stands: how do I calculate average speed values for various parts of the newly created and simplified multiline?


texts: Export from earthengine API using a Gservice account?

texts:
In earthengine Python API, I connect using a Google service account :
service_account = 'test-sepal-ui@sepal-ui.iam.gserviceaccount.com'
credentials = ee.ServiceAccountCredentials(service_account, 'test.json')
ee.Initialize(credentials)

I would like to launch an exporting task on on one of my personnal folders.
I then created a folder in my assets test and authorized the test-sepal-ui@sepal-ui.iam.gserviceaccount.com as a writer.
Then I try to launch the following code :
toto = ee.FeatureCollection("FAO/GAUL/2015/level0").filter(ee.Filter.eq('ADM0_NAME', 'France'))

task_config = {
    'collection': toto,
    'description':'toto',
    'assetId': 'projects/earthengine-legacy/assets/users/bornToBeAlive/sepal_ui_test/toto'
}
task = ee.batch.Export.table.toAsset(**task_config)
task.start()

but when I go to the codeEditor, no task is launched and if I print the task in my notebook I get :

<Task 4TX2KTMGHHWJZE5K2JQ4XCVO EXPORT_FEATURES: toto (UNSUBMITTED)>

Is there a way to perform this trick ?
context
I try to build testing function and I need to import and export assets. I wanted to use the public one but it's impossible so I tried to use a personal one instead.


texts: MapServer Mapfile: One layer, disable the label but not the point for first two zoom levels

texts:
This map is not using OpenLayers, just MapServer with a mapfile.

There is a layer that displays both the point and the label.

Looking through the documentation, I'm not entirely sure what the best way is.

Edit: Ideally (not even sure if this is possible), I would like to disable the label for the first two zoom levels but display the associated point. Is this even possible if they are tied together as one layer?

Edit2: Looks like LABELMAXSCALEDENOM is what I'm looking for, but how to configure correctly? I know it's

 LABELMAXSCALEDENOM integer

but how do you exactly know the scale of your map?

Currently:

MAP
  SIZE 128 768
    EXTENT -180 -90 180 90

I tried LABELMAXSCALEDENOM equal to some different integers. Didn't notice any difference: the labels were just gone completely, no matter what the zoom.

Any input is appreciated.


texts: Selecting features that do NOT intersect in PostGIS

texts:
Using PostGIS, I just want to select points that fall outside of polygons.  Ultimately this is the inverse of the ST_Intersects, as far as I can see it.
Example: I have a taxlot layer and a address point layer.  I assume I should be using the ST_Intersects, but how do I tell it to do the reverse selection?  I thought maybe adding a NOT statement in front of the code below, but that did not work.
CREATE table t_intersect AS
SELECT
  hp.gid,
  hp.st_address,
  hp.city,
  hp.st_num,
  hp.the_geom
FROM
  public.parcel as par,
  public.housepoints as hp
WHERE
  ST_Intersects(hp.the_geom,par.the_geom);


texts:
I'm using QGIS 2.8.

I need to calculate the distance between points depending on points ID (not nearest point). I have a CSV table including points ID and coordinates. Each point is referenced on 2 different lines with a different pair of coordinates. I need to know the distance between the 2 positions of each point.

I have tried fTools plugin but the Distance matrix function returns nearest points measures (I had beforehand splitted the CSV table in 2 sets of unique points ID).

Any idea about a plugin or function that would enable to do this?


texts: Copy Raster for a Raster Mosiac with Function produces corrupted data

texts:
I am trying to export a ArcGIS raster mosaic dataset to a new raster that has a "Elevation Void Fill" function on it to fill holes in the data, but I am finding the function data is getting corrupted when exporting to a new raster.

Any ideas on how to remedy this?

Orignal .tif file

Mosaic dataset raster with function to fill gaps

Exported to new raster by Copy Raster Tool


texts: I have a library of .bmp logos that I would like to use to symbolize points in QGIS. Does QGIS only support .svg logos?


texts: Modifying a Python script called in a shell command in Google Colab

texts:
I am relatively new to coding and am working on a script to average 24 hourly GRIB2 files to daily files and output to a TIF. I am currently doing this in Google Colab by calling a GDAL python script called gdal_calc.py (downloadable from GDAL), which I have uploaded to my drive.
My code looks like this (for November 29):
import os

! python /content/drive/MyDrive/NOAA/grib2_files_test/gdal_calc.py \
    -A  /content/drive/MyDrive/NOAA/grib2_files/202011/2020112900/postprd/nearsfc_smoke_00.grib2 \
    -B  /content/drive/MyDrive/NOAA/grib2_files/202011/2020112900/postprd/nearsfc_smoke_01.grib2 \
    -C  /content/drive/MyDrive/NOAA/grib2_files/202011/2020112900/postprd/nearsfc_smoke_02.grib2 \

etc (D-V)

    -W  /content/drive/MyDrive/NOAA/grib2_files/202011/2020112900/postprd/nearsfc_smoke_22.grib2 \
    -X  /content/drive/MyDrive/NOAA/grib2_files/202011/2020112900/postprd/nearsfc_smoke_23.grib2 \
    --outfile=/content/drive/MyDrive/NOAA/grib2_files/averaged2/20201129.tiff  \
    --calc="(A+B+C+D+E+F+G+H+I+J+K+L+M+N+O+P+Q+R+S+T+U+V+W+X)/24"

I have data for all of 2020, but there are missing days so I haven't tried to make a monthly loop. I would however love to create a "dd" variable that I could change for each day instead of manually changing the day value in each of the file names and output name. I think my attempts to do this have not worked because I am calling an external script. Any ideas on how to create this variable?


texts: Black map after attempt to style by column in QGIS

texts:
I'm very new to QGIS.

I have a shapefile of NYC census tracts and 2 csv info files with tract info joined with it. I'm trying to style the individual tracts to be filled with a different color based on the income of that tract (that info is in one of the columns of the csv file).
In order to do that, I went to Layer Properties > Style and chose Categorized from the top drop-down. I then selected the income column and hit classify to display the color-coded contents. Once everything looked good, I pressed okay, and suddenly the map turned black, momentarily flashing random tracts until stopping.

It looks like this:

Any idea what I'm doing wrong?


texts:
I would like to incorporate MS Office into my GIS program.  I have been using QGIS and was wondering if this would be possible?

I am trying to create a one stop shop for my staff to be able to write reports and print maps from one program instead of going between programs.

An example would be I have a complete form with blank fields.  The fields would be pulled from the map attributes.


texts: Set a rule for grow region in eCognition

texts:
I am using the grow region algorithm in eCognition sofrware in order to merge my candidate class.

I want to set the rule : two candidated objects are merged only if the length/width ratio of the target object increases after the merging.

How can I perform this rule?


texts: Should Endmembers be Normalized for Linear Unmixing?

texts:
In biological spectral microscopy, the standard in the field is to use normalized endmembers in linear unmixing algorithms.  In most cases, the endmembers are normalized to the maximum value.  I have not found a reference that explains why endmembers are normalized prior to linear unmixing.  By normalizing, are we not removing information about relative intensity of objects?

I am asking my question here because I am curious if endmember normalization is routine in remote sensing as well.


texts:
Does anybody know how to plot any attribute of a shapefile with R?

For example, like this.


texts: Projecting Data with Transformations - Alignment Issue

texts:
Problem:

I am having problems getting data to align properly with basemaps layers.  The image below shows the difference between roads.

The centerlines shapefile's coordinate system is GCS_North_American_1983 while the basemap is WGS_1984_Web_Mercator_Auxiliary_Sphere.

Things I have tried:

I have tried projecting the centerlines shapefile into WGS_84 to get it to better alignment.  I have tried using WGS_1984_(ITRF00)_To_NAD_1983 and NAD_1983_To_WGS_1984_5 transformations without much luck. I have tried changing the projection of the data to different ones thinking it may have been wrong at the start.

Shapefile Info:

Comparing the data and the basemaps, it looks like there is about a 40 meter difference.  I am assuming this is a projection issue because it is hard to imagine the data being off by that much when it was created.

My data is in western Pennsylvania.  The shapefile's extent is: top 41.892328dd, right -77.075914dd, bottom 39.625013dd, left -80.565890dd.  The basemap is the photo is from ArcMap, but the same issue is there is another program I am using with Google Maps as basemaps.

These shapefiles contain M and Z values, but I am not sure if that would really factor into the alignment.

We are using this data for offline field collection, so it needs to match (more closely) the basemaps we will show it on.  All of the shapefiles I have are off according to the basemaps, too.


texts: Apply transparency on pixels in an image with PyQGIS

texts:
I'm trying to grab the pixels that intersect with my polygon shapefile and apply transparency to them (that is, I'm trying to create a hole in my raster)

For that, I get all the pixels that are in the footprint of the shapefile in a list (see code below)

Then I used the QgsRasterTransparency.TransparentThreeValuePixel () class to apply transparency, but I get an error message:

  Error message: AttributeError: 'Dataset' object has no attribute
  'renderer'

Code:

from PyQt4.QtCore import *
from osgeo import gdal, ogr
from osgeo.gdalconst import *
import numpy as np

iface = qgis.utils.iface

layers = iface.legendInterface().layers()
raster = gdal.Open("F:\\MASTER GEOM\\MasterGeomatique_data\\PLEIADES\\PLEIADES_20130415_COLNAT.TIF")
print(raster)

transform = raster.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]
print xOrigin, yOrigin

shp = ogr.Open("F:\essai.shp")
lyr = shp.GetLayer()
feat = lyr.GetNextFeature()
geom = feat.GetGeometryRef()

if (geom.GetGeometryName() == 'POLYGON'):
    ring = geom.GetGeometryRef(0)
    numpoints = ring.GetPointCount()
    pointsX = []; pointsY = []
    for p in range(numpoints):
                lon, lat, z = ring.GetPoint(p)
                pointsX.append(lon)
                pointsY.append(lat)
                #print (pointsX)
xmin = min(pointsX)
xmax = max(pointsX)
ymin = min(pointsY)
ymax = max(pointsY)
print xmin, ymin
print xmax, ymax
 # Specify offset and rows and columns to read
xoff = int((xmin - xOrigin)/pixelWidth)
yoff = int((yOrigin - ymax)/pixelWidth)
xcount = int((xmax - xmin)/pixelWidth)+1
ycount = int((ymax - ymin)/pixelWidth)+1
print xoff, yoff
print xcount, ycount

bandList = []
# retrieve pixel on region of interest
for i in range(raster.RasterCount):
    band = raster.GetRasterBand(i+1)
    data = band.ReadAsArray(xoff, yoff, xcount, ycount)
    bandList.append(data)

# Apply transparency using QgsRasterTransparency.TransparentThreeValuePixel()
active_layer = raster
raster_transparency  = active_layer.renderer().rasterTransparency()
ltr = QgsRasterTransparency.TransparentThreeValuePixel()
tr_list = []
ltr.min = np.min(bandList)
ltr.max = np.max(bandList)
ltr.percentTransparent = 100
tr_list.append(ltr)
active_layer.renderer().rasterTransparency().setTransparentThreeValuePixelList(tr_list)
active_layer.triggerRepaint()

# Refresh canvas
QgsMapLayerRegistry.instance().addMapLayer(raster)

# Message d'erreur : AttributeError: 'Dataset' object has no attribute 'renderer'


texts: Reclassify on land cover raster won't work

texts:
I'm new to GIS so I apologize if this question is basic but I couldn't find my answer online.

I'm trying to reclassify a land cover raster from the National Land Cover Database. I want to reclassify into a three category raster: 1 - natural areas, 2 - agricultural areas, and 3 - developed areas.

I am going to Spatial Analyst > Reclassify and then I reclassify by unique values and enter my new values into their respective places. When I reclassify, I only get one value: 2. And it doesn't even correspond to what I've classified as 2.

What can I do to fix this? I have no clue what I'm doing wrong.


texts: Extract the outer points of a Point Shape into a Line Shape?

texts:
I've installed QGIS 2.8.2 since a few days and have no experience with this application. I try to extract the outer points of a Point Shape into a Line Shape, like a shell for the inner points. All points of the shape has to be in or on the border of the line.

Is there any possibility or function in QGIS, to do this ? (see figure 1 below)


texts:
I use OpenLayers to display an OpenStreetMap map and on top of that I want to add a vector layer with features. But it seems that I have a problem with the projection of the vector layer or the features which all use GPS data.

The features are only correctly applied to the map if the map is maximally zoomed. As soon as I zoom out the features move from the original point on the map and if I zoom out the map long enough, they form a circle.

Here is how I set the features:

var iconFeature = new ol.Feature({
geometry: new ol.geom.Point(ol.proj.transform([<somelonvalue>, <somelanvalue>], 'EPSG:4326', 'EPSG:3857')),
popuptext: '<sometext>',
});

iconFeature.setId('$featureID');

And here is how I created the vector layer:

vectorSource = new ol.source.Vector({
features: iconFeatures
});

var vectorLayer = new ol.layer.Vector({
source: vectorSource
});

map.addLayer(vectorLayer);

Here is the outcome if I zoom out the map:

I believe I did the transformation correctly.

Does anybody have a clue as to what I am doing wrong?


texts:
I'm having trouble getting data from a createLayer layer to populate a CARTO map when hosted on GitHub. The map itself draws but the data is missing, although, it works when hosted locally. Also, I am able to view data using a CARTO vis.json URL hosted on GitHub. Below is the code, which is a little messy since I'm fairly new to CARTO.

Has anyone encountered a similar problem in the past and found a solution?

      var tableName = "powerplants_us_201603";

      // public layer in CARTO account for data queries
      var layerSource = {
        user_name: 'mishmashmaps',
        type: 'cartodb',
        sublayers: [{
          sql: "SELECT * FROM " + tableName, // All power plants in dataset
          cartocss: $("#simple").text() // Simple visualization
        }]
      }

      var sublayer;

      var map_object = new L.Map('map', {
        center: [42.6526,-73.7562], // Albany
        zoom: 7
      });

      L.tileLayer('http://{s}.basemaps.cartocdn.com/dark_nolabels/{z}/{x}/{y}.png',attribution: '&copy; OpenStreetMap contributors'
      }).addTo(map_object);

      cartodb.createLayer(map_object,layerSource)
        .addTo(map_object)
        .done(function(layer) {
          sublayer = layer.getSubLayer(0);
          createSelector(sublayer);
        })
        .error(function(err) {
          console.log("error: " + err);
        });
    }


texts: SLD file GeoServer: Avoiding duplicate labels

texts:
I have an SLD file for a road layer in GeoServer. I wish only to show the labels for that layer (no symbol for the line objects) and I wish to only see them once for each road when I view my layer in for instance an OpenLayers map viewer. Right now I have many duplicate road names showing when I zoom in close.

I tried to use functions such as: VendorOption: name="group"true and VendorOption name="repeat">0< etc. But that didn't help me.

To show how the problem looks:

The SLD I got now is this:


texts: I own a Garmin Oregon 600. I know that it has some communication capabilities, but I have not found on user guide how to use it. Does anyone of you know how to connect this device by USB/Serial in order to fetch its stream of data? I am confident with data stream and serial communication. I just have not found how to setup my device in order to turn it as a recognizable COM device and make it stream its data. I use windows and linux indifferently.


texts: How to figure out the SRID of a point by comparing it to another point with known SRID?

texts:
I have a set of store locations defined in the coordinate system of Bottincarto, which I need to display on Google Maps iOS.
Here is an example of coordinates for Notre-Dame-de-Paris:

600926 2428323 (Notre-Dame-de-Paris on http://maps.bottincarto.com/...)
48.852968 2.349902 (Notre-Dame-de-Paris on http://www.google.fr/maps/...)

Here is an example of coordinates for Trafalgar Square:

428285 2726837 (Trafalgar Square on http://maps.bottincarto.com/...)
51.508039 -0.128069 (Trafalgar Square on http://www.google.fr/maps/...)

I already tried a linear regression, but it is not even good as an approximation.
How can I reverse engineer BottinCarto projection system?


texts: Is there any known reason for QGIS recurrently crashing in windows7?

texts:
I use (and been a big fan of) QGIS daily for 7 months now, and only recently - past couple weeks - been experiencing problems with QGIS recurrently crashing. This is even more serius cause it happens when i try to save the editions I'm doing.
I'm using a windows7 machine and i'm performing the same task for 2 months now which consists on visually interpreting a satellite image to digitalize features to a shapefile. I haven't change anything and I don't understand why this occurs but it's very annoying to have QGIS crashing every half an hour...
Any ideas? similar experiences?

Thanks,
Vasco Nunes

PS: I'm using version 1.7.3 but this has been happening since I had 1.7.2


texts: QGIS3 python: right(QTextStream) - getting unexpected type ' str'

texts:
I have a field type double that I need to reduce from 9 digits to 8 ie 123456789 -> 23456789. If I do it as an integer I loose leading zeros ie 10002345 -> 2345. I am attempting to use refactor fields to make the field a string then use right(x,8). Refactor is correctly returning it as a string field, however I continue to get the 'TypeError: right(QTextStream): argument 1 has unexpected type 'str'' error.

processing.run("qgis:refactorfields", {'INPUT':CMLot,'FIELDS_MAPPING':[{'expression': '"CADID"', 'length': 10, 'name': 'CADID', 'precision': 0, 'type': 10}],'OUTPUT':'L:\CMLot_test.tab'})
layer = QgsVectorLayer('L:\CMLot_test.tab', "CMLot_test", "ogr")
if not layer.isValid():
    print ("Layer failed to load!")
QgsProject.instance().addMapLayer(layer)
CM_test = QgsProject.instance().mapLayersByName("CMLot_test")[0]
iface.setActiveLayer( CM_test )
layer = qgis.utils.iface.activeLayer()
features=layer.getFeatures()
layer.startEditing()
for f in features:
    m=(f.attributes()[0])
    b=right(m,8)
    f.setAttribute(f.fieldNameIndex('CADID'),b)
    layer.updateFeature(f)

I thought the right function required a string field so I am at a bit of a loss of how to return the 8 right hand digits from the 9 digit integer without loosing the leading 0s.

Am I refactoring to the wrong field type or is my method for reducing 9 to 8 completely off the mark?


texts: QGIS geotiff extract world fails with AttributeError: 'ExtractThread' object has no attribute 'fileProcessed'

texts:
I am trying extract the projection and create the world file and prj file from  a geotiff but I get a failure.

I have the GDAL plugin enabled and I know it is a geotiff with the projection embedded in the tif file. When I open it with qgis browser, it displays all of the information including the layer original source projection. When geotiff is opened in qgis it shows the coordinates correctly.

I actually can overlay another geotiff and they match up perfectly. But I get an error trying to extract the projection. Any ideas?

Error is as follows:

An error has occurred while executing Python code:

AttributeError: 'ExtractThread' object has no attribute 'fileProcessed'
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS2~1.16/apps/qgis/./python/plugins\GdalTools\tools\doExtractProj.py", line 111, in accept
    self.extractor.fileProcessed.connect(self.updateProgress)
AttributeError: 'ExtractThread' object has no attribute 'fileProcessed'

Python version: 2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]
QGIS version: 2.16.3 Nødebo, 7d00cae

Python Path:
C:/PROGRA~1/QGIS2~1.16/apps/qgis/./python/plugins\processing
C:/PROGRA~1/QGIS2~1.16/apps/qgis/./python
C:/Users/john/.qgis2/python
C:/Users/john/.qgis2/python/plugins
C:/PROGRA~1/QGIS2~1.16/apps/qgis/./python/plugins
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win32.egg
C:\PROGRA~1\QGIS2~1.16\bin\python27.zip
C:\PROGRA~1\QGIS2~1.16\apps\Python27\DLLs
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\plat-win
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGIS2~1.16\bin
C:\PROGRA~1\QGIS2~1.16\apps\Python27
c:\progra~1\qgis2~1.16\apps\python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
c:\progra~1\qgis2~1.16\apps\python27\lib\site-packages\markupsafe-0.23-py2.7-win32.egg
c:\progra~1\qgis2~1.16\apps\python27\lib\site-packages\python_dateutil-2.2-py2.7.egg
c:\progra~1\qgis2~1.16\apps\python27\lib\site-packages\pytz-2014.2-py2.7.egg
c:\progra~1\qgis2~1.16\apps\python27\lib\site-packages\shapely-1.2.18-py2.7-win32.egg
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:/Users/john/.qgis2//python
C:\Users\john\.qgis2\python\plugins\gearthview\ext-libs


texts: Grass Processing tools via QGIS : v.clean snap, chdangle, and rmdangle

texts:
Using QGIS 2.6 Grass Processing tools, I'm trying to clean up a topology on a water network layer.

I noticed topological errors when first trying to break all lines on intersections using v.clean break. Some didn't break, which makes me realize that some are not completely touching.

Which is the best to use: v.clean snap, v.clean chdangle, or v.clean rmdangle?

Lines seem to be moving for all cases tested, but I cannot tell whether the endpoints are extending towards the other feature or the feature is moving to meet the endpoint. I'm looking to do the former, if possible.

I've looked at the GRASS manual, but am still unsure:

  The chdangle tool is similar to the rmdangle tool but works only on
  boundaries and changes dangling boundaries to lines instead of
  removing them.

I guess I'm unable to visualise this.


texts:
i'm using leaflet to render a map on my site.

The problem is, i use mongoDB to perform some geospatial queries serverside, which use as standard the infamous epsg4326.

It appears that on leaflet there is support for this format, though i wasn't able to find a tile server with this specification.

Can anyone suggest one ?


texts: Filtering and updating manually selected features in MapInfo?

texts:
I'm new to MapInfo having spent most of my time with ArcGIS and am slowly getting used to how to do various common tasks.

Is there a way to filter/show only those features you have manually selected from the map window in the browser window and then make group edits to only those selected records (akin to the "show only selected features" and "field calculator" functions in the attribute table in ArcGIS)?


texts:
I would like to use Google's S2 library to convert geospatial data (for instance shapefiles or geojson files) to S2 cell IDs.
The data I want to convert is made of polygons.
I saw there exists this repository which enables the use of the S2 library in Python, but the S2Polygon object is not implemented.
Otherwise, I found this very nice tool which converts shapefiles to S2 IDs, but can't make it work on files other than those given in the example, it might be due to my files.
I have also found a Go implementation of the library, but it doesn't fully implement the S2Polygon object either.

Is there such a tool or function, other than the C++ S2 library ?


texts:
I need some help on how to go about it. I do not have much knowledge about GIS.

For a project we need a complete set of city data to perform calculations on.

In terms of generating a map of the city, is there a way that I can get the details of the city - i.e. all streets, landmarks, lanes etc. etc. in any format? Once I have these data, then I can insert them into the database and proceed.


texts: Finding duplicates of same point in feature class and only keep maximum of field using ArcPy?

texts:
I have a feature class with fields: Latitude, Longitude, Type, Magnitude, City, State, Comment.

If Latitude, Longitude , and Type are the same, I only want to keep the record with the highest Magnitude value (magnitude will vary usually).

However, sometimes there are duplicates in the data (all fields have duplicate values), so I would only want to keep one record.

I am running ArcGIS 10.3 for Desktop and Python 2.7.8.

After creating a table of max magnitude would I then join this table back to my feature class to create a new feature class?

arcpy.Statistics_analysis("C:/csv/all_reports.gdb/reportsSel", "C:/csv/all_reports.gdb/reportsSel2", [["MAG", "MAX"]], ["Type", "Longitude", "Latitude"])


texts:
I am comparing 2 lists of census data for the state of IL.

One from the FCC.GOV
where I download all the census tracts in CSV format and the data looks like this

state,county,cnamelong,tract,tractname,block,tractcode,blockcode
17,001,Adams County,000100,1,1000,17001000100,170010001001000
17,001,Adams County,000100,1,1001,17001000100,170010001001001
17,001,Adams County,000100,1,1002,17001000100,170010001001002
...

The 2nd comes from CENSUS.GOV which comes as a KML file. A sample of the data looks like this

GEOID
170119649001

I need to cross reference the data, but there is a problem.
The CSV tractcode is 11 digits, and the FCC geocode is 12 digits.
The first 11 digits are the same.
This site says the last number is the "block group" number. I don't see where that comes from. I need it in the CSV.

My questions are:

What data is apparent for me to cross reference these 2 data sources?
Where can I find the "block group" number where I can append it to the tractcode in the CSV for use of cross-reference?
Can I just use the first 12 digits of the "blockcode" in the CSV to compare?


texts: Insert Dataframe data into a map using Geopandas

texts:
I posted my question in stackoverflow previously, howeverI think this site is a better place to find some clues (this is my first project with spacial data).

I need to insert some coordinates from a dataframe to a NYC map. For this task I am using geopandas. This is my code:

gdf = gpd.GeoDataFrame(DF, geometry=gpd.points_from_xy(DF.pickup_longitude, DF.pickup_latitude))

I got the following:

 Unnamed: 0  DOLocationID  ...  time                    geometry
0           0           NaN  ...     0.380278  POINT (-74.01576 40.71555)
1           1           NaN  ...     0.390833  POINT (-73.98205 40.76276)
2           2           NaN  ...     0.350556  POINT (-73.97941 40.75543)
3           3           NaN  ...     0.503056  POINT (-73.77679 40.64509)
4           4           NaN  ...     0.618333  POINT (-73.78806 40.64167)

[5 rows x 31 columns]

Then I create the map using this code:

First

gdf.crs = 4326
gdf = gdf.to_crs(3857)

Then:

df = gpd.read_file(gpd.datasets.get_path('nybb'))
df = df.to_crs(epsg=3857)
ax = df.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
ctx.add_basemap(ax)

As you see, NYC becomes into the entire continent.

My question is: how do I plot the geometry column into the map accurately?


texts: Exception error when enabling workspace specific services in Geoserver 2.5 OpenGeo suite

texts:
When I try and enable a workspace specific WMS in the bundled GeoServer app from boundless I get the following error.

java.lang.reflect.InvocationTargetException

This is running on OSX 10.9.5 with java version

java version "1.6.0_65"
Java(TM) SE Runtime Environment (build 1.6.0_65-b14-462-11M4609)
Java HotSpot(TM) 64-Bit Server VM (build 20.65-b04-462, mixed mode)

The WMS works in global mode.


texts: Algorithm for all intersecting geometry parts

texts:
I'm trying to understand the ArcGIS "union" method so I can replicate it with an algorithm using Python and Shapely.

I do not want the cumulative intersection of all features where all of them intersect at the same time (ABC):

(Image is from another SO post)

Instead, I want each sub intersection that intersects at least once, repeated how many times it occurs. So in the example image I want ABC returned three times, and AB, AC, and BC returned twice. While the ArcGis method also returns the unique parts A, B, and C that have no intersections, I do not need this.

This is a pretty standard GIS method also available in QGIS, so what is the accepted set theoretic way of achieving this?

Code

In the ArcGis docs for union, they write that they use a "cracks and clusters" method, which is likely a low level modification of the intersection/union clipping algorithm. Since I am not implementing the underlying clipping algorithm, I need a different approach.

Here is pseudocode that I will be trying to implement. It basically means for each feature i find all intersections with others, and then recursively find all intersections between the intersecting parts, adding only the "node" parts that have no intersections with other intersections (i.e. we have cut it up as much as possible):

def isecs(g, geoms):
    for og in geoms:
        if og != g:
            if og.crosses(g):
                yield g.intersection(og)

def process(isecs):
    parts = []
    for g in isecs:
        isecs = getisecs(g, isecs)
        if no isecs:
            # "node" reached
            parts += g
        else:
            # this is the recursive part
            parts += process(isecs)
    return parts

for f in features:
    top_isecs = isecs(f.geom, features.bbox_overlaps(f.bbox))
    parts = process(top_isecs)
    for g in parts:
        addfeat(g)

Correction

@Vince correctly pointed out that the ArcGis operation is called union, instead of intersection as I first wrote. This is not to be confused that we are looking for geometrical unions, ArcGis is simply referring to the fact that they are returning all  (hence union) geometrical intersections: "Union calculates the geometric intersection of any number of feature class..." They also allow returning an actual geometrical union through a dissolve option.


texts: Extracting column only from R spatial dataframe without $ notation

texts:
I have a SpatialPolygonsDataFrame in R, from which I would like to extract a vector of one of the columns. This can be done easily by using the $ notation, e.g. dataFrame$columnName. My code however needs to pass the column names as strings, in which case I would normally use the bracket notation dataFrame[,'columnName'] to get the same result. R-spatial however, per the documentation, doesn't seem to allow this, returning instead a full SpatialPolygonsDataFrame object with geometry added to the single requested column.

How can I go from a string of the column name to just a vector of the data in that column? I'd prefer to do this without creating a separate non-spatial dataframe if possible.


texts: Reading ephemeris data from GPS receiver?

texts:
I want to filter the GPS data. when I use the low cost GPS receiver, it gives erroneous data. therefore I want to apply some filtering technique.

How can I read ephemeris data from GPS receiver?

I need to use ephemeris data for filtering.


texts: Keeping custom tool active while using edit tools via ArcObjects?

texts:
I'm developing a tool that, when the user creates a feature in a certain feature class, prompts the user with a form to populate one of the fields in its attribute table. I'd like to have the user use the Editor sketch tool to create the feature, rather than creating my own implementation of the sketch tool, and the form should appear only when the tool is active.

How do I keep my custom tool activated and allow the user to use the sketch tool at the same time?

Normally, a tool is deactivated when you switch to using another tool.

Following Preventing user from selecting different tool when custom tool is selected?, I looked at the ITool.Deactivate method, but according to the ArcObjects documentation, setting Deactivate to False makes it so that the user cannot select another tool.

I'm using and developing against ArcGIS 10.2.


texts: Slow performance when using workspace in URL

texts:
I have one workspaces with 30.000 raster layers in each one.
This request take in average 600ms
https://myserver/geoserver/WORKSPACE_A/wms?SERVICE=WMS&VERSION=1.1.0

This one take in average 180ms
https://myserver/geoserver/wms?SERVICE=WMS&VERSION=1.1.0....

The response is exactly the same image, same size, same interpolation. The only difference is that the slower one I pass the workspace in the route and the other one no.
Any idea what can be happening?


texts: Getting adjacent polygons to merge in QGIS?

texts:
I was given a shapefile made up of a few thousand squares.  The data in the squares has 5 different values.  I wanted to merge all the common values into the same polygon to simplify things.  This is what it looked like.

I used the the Select by Expression tool to select all the polygons of the same values and then used the Merge Selected Features tool to merge the polygons together.  I was hoping I would be left with 5 geometries that held each value.  However after the merge I had vertical stripes through each of the zones as shown here.

I suspect that I'm looking at polygons that have points that are not exactly at the same coordinates.  I also suspect that it is a problem created by the software program that generated the shapefile.  I was able to confirm this by using the Snap tool to move the points and snap it back over the other one and then I could merge the polygons together.  I don't have the time nor the patience to do that for every point on the map.

I've tried to use the GRASS clean tools to get rid of these points as well as a few other tools I could find and have searched the forum for solutions, but after an afternoon of not getting anywhere, I've decided to see if anyone could help me figure this out.  I've provided the (Shapefile Data) here if someone wants to try to do it.

How do I clean this up?

I'm running 2.18.14.


texts: Forcing "snap-to-road" when reshaping a line via basemap reference in ArcGIS desktop?

texts:
Sometimes, when I use "reshape a feature" in concert with a basemap, I'll find that the resultant line doesn't automatically snap to roads if a pair of verticies are far apart:

But at the same time, having to place closely spaced vertices along twisty roads is a pain.

Is there a balance that can be struck when it comes to ArcGIS's snap-to-road tolerance? Or do I just have to suck it up and put in closely-spaced vertices as the only way to guarantee that the resultant line will trace the path that I want?


texts:
Has anyone tried to store nation-wide lat/long ranges as a look-up table?

For example: If you had a polygon dataset that covered the entire US (for example, county polygons). And you wanted to look up any lat/long combination and return which county it is located within. Wouldn't it be possible to break the county data into say, 50m x 50m cells and assign lat/long ranges to each to make a giant matrix of all lat/longs in the US with assigned county attributes?

Has this been done?


texts: Why is iSpatialFilter in an "Add-in" crashing ArcMap on Partial or Full refresh?

texts:
I am building a tool in my Add-in that is partially based on a .Net sample and partially on my former VBA tool. (Both work)

My problem is that I would like to be able to only select features from the current layer in the TOC. For this i have to use iFeatureSelection:SelectFeatures with iSpatialFilter. This works flawlessly in VBA but not in my .Net Add-in

ArcMap crashes when I refresh / partial refresh the map.

Any ideas?

pGeoFeatureLayer.Selectable = True

pFeatureSelection = pGeoFeatureLayer

pEnv = m_pFeedbackEnv.Stop()

With pSpatialFilter
    .GeometryField = "Shape"
    .Geometry = pEnv
    .SpatialRel = ESRI.ArcGIS.Geodatabase.esriSpatialRelEnum.esriSpatialRelIntersects
End With

If Not arg.Shift Then
    pFeatureSelection.SelectFeatures(pSpatialFilter, ESRI.ArcGIS.Carto.esriSelectionResultEnum.esriSelectionResultNew, False)
Else
    pFeatureSelection.SelectFeatures(pSpatialFilter, ESRI.ArcGIS.Carto.esriSelectionResultEnum.esriSelectionResultAdd, False)
End If

m_focusMap.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, Nothing, Nothing)


texts: Can ogr2ogr read from PostgreSQL materialized view?

texts:
I'm writing a quick Python script to bulk export from PosgtreSQL materialized view (PostGIS geometry) to ESRI FileGDB and ogr2ogr (I'm using the one in OSGeo4W through Python subprocess module) is throwing an error

ERROR 1: ERROR:  relation "my_schema_name.my_TABLE_name" does not exist

The name is just an example. I want to underline that the view has upper letters, but I'm passing the argument double quoted like this table_name = '"my_schema_name.my_TABLE_name"'. Is that correct?

I have no problem exporting from normal tables that the view is based on.

Does ogr2ogr read materialized view at all?


texts: Add layer in Openlayers using layer variable passed to a function with setVisibility

texts:
I am trying to add and remove layers by passing the name of the layer variable to two separate functions and I can't seem to find any place that describes what I am trying to do. Here is an abridged version of my code:

function init(){
        var options = {
            projection: new OpenLayers.Projection("EPSG:900913"),
            displayProjection: new OpenLayers.Projection("EPSG:4326"),
            units: "m"
        };

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

        var drainage = new OpenLayers.Layer.Vector("Drainage Basins", {
            extractAttributes: true,
            protocol: new OpenLayers.Protocol.HTTP({
                url: "drainage_basins.gml",
                format: new OpenLayers.Format.GML()
            }),
            strategies: [new OpenLayers.Strategy.Fixed()]
        });

        map.addLayers([drainage]);
}

I would like to call showLayer(drainage); to turn the layer on but I am not what to do:

function showLayer(layer){
    layer.setVisibility(true);
}

and hideLayer(drainage); to turn the layer off as:

function hideLayer(layer){
    layer.setVisibility(false);
}

Thanks for your help.


texts:
I would like to apply an inverted polygon style, including opacity and inner shadow using PyQGIS.
So far, everything works except for the inner shadow effect. I found in the QGIS API that there is a class named QgsShadowEffect, but I do not get it to work.

The code looks like this:

#Open layer
poly_1_link = r"...\Polygon_0.shp"
poly_1 = iface.addVectorLayer(poly_1_link, "", "ogr")

##Edit layer style##
#Define layer renderer
rend = poly_1.renderer()
print(rend.type())

#Define symbol from renderer
some_symbol = rend.symbol()

#Define symbology
some_symbol.setColor(QColor.fromRgb(255,0,0))   #Colour
some_symbol.setOpacity(0.5)                     #Opacity

#Apply shadow
#...

#Create an inverted polygon
sub_rend = QgsInvertedPolygonRenderer().convertFromRenderer(rend)
poly_1.setRenderer(sub_rend)

#Update layer
poly_1.triggerRepaint()


texts:
I know the coordinates below are in UTM 36N (due to the EPSG code), and I can convert them with several websites online to wgs 84 and see them in Google Maps. However, been helplessly looking the internet for the conversion function itself, but found none.

<EPSGCode>32636</EPSGCode>
<Eastings>750400</Eastings>
<Northings>3000500</Northings>

Any idea how to convert those to WGS 84?


texts: Creating feature class in file geodatabase in ArcGIS Desktop with arcgisscripting?

texts:
I am using ArcGIS Desktop 9.3.

What I am looking to do is create a file geodatabase and add a feature class directly into that file geodatabase, as opposed to creating a feature class outside the file geodatabase and then importing it into the file geodatabase.  The reason I need to be able to create the feature class directly in a file geodatabase is because I need to have attribute table columns that are nullable, whereas if I create a feature class outside the file geodatabase and then import it into the file geodatabase, nulls will not be allowed. Please correct me if I am wrong.  Currently, I am able to create the file geodatabase, but I can't add a feature class directly to the file geodatabase.
The code I currently have is:

 import os, arcgisscripting
 gp = arcgisscripting.create(9.3)
 gp.Overwriteoutput = 1
 cwd = "C:\\Path\\To\\directory"
 gp.workspace = cwd
 gp.toolbox = "management"
 number = '1'
 fileGDB = "test_%s.gdb" % number
 shpFile = "file_%s.shp" % number
 gp.CreateFileGDB(cwd, fileGDB)
 gp.CreateFeatureclass(fileGDB, shpFile, "POINT", '#', '#', '#', '#')
 gp.addfield ( shpFile, fieldName, "FLOAT", "20","20", "#", "#", "NULLABLE", "#", "#")
 ...the rest of my code here


texts: ArcGIS Exporting (ALL) layers from Table of Contents

texts:
I'm trying to export ALL layers from the table of contents from a map document but I get the following error:

Runtime error  Traceback (most recent call last):   File "", line 14,
   in    File "c:\program files
   (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 2434, in
   CopyFeatures     raise e ExecuteError: ERROR 000210: Cannot create output
   C:\Users\jnmiller\Desktop\Global_Map.mxd\New_Shapefile.shp Failed to execute
   (CopyFeatures).

Python Code:

import os, arcpy

infolder = "C:\Users\jnmiller\Desktop\Global_Map.mxd"
outfolder = "C:\EnterpriseFolder\NEO"
#Get name of vector layers in the TOC
mxd = arcpy.mapping.MapDocument("current")
df = arcpy.mapping.ListDataFrames(mxd)[0]
layers = [f.name for f in arcpy.mapping.ListLayers(mxd, "*", df) if
f.isFeatureLayer]

for layer in layers:
    #Join output infolder to layer name and append .shp
    outfolder = os.path.join(infolder, "{}.shp".format(layer))
    arcpy.CopyFeatures_management(layer, outfolder)


texts: Show coordinates in other projections (ArcGIS Javascript API)

texts:
Is it possible to change this sample that it will be show x,y coordinates in wkid:5514 ?

<meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
<title>Create Map Display Mouse Coordinates</title>
<link rel="stylesheet" href="https://js.arcgis.com/3.15/esri/css/esri.css">

<script src="https://js.arcgis.com/3.15/"></script>
<script>
  var map;
  require([
    "esri/map", "esri/geometry/webMercatorUtils", "dojo/dom",
    "dojo/domReady!"
  ], function(
    Map, webMercatorUtils, dom
  ) {
    map = new Map("map", {
      basemap: "streets",
      center: [-47.109, 14.945],
      zoom: 2
    });
    map.on("load", function() {
      //after map loads, connect to listen to mouse move & drag events
      map.on("mouse-move", showCoordinates);
      map.on("mouse-drag", showCoordinates);
    });

    function showCoordinates(evt) {
      //the map is in web mercator but display coordinates in geographic (lat, long)
      var mp = webMercatorUtils.webMercatorToGeographic(evt.mapPoint);
      //display mouse coordinates
      dom.byId("info").innerHTML = mp.x.toFixed(3) + ", " + mp.y.toFixed(3);
    }
  });
</script>


texts:
When I try to install OpenLayers plugin for QGIS 2.2.0 Valmiera I get the error message
the plugin is broken; python said:
'QPyNULLVARIANT'object has no attribute '_getitem'

Does anyone have any ideas; i am a complete beginner, so I hope for a simple fix.
Or, where else can I get base maps.
I am using mac OS X 10.9.2

Thanks.


texts: Merging multiple tables to shape using for loop with ArcPy?

texts:
I have to join multiple tables to a shapefile using a common field. I want to script this in python (i am using ArcGIS 10.1). Below is the script i am using.

I am calling the tables one by one and merging them to the shape. I would like to put this in a for loop so the script reads all the tables (that have prefix dTable) and merge them to the shape.

#Redefine the workspace to the geodatabase

env.workspace = os.path.dirname(os.path.dirname(__file__))+'\\'+ out_name +'.gdb'
print "Workspace for GDB:"+ env.workspace

################################################################################
#STEP 5: Start Joining Tables
# Set local variables

featureName = "temp_Streets"
layerName="Streets2"

outFeature = "temp_Streets2"
expression = "temp_Streets.FUNC_CLASS < 99"

################################################################################
# Exclude centriod connectors from the shapefile before joining
try:
     # Make a layer from the feature class
    arcpy.MakeFeatureLayer_management(featureName,layerName)
    print "Converting to Layer Done !"

    # Select desired features from street layer
    arcpy.SelectLayerByAttribute_management(layerName, "NEW_SELECTION", expression)
    print "Selection Done !"

    # Copy the layer to a new permanent feature class
    arcpy.CopyFeatures_management(layerName, outFeature)
    print "Copy to new layer Done !"

    print "Step 5: Creating Streets file Done !"
    print "#################################################"
    arcpy.AddMessage("Creating Streets file Done !")
except:
    print arcpy.GetMessages()
################################################################################
joinTable1 = "dTable_HBW_csv"
layername1="j_Streets1"
outFeature1 = "temp_j_Streets1"
joinField = "CONNECT"

# Join the first CSV file
try:
     # Make a layer from the feature class
    arcpy.MakeFeatureLayer_management(outFeature,layername1)
    print "Converting to layername1 Done !"

    # Join the feature layer to a table
    arcpy.AddJoin_management(layername1, joinField, joinTable1, joinField, "KEEP_ALL")
    print "Join 1 Done !"

    # Copy the layer to a new permanent feature class
    arcpy.CopyFeatures_management(layername1, outFeature1)
    print "Copy to new layer 1 Done !"

    print "Step 5a: Joining Streets 1 file Done !"
    print "#################################################"
    arcpy.AddMessage("Joining Streets 1 file Done !")

except:
    print arcpy.GetMessages()
###############################################################################
joinTable2 = "dTable_HBO_csv"
layername2="j_Streets2"
outFeature2 = "temp_j_Streets2"
joinField2 = joinTable1+"_CONNECT"

# Join Second CSV file
try:
     # Make a layer from the feature class
    arcpy.MakeFeatureLayer_management(outFeature1,layername2)
    print "Converting to layername2 Done !"

    # Join the feature layer to a table
    arcpy.AddJoin_management(layername2, joinField2, joinTable2, joinField, "KEEP_ALL")
    print "Join 2 Done !"

    # Copy the layer to a new permanent feature class
    arcpy.CopyFeatures_management(layername2, outFeature2)
    print "Copy to new layer 2 Done !"

    print "Step 5b: Joining Streets 2 file Done !"
    print "#################################################"
    arcpy.AddMessage("Joining Streets 2 file Done !")

except:
    print arcpy.GetMessages()


texts: Combine two parallel lines to one with postgis

texts:
In OpenStreetMap, a highway=motorway consists out of two parallel lines, each representing one direction:

On a professional map the two lanes get combined to one single line and drawn with a middle line:

Is there an easy way with postgis functions to combine that two lines to one in between to be able to render the highway accordingly? It has to be considered that the two lanes could split somewhere and every direction has to be rendered separately again (e.g. for going around a mountain). So there needs to be a parameter to set the maximum distance between both lines.


texts: What is correct term for map with shaded country outlines?

texts:
Many world maps, especially old ones, often have shaded country outlines in a particular color. The saturation becomes weaker towards the inside of the country and is strongest at its borders. Here is an example of what I mean:

Is there a correct term for such a map style?

This style does not just appear in historical maps, but for example also in this National Geographic map of Europe from 1962.


texts:
i want to know how to name maps, put an id or something to find it later on openlayers, for example:

  div id="map1">

  div id="map2">

adition all code to set a map for fill it like:

  map1 = new ol.Map({
          id: 'Maps 1',
          target: 'map1' });

  map2 = new ol.Map({
          id: 'Maps 1',
          target: 'map2' });

then i could find the div for map 1 or 2 with

  mapDiv1 = document.getElementbyId('map1');

  mapDiv2 = document.getElementbyId('map2');

but how can i Access the map in mapDiv1 o mapDiv2 ?

also i would like how to do this:

  $.getJSON(theme.url, null, function (data)

but with pure JavaScript, that is couse data only exist inside this function, even if i asign it to a global variable it gets me null, also i don't use jquery or other frameworks just JavaScript.


texts:
Can I use Google Maps as a layer in ArcMap?

Is it possible?


texts: Geodjango PostgreSQL/PostGIS database question

texts:
I have a simple django app using sqlite3. My app is simple I have django user authentication system where the user can be upload spatial raster images and can see that images in user personal HTML page filter by use id.

Here is my model :

class nomodel(models.Model):
    name = models.TextField()
    user = models.ForeignKey(User, unique=True)
    upload = models.ImageField(upload_to=user_directory_path)

sqlite3 database :

DATABASES = {
    'default': {
        'ENGINE': 'django.db.backends.sqlite3',
        'NAME': os.path.join(BASE_DIR, 'db.sqlite3'),
    }
}

But that way is not GIS way to work with leaflet to add this raster in open layers.

If I change my code to :

from django.contrib.gis.db import models

class testRaster(models.Model):
    raster_id = models.TextField(primary_key=True, default=1)
    name = models.TextField(default='test_layer')
    raster = models.RasterField()
    #user = models.ForeignKey(User, unique=True)
    #raster = models.RasterField(upload_to=user_directory_path)

DATABASES = {
    'default': {
         'ENGINE': 'django.contrib.gis.db.backends.postgis',
         'NAME': 'geodjango',
         'USER': 'geo',
    },
}

but my question is because never use PostGIS in the past is in sqlite3 in upload field save the image path from my nomodel model but in the raster  field from testRaster model save again the path of images ?

in testRaster ant PostGIS database can I use again user authentication like sqlite3 model ?


texts: Optimum Path finding (for Cross Country Movement) using TIN

texts:
I am working on path finding algorithm (using A*) for cross country movement (not for road movement) on map using grid based approach. Algorithm supports path calculation considering slope of terrain.

This algorithm is based on A* grid approach and it expands more in case of long distance with obstacle and highly dense non uniform cost weighted region. I have optimized heuristic and type of data structure of my open and close list and applied many other techniques. But as I said it is taking long time as it is working on grid. I can't increase the resolution of the grid from certain limit as area is having dense non uniform cost.

I Googled about the issue and found that if I can use navigation mesh instead of grid then it will be faster. I came across this paper which proposed a solution of using a Triangulated Irregular Network (TIN) for finding the shortest path. I am currently studying this approach.

Is there anyone who had already worked on the solution proposed in the paper for finding shortest path using TIN model. Will it be useful for my application?

Is there any other technique you think will be useful for my algorithm?


texts: Get a list of leaflet marker clusters within current map bounds in Leaflet

texts:
I am using leaflet.markercluster plugin to show clusters on my map and getting information about each cluster on hover/click event.
What I'd like to do now is to write a method that if a user zooms or pans the map, it returns list of clusters within the current map bounds.
var option = {
    useCanvas: true,
    anchorWidth: 20,
    anchorHeight: 40,
    canvasPadding: 0.0,
    map: map
};
var clusterMarkers = L.markerClusterGroup(options, { chunkedLoading: true });

for(var i = 0;i<pakData.length;i++){
  option.lat = pakData[i][0];
  option.lng = pakData[i][1];
var marker = TPLMaps.overlays.point(option);
clusterMarkers.addLayer(marker);

}
map.addLayer(clusterMarkers);
clusterMarkers.addEventListener('clusterclick', (e)=>{
    console.log(e); // I want list of these e.target within current bounds
})


texts: Python script to read coordinate system of an input feature class?

texts:
How to read the coordinate system of an input feature class in python script.
How we can use "Describe()" function and tool validator for that purpose


texts: QGIS imported geotiff rendered as perspective view, how to change to top

texts:
I've created a GeoTiff that should have GPS coordinates in it using AgiSoft.  When I load it into QGIS, it renders the image in a perspective and I want to always see the image from the "top" down.  When I split the GeoTiff into bands, I get a top down view of the raster.

Is this an issue with how I exported the GeoTiff or is there a way to tell QGIS to "view from top" (change the projection?).


texts: Open source/free software for generating a 3D view of mapping?

texts: Is there any software that will generate a 3D view given a DTM/height data and raster mapping, and output an image in a printable resolution? I am looking for something similar to what you'd find in a commercial hike planning application, but free or open source preferably.


texts: How to add values of a mask image as one of the classification classes

texts:
I have performed supervised classification on GEE in the Russian Lena Delta to map different vegetation classes using sentinel 2 images (BUT i had to mask out the water and sandbank pixels) to improve the classification results.
However, I would like now (after the classification) to include the sand and water masks as additional classes in the classified layer (for visualization reasons) but somehow i cannot achieve this. Does anyone have any idea how?
beloow is my work flow
var roiDelta = ee.FeatureCollection("*********"),
    HC_Heim = ee.Image("*******"),
    geoms =
    /* color: #98ff00 */
    /* shown: false */
    ee.Geometry({
      "type": "GeometryCollection",
      "geometries": [
        {
          "type": "Point",
          "coordinates": [
            126.4953564453125,
            72.3816136993482
          ]
        },
        {
          "type": "Polygon",
          "coordinates": [
            [
              [
                126.09710205078125,
                72.19780053250948
              ],
              [
                126.676630859375,
                72.20787438492015
              ],
              [
                126.654658203125,
                72.44551315285538
              ],
              [
                126.061396484375,
                72.43639841991575
              ]
            ]
          ],
          "geodesic": true,
          "evenOdd": true
        }
      ],
      "coordinates": []
    });

var roi = geoms.buffer(25000); // to test on smaller subset
// var roi = roiDelta;

// cloud maks
function cloudMask(image){
  var QA60 = image.select('QA60');
  return image.updateMask(QA60.lt(1));
}

// compute ndvi across images in collection
function getNDVI(image){
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi)
}

// load Sentinel2 data for summer 2018
var S2dat = ee.ImageCollection("COPERNICUS/S2")
        .filterDate('2018-06-01', '2018-09-15')
        .filterBounds(roi)
        .map(cloudMask).map(getNDVI);

//Global surface water
var gsw = ee.Image('JRC/GSW1_0/GlobalSurfaceWater')
            .select('occurrence')
            .unmask(0);

var WaterMask = gsw.lt(60);

var ndviComp  = S2dat.qualityMosaic("NDVI").clip(roi);
var SandMask  = ndviComp.select("NDVI").gt(0.21).updateMask(WaterMask);

var ndviPalette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];

var ndviCompMask = ndviComp.updateMask(WaterMask).updateMask(SandMask);

//classification

var samplepoints = ee.FeatureCollection.randomPoints(geoms, 2000, 0, 10);
//Map.addLayer(samplepoints, {}, 'sample points', false);

var tablePts = HC_Heim.reduceRegions(samplepoints, ee.Reducer.first(), 10);
var tablePtsID = ee.FeatureCollection(tablePts.map(function(feat){
  var id = ee.Feature(feat).get('id');
  return ee.Feature(feat).set("Class_ID", id);
}));

// This property of the table stores the land cover labels.
var bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7','B8','B9','B10','B11','B12', 'NDVI'];
var label = 'first';

// Overlay the points on the imagery to get training.
var training = ndviCompMask.select(bands).sampleRegions({
  'collection': tablePts,
  'properties': [label],
  'scale': 10
});
//print(training);

var trained = ee.Classifier.smileRandomForest(10).train(training, label, bands);
print(trained);

var classified = ndviCompMask.select(bands).classify(trained).clip(roi);
// print(classified);

//Map.addLayer(classified,
             {min: 0, max: 11, palette: [
                  '68ab5f', '1c5f2c', 'a3cc51', '82ba9e', 'ccb879',
                  'dcd939', 'b5c58f', 'b8d9eb', 'af963c', '6c9fb8',
                  'f19414', '202cff']}, 'classes', false);

/// legend


texts: gdalwarp cutline without using a shapefile

texts:
I'm currently working with gdalwarp to split some tiffs.

gdalwarp -cutline file.shp -dstalpha -of GTiff image.tif split_image.tif

I wanted to know if I can use an other source rather than use a shapefile, for example I'd like to use a WKT entry, something like POLYGON ((0 0, 0 1, 1 1 , 1 0, 0 0))


texts: Installing qgis 2.0 and saga 2.1 on ubuntu 13.04

texts:
I have been trying to get this working for a few weeks now. I am trying to install QGIS Dufor (2.0.1) and Saga GIS from the Unstable branch of the UbuntuGIS PPA. When I attempt to do this I get the following error

The following packages have unmet dependencies:
 saga : Depends: libsaga (= 2.0.8+dfsg-4~raring) but 2.1.0-0~raring2 is to be installed
E: Unable to correct problems, you have held broken packages.

It appears that the dependencies have not been updated. When I try to install saga from the Johanvdw/saga-gis ppa i am told that two dependencies are too old and required by QGIS. The only way I can get both programs installed at the same time is to install the QGIS nightly build. When I do things this way Sextante can not seem to locate the saga routines. With the changes in sextante now that it has been added to trunk, there is no longer a field to specify the directory that SAGA resides in. I have a similar topic open HERE. I am out of ideas at this point.


texts:
ArcMap edit annotation tool button grey out?
I try to edit annotation, but can't because the select annotation button was grey out see image below

My annotation  was properly convert from label, store in Map.

Why I can't edit annotation? I was be able to edit annotation before, but I can't this time, why?


texts:
Is there a way to colorimetrically and automatically vectorize georeferenced rasters that came from manually colored, scanned maps? Currently we Print the  topography maps 20k scale --> hand-coloring work (hazard mapping) --> georeferencing (perfectly, if I may say so) ---(current stage)--> digitization of manually colored map.

If a solution exists i know some eye-to-screen adjustment would still be in order though I believe this path would be easier in the long run.


texts: Query Postgis to convert UTM to Lat Long in decimal degrees

texts: I want the Postgres / PostGIS queries to convert the UTM coordinates of the SIRGAS 2000 / UTM zone 20N (EPSG: 31974) projection to the Decimal Degrees projection. Separating longitude and latitude into distinct columns.


texts:
Despite having

<title>Threatened Ecosystems of South Africa</title>

in the html file, the Open Layers map does not include a title like Leaflet does. Is there some way around this, or is this just how it is unfortunately?


texts: Is it possible to create an addin for a specifc MXD?

texts: I am currently in the process of creating Addin's and also converting our VBA projects to .NET addins. With the VBA applications, they would only work with specific MXD's. Is this possible to do with Addin's as well?