texts: how to display simple arrow on the map in R using Leaflet

texts:
I try to display an  arrow on the map using R and the leaflet package,

code of map see my related question:

# install packages
install.packages(c("rgdal", "maptools", "leaflet", "htmlwidgets"), dependencies = TRUE)

# load libraries
library("rgdal")
library("maptools")
library("leaflet")
library("htmlwidgets")

# load france as shape
FRA=readShapePoly("data/france.shp")

# Initialsation
m <- leaflet(padding = 0)

# Add country
m <- addPolygons(map = m, data = FRA, opacity = 100,
                 color = "#FAFCFA",
                 weight = 0.25,popup = NULL,
                 options = list(clickable = FALSE),
                 fill = T, fillColor = "#B3C4B3",
                 fillOpacity = 100)

# Dimention of the map
m$width <- 874
m$height <- 700

# Export as HTML file
saveWidget(m, 'mapfrance.html', selfcontained = FALSE)

code of arrow :

x<-c(2.484722,-8.445833)
y<-c(51.5085300,10.50306)
i <- order(x, y); x <- x[i]; y <- y[i]
## draw arrows from point to point :
s <- seq(length(x)-1)  # one shorter than data
plot(FRA)
points(x,y)
arrows(x[s], y[s], x[s+1], y[s+1], col = 1:3,angle = 20)

how can I display a simple arrow on the map?


texts: How to change data type in Excel file before creating feature class?

texts:
I am trying to import data from Excel spreadsheets where each point coordinate has a 'datetime' timestamp (e.g. "2012-06-01 14:45:32") associated to it. When I check the table properties ArcGIS 10 always recognizes the timestamp field as 'date' (i.e. "2012-06-01") and drops the time-of-the-day info.

Is there a way to import the timestamp field as a full datetime type? And if yes, how do I go about it?


texts: Troubleshooting unsupported error message from ArcGIS Pro

texts: I keep receiving an error message for trying to upload a temperature raster layer onto ArcGIS Pro. I was wondering what would be the best course of action to troubleshoot?


texts: CartoDB only georeferences half of my table

texts:
I imported a CSV file into CartoDB, with over a million points, all with a clear (and correct) latitude and longitude. For many of the collected data points the georeference column says ‘null’, while the latitude & longitude are correct.

So CartoDB only georeferenced half of my table and I am missing a lot of points.

I do have a payed account and did not reach my maximum of points I can georeference.

Anybody know how I can get my whole CSV georeferenced?


texts:
I followed Mike Bostock's 'Let's Make a Map' tutorial to make a GeoJSON map of the EU starting from a 1:10M cultural admin0 shapefile from Natural Earth. The problem that I'm facing if that some of the countries are split into multiple features, i.e.:

{ "type": "Feature", "properties": { "SOVEREIGNT": "Italy", "ADMIN": "Italy", "ADM0_A3": "ITA", "SUBUNIT": "Sardinia" },
{ "type": "Feature", "properties": { "SOVEREIGNT": "Italy", "ADMIN": "Italy", "ADM0_A3": "ITA", "SUBUNIT": "Isole Pelagie" },
{ "type": "Feature", "properties": { "SOVEREIGNT": "Italy", "ADMIN": "Italy", "ADM0_A3": "ITA", "SUBUNIT": "Pantelleria" },
{ "type": "Feature", "properties": { "SOVEREIGNT": "Italy", "ADMIN": "Italy", "ADM0_A3": "ITA", "SUBUNIT": "Italy" },
{ "type": "Feature", "properties": { "SOVEREIGNT": "Italy", "ADMIN": "Italy", "ADM0_A3": "ITA", "SUBUNIT": "Sicily" }

N.B. I stripped most properties from the above the example.

The command I used to generate the GeoJSON looked something like:

ogr2ogr   -f GeoJSON   -where "ADM0_A3 IN ('PRT', 'ESP', 'ITA', ...)"   myMap.json   ne_10m_admin_0_map_subunits.shp

My question is:

Is there a way that I can have all the feature geometries that have a common ADMIN/SOVEREIGNT merged so that the GeoJSON output will only contain a single entry per country?


texts:
I'm am having some trouble plotting my spatial data using ggplot2. The map looks fine when plotted using spplot, so I'm assuming that the tearing occurs at the fortify stage.

The code is as follows:

#install the packages
library(rgdal)
library(mapproj)
library(raster)
library(rgeos)
library(ggplot2)
library(plyr)

if (!require(gpclib)) install.packages("gpclib", type="source")
gpclibPermit()

setwd("C:/Users/My Documents")

#read in laa to regional mapping
#must aggregate to higher level regions as data is provided at this higher level
laa_region_mapping <- read.csv("laa_region.csv", header = TRUE)

#read in LAA polygons
laa_polygons <- readOGR("ctyua_ew_generalised_WGS84.json", "OGRGeoJSON")

#merge by laa to add region column to polygon data
laa_polygons_with_region_data <- merge(laa_polygons, laa_region_mapping,
                                by.x = "CTYUA13NM", by.y = "LAA",
                                all.x = TRUE, all.y = TRUE)

# aggregate laa polygons by the 21 regions (aggregate by regoin_code)
region_polygons <- raster::aggregate(laa_polygons_with_region_data, "region_code")

The aggregate has worked, as can be seen by the spplot (note: I found how to aggregate by the regions from this SE post: Join spatial polygons by code in R)

#plot the resulting polygons using spplot
spplot(region_polygons)

But when I fortify the spatial data so that I can use ggplot, there is tearing around the edges.

#fortify and merge to create the data frame ggplot will show on the map
region_polygons@data$id <- rownames(region_polygons@data)
region_polygons.points <- fortify(region_polygons, region = "id")

# plot the fortified df using ggplot
ggplot(data = region_polygons.points, aes(x= long, y = lat, group = id, fill=id)) + geom_polygon()

How can I stop this tearing?

I've looked at similar responses on SE, but the responses suggest that tearing occurs during a merge (What is the cause of 'tearing' of polygons (artifacts) using R, ggplot and geom_polygon?). I think my tearing occurs at the fortify stage as the spplot before fortifying looks fine.


texts: QGIS, add a layer (shapefile) to a running QGIS from an external application

texts:
I'm trying to load a shapefile as a new layer into a running QGIS application using python within a Windows command shell.  The python is basically what is explained in this post:

How to open a layer in QGIS from another application?

I have worked through various errors and now I get none, but the layer does not appear in my QGIS.  And I get no feedback.  I can drag the shapefile onto the map and it renders fine.  I'm using QGIS 2.0.1.

I don't know if I'm doing something fundamentally wrong.  Can I communicate with a running QGIS using the mechanism described?


texts: How to avoid having the layer name in the field name?

texts:
I'm working with ArcMap 10. I might be doing something wrong but can't see what:
I want to add a field in an attribute table. Problem is that as I name it "FieldName" the resulting name in the header is "LayerName.FieldName".
Why? How to avoid the layer name to be added in the field name?


texts: How to make band composite image in QGIS

texts: I am wondering how I can make a band composite from three large bands files from Landsat? I was using option Raster>Miscellaneous>Merge but when I try to do it it crashes.


texts: Want to dissolve adjacent polygons but no common edges

texts:
I have a bunch of parcel data, I want to join the polygons basing on the same owner name.
Even they are not touched, I want merge them together, such a 100 feet tolerance. The purpose is to show the somebody's properties without too much separated lines.
Like this
I have some programming experience, but just start to try to get the tool by python script, and tried arcgis add_in wizard.
If there is some idea about the script work flow, it will be really helpful!


texts: How can i include program form (created in Designer ) to my QGIS plugin?

texts:
How can I include a form created in Designer in my QGIS plugin?

And how to get controls values (textbox,label,etc)?


texts: googleearthengine - Code to calculate EVI having rendering issues

texts:
Below is the link of my code to calculate EVI in Google Earth engine. It is not giving/rendering color according to the color palette! Please help me in this regard.
https://code.earthengine.google.com/995504963211f6c4d300ac364dc6697c

function maskS2clouds(image) {

  var qa = image.select('QA60');

  var cloudBitMask = ee.Number(2).pow(10).int();
  var cirrusBitMask = ee.Number(2).pow(11).int();

  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
}

var composite = s2.filterDate('2018-03-01', '2018-03-15')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .filterBounds(ft)
                  .map(maskS2clouds)
                  //.multiply(100)
                  .median();

Map.setCenter(87.32,22.40, 6);
Map.addLayer(composite, {bands: ['B6', 'B4', 'B3'], min: 0, max: 0.3},'fcc');

var landcover_roi = composite.clip(ft);
var evi =landcover_roi.expression(
    '2.5 * ((NIR-RED) / (NIR +6 * RED -7.5* BLUE))', {
      'NIR':landcover_roi.select('B8'),
      'RED':landcover_roi.select('B4'),
      'BLUE':landcover_roi.select('B2')
    });
Map.addLayer(ft, {}, 'From Fusion Table');

var EVI=evi.multiply(100).uint8()
var meanDictionary = EVI.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: ft.geometry(),
  scale: 20,
  maxPixels: 1e9
});

print(meanDictionary)

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

// Center the map
Map.setCenter(87.32,22.40, 6);

print(EVI);
 Map.addLayer(EVI,{min: -1, max: 1, palette: palette},'Sentinel-2 evi')
Export.image.toDrive({
  image:EVI,
  description: 'evi',
  scale: 30,
 region: ft,
  maxPixels:34089663741
  });


texts: How to create multiple dockable windows?

texts:
I am building an add-in and I have several different windows, which I open with the code below. This works fine. However, how can I create the same window n times?

    UID dockWinID = new UIDClass();
    dockWinID.Value = ThisAddIn.IDs.MyWindow;

    IDockableWindow dockWindow = ArcMap.DockableWindowManager.GetDockableWindow(dockWinID);

    dockWindow.Show(true);


texts: ArcMap 10.5: Point layer's custom symbols not all exporting with geoprocessing

texts:
I am using a Geoprocessing Service to export (package and send email to download) an mxd with it's sourced file geodatabase. One of the point layers in my mxd has custom symbols. Once I export my files and download the mxd with it's gdb, sometimes my point layer has all the custom symbols, and sometimes it only has a few. The features of my layer do exist after I download, but the custom symbols which are images are lost.

Is there a specific way to save the custom symbol images, so they are not lost upon exporting or is there any reason why the symbol images would be lost while being processed with the Geoprocessing service?


texts:
Hopefully someone can help me figure this issue out. I am trying to system call gdal functions inside of R using the system / system2 functions. Here is an example of the code:

> system("gdalinfo --version")
sh: gdalinfo: command not found

> system2("gdalinfo --version")
sh: gdalinfo --version: command not found

As you can see system cannot call the gdal functions. Anyone know how to fix this?

Inside of the Mac terminal, I can call any of the gdal command line functions:

user-MacBook-Pro:~ user$ gdalinfo --version
GDAL 1.10.1, released 2013/08/26

Here is a Sys.info -- I am running Yosemite on Mac 10.10.1:

Sys.info()
# sysname        "Darwin"
# release        "14.0.0"
# version        "Darwin Kernel Version 14.0.0: Fri Sep 19 00:26:44 PDT 2014; root:xnu-2782.1.97~2/RELEASE_X86_64"
# nodename       "user-MacBook-Pro.local"
# machine        "x86_64"
# login          "user"
# user           "user"
# effective_user "user"

Also see the difference between these two calls. First in R:

  system("echo $PATH")
  /usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin

And now in terminal, where GDAL is explicitly mentioned:

  user-Pro:~ user$ echo $PATH
  /Library/Frameworks/GDAL.framework/Programs:/Library/Frameworks/GDAL.framework/Programs:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin


texts: Cannot bring WFS layers from GeoServer to OpenLayers

texts:
So I found this course on Udemy and I made it to the part where we start adding vector layers into an OpenLayers map. However, the layer we add does not appear on my map! The code for this is here:

wfs_url = 'http://localhost:8080/geoserver/richland_postgres/wfs?service=WFS&' +
          'version=1.1.0' +
          '&request=GetFeature' +
          '&typename=richland_postgres:schools&' +
          '&CQL_FILTER=' + cqlfilter + '&' +
          'outputFormat=application/json&srsname=EPSG:3857&' +
          ',EPSG:3857';
        prompt('',wfs_url);

        var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: function(extent) {
      return 'http://localhost:8080/geoserver/richland_postgres/wfs?service=WFS&' +
          'version=1.1.0' +
          '&request=GetFeature' +
          '&typename=richland_postgres:schools&' +
          '&CQL_FILTER=' + cqlfilter + '&' +
          'outputFormat=application/json' +
          '&srsname=EPSG:3857&,EPSG:3857';
    },
    strategy: ol.loadingstrategy.bbox
  });

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

I copied and pasted the URL from the prompt into GeoServer and it gave me this:

I think that is a GeoServer issue and I am about to have it re-installed anyway. However, I took out the CQL from my code and while the URL worked perfectly, however the layer still doesn't show up.

Looking around this site I did find another student from this class who had the same problem and someone suggested adding a BBOX. I tried that myself and it didn't work. I don't know whether it worked for them. I've tried a variety of tutorials but I have never been able to get a WFS into OpenLayers, however I can get it into QGIS. I wonder if this is a GeoServer issue.


texts: Convert a KML ImageOverlay to World Files with GDAL

texts:
I have many raster outputs from an Open Source radio transmission topographic path tool called SPLAT, georeferenced in KML image overlay format.

I have figured out an easy way to use tools to convert from KML Image Overlays to PGW World Files using this small Java tool (which has fallen off the internet, I rebuilt it from Subversion repo):

WorldFileTool
0.3.7 jar download link, my personal build

Is there a way to use GDAL to make this conversion more script-able?

If I were to implement a simple script to write out world files would there be any gotchas?

I realized I had originally referenced Create world file from kml a few weeks ago, don't miss that one.


texts: Create FeatureClass in Personal Geodatabase and store data

texts:
I have to create a new featureClass in my existing personal GeoDatabase.
And save some polylines into this featureClass.

How can i do it?

I need c# code to do this.

Edit: I have a problem during Geometry type featureClass creation.
Here the code i'm using:

IFields fields = objectClassDescription.RequiredFields;
ESRI.ArcGIS.Geodatabase.IFieldsEdit fieldsEdit = ESRI.ArcGIS.Geodatabase.IFieldsEdit)fields; // Explicit Cast
ESRI.ArcGIS.Geodatabase.IField field = new ESRI.ArcGIS.Geodatabase.Field();
ESRI.ArcGIS.Geodatabase.IFieldEdit fieldEdit = (ESRI.ArcGIS.Geodatabase.IFieldEdit)field;
fieldEdit.Name_2 = "Shape";
fieldEdit.Type_2 = ESRI.ArcGIS.Geodatabase.esriFieldType.esriFieldTypeGeometry;
fieldEdit.IsNullable_2 = true;
fieldEdit.AliasName_2 = "SampleShap";
fieldEdit.DefaultValue_2 = "test";
fieldsEdit.AddField(field);
fields = (ESRI.ArcGIS.Geodatabase.IFields)fieldsEdit;
ESRI.ArcGIS.Geodatabase.IFieldChecker fieldChecker = new ESRI.ArcGIS.Geodatabase.FieldChecker();
ESRI.ArcGIS.Geodatabase.IEnumFieldError enumFieldError = null;
ESRI.ArcGIS.Geodatabase.IFields validatedFields = null;
fieldChecker.ValidateWorkspace = (ESRI.ArcGIS.Geodatabase.IWorkspace)workspace;
fieldChecker.Validate(fields, out enumFieldError, out validatedFields);
IFeatureClass featureClass = featureWorkspace.CreateFeatureClass("Test01", validatedFields, null, null, ESRI.ArcGIS.Geodatabase.esriFeatureType.esriFTSimple, strShapeField, "");

I get this error:

The module of database management can't find the table or the input query 'GDB_DataChanges'

How can i solve it?


texts: ArcGIS Buffer tool leaves out zero width buffers?

texts:
I have a polygon shapefile with 93 objects in it. They have a column with values for the buffer distance. 90 of these objects have 5 as a value, while the remaining 3 have 0.

When I use the buffer tool, the shapefile created leave out the 3 records containing the 0-values, so only a shapefile with 90 objects are in it.

Why is this? I have managed to do this before with 0-values, is there something I have forgotten?


texts: What are the limitations of using MySQL with ArcGIS?

texts:
I recently connected a MySQL database to ArcGIS by adding an OLE DB Connection. I can view the tables just fine; however, my XY data will not show up on the map directly from the MySQL table, I have to export the data to a standalone File Geodatabase table first, and then display the XY Data. Is this normal?

It seems like it sort of defeats the purpose of dynamic data tables. I wanted to directly manipulate the MySQL tables in ArcGIS so that changes I make to the attributes in ArcGIS are reflected on the MySQL database and vice versa.

This brings me to my real question. I was unable to locate any detailed documentation concerning:

What are the limitations of using MySQL through an OLE DB Connection with ArcGIS?

I found an inkling on the limitations from this post.

Establishing a MySQL connection in ArcGIS 10
". . geometries from the Spatial extension will not come through to the GIS."

Can anyone tell me where I can find out more? Should I just switch to Microsoft SQL for use with ArcGIS?


texts: feature store get NullReferenceException error

texts:
I am trying to update a field value according to the selected OBJECTID. But get System.NullReferenceException error. Don't know what's wrong with the code.

My code is as follows:

        pActiveView = pMxDoc.FocusMap
        pMxDoc = My.ArcMap.Application.Document
        pMap = pMxDoc.FocusMap

        Dim layerNum = GetIndexNumberFromLayerName(pActiveView, "dev_db.DBO.Segment")

        Dim pFLayer As IFeatureLayer
        pFLayer = pMap.Layer(layerNum)

        Dim featureClass As ESRI.ArcGIS.Geodatabase.IFeatureClass = pFLayer.FeatureClass

        Dim feature As IFeature = pFLayer.GetFeature(1)
        Dim nameFieldIndex As Integer = pFLayer.FindField("Seg_ID")
        feature.Value(nameFieldIndex) = "12/11"
        feature.Store()

The method is called in an editing session, and confirm has the OBJECTID 1 in the specified layer.

I am using arcobjects 10.2.2 and vb.net arcmap add-in.


texts: Great circle rendering of satellite ground tracks reversed in southern hemisphere

texts:
I'm working on rendering some satellite ground tracks from postgis, through geoserver, and onto openlayers. The python library I'm using to get lat, lon from TLE data isn't perfect, so there are distortions in the lines.

To smooth the lines out, I wanted to use the great circle projection solution that I've used before. https://anitagraser.com/2011/08/20/visualizing-global-connections/

Unfortunately, the great circle lines drawn in the southern hemisphere are being projected in the wrong direction(I think). Is the 953027 projection only valid in the northern hemisphere?

The data is currently stored as a multilinestring in the DB, and the reprojection is being done by geoserver in a view.

Any ideas on how I can correct the great circle projection in the southern hemisphere?

EDIT: Here's some additional info if it helps.

The multilinestrings I'm working with are comprised of 180(one per minute for 3 hrs) points. The density of the points can be adjusted fairly easily, just takes a bit to run the script. They span from -180 to 180 and are split into multi to avoid the wrapping issue. This is an example if you want to look at one. https://pastebin.com/nfmRDkwm

My transform :

select  ST_Transform(ST_Segmentize(ST_Transform(ST_Force2D(sat_lines::geometry), 953027), 100000), 4326) as sat_lines LIMIT 5

953027 Projection

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 953027, 'esri', 53027, '+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=60 +lat_2=60 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs ', 'PROJCS["Sphere_Equidistant_Conic",GEOGCS["GCS_Sphere",DATUM["Not_specified_based_on_Authalic_Sphere",SPHEROID["Sphere",6371000,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Equidistant_Conic"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",60],PARAMETER["Standard_Parallel_2",60],PARAMETER["Latitude_Of_Origin",0],UNIT["Meter",1],AUTHORITY["EPSG","53027"]]');


texts:
I am trying to implement a method for performing a task at regular intervals in QGIS. This is to be triggered by a timer (currently 'threading.Timer', although it could be something else) and should utilise a thread separate from the GUI as well as maintain communication with the user to allow for feedback and for cancelling of the task.

I am basing my method on this QGIS multithreading tutorial. It seems an excellent source of information but it is slightly out of date (i.e. PyQt4, not PyQt5) and assumes a single action not a recurring one.

What I have got so far is below. The main thing that confuses me here is the double use of threading. I don't presumably want to create a thread which creates another thread. Since the Worker class handles the threading does that mean that the timer should operate on the original thread and therefore 'threading.Timer' is not the most appropriate option?

My objective is to get this working as 'template' code which I can then implement to do something useful.

(And if there is a totally different way of doing it, of course, I am open to suggestions. )

from qgis.core import *
from PyQt5 import QtCore, QtGui, QtWidgets
import traceback
import time
import threading
from PyQt5.QtCore import QTimer

class Worker():

    def __init__(self):
        self.killed = False

    def run(self):
        ret = False
        try:

            # Do something that takes ten seconds...
            for i in range(0,10): # loop 0 to 9
                time.sleep(1) # wait for a second
                self.progress.emit(100 * i/9) # update progress bar

            if self.killed is False:
                self.progress.emit(100)
                ret = True # successfully completed
        except Exception as e:
            # forward the exception upstream
            #self.error.emit(e, traceback.format_exc())
            pass
        self.finished.emit(ret)

    def kill(self):
        self.killed = True

    finished = QtCore.pyqtSignal(object)
    error = QtCore.pyqtSignal(Exception, basestring)
    progress = QtCore.pyqtSignal(float)

# end class

def workerError(self, e, exception_string):
    QgsMessageLog.logMessage('Worker thread raised an exception:\n'.format(exception_string), level=QgsMessageLog.CRITICAL)

def startWorker(self):

    #import QgsMessageLog
    QgsMessageLog.logMessage('About to create Worker')

    # create a new worker instance
    worker = Worker()

    # configure the QgsMessageBar
    messageBar = self.iface.messageBar().createMessage('Doing something...', )
    progressBar = QtWidgets.QProgressBar()
    progressBar.setAlignment(QtCore.Qt.AlignLeft|QtCore.Qt.AlignVCenter)
    cancelButton = QtWidgets.QPushButton()
    cancelButton.setText('Cancel')
    cancelButton.clicked.connect(worker.kill)
    messageBar.layout().addWidget(progressBar)
    messageBar.layout().addWidget(cancelButton)
    self.iface.messageBar().pushWidget(messageBar, Qgis.Info)
    self.messageBar = messageBar

    # start the worker in a new thread
    thread = QtCore.QThread()
    worker.moveToThread(thread)
    worker.finished.connect(workerFinished)
    worker.error.connect(workerError)
    worker.progress.connect(progressBar.setValue)
    thread.started.connect(worker.run)
    thread.start()
    self.thread = thread
    self.worker = worker

def workerFinished(self, ret):
    # clean up the worker and thread
    self.worker.deleteLater()
    self.thread.quit()
    self.thread.wait()
    self.thread.deleteLater()
    # remove widget from message bar
    self.iface.messageBar().popWidget(self.messageBar)
    if ret: # ret = True - completed successfully
        self.iface.messageBar().pushMessage('Task completed')
    else:
        # notify the user that something went wrong
        self.iface.messageBar().pushMessage('Something went wrong! See the message log for more information.', level=QgsMessageBar.CRITICAL, duration=3)

QTimer.singleShot(60, lambda: startWorker(qgis.utils)) # execute every 60 seconds


texts:
I have a shapefile with 52K points of road accidents, between 2007 and 2012.

Each accident has a severity index. I want show clusters of 5 accidents (at least), for each year, that have a severity index equal or above 20, and that are in a 100 meters radius of the same road.

I tried already create buffers of 100 meters and intersect the buffers, but I was unable.

I'm working on QGIS 2.0.1.

Can you help me?


texts: Getting Latitude & longitude of a Mollweide projection

texts:
I am very new to geotools and map-projections and after researching a hole day I still don't know where to start.

I have a map with a size of X,Y Pixels and it is in Mollweide projection.
How can I manage in Javacode to find the long/latitude of a specific point ( X1,X2 ) in my map.

Can anyone help with a starting hint, or the best would be a code example ( link )

Thanks for mentoining the 'Java Map Projection Library'. I tried to do the transformation with JMPL, but i guess I have a lack of understanding. How can I tell the transformation, how is my map sized ( X,Y ). Or how shall I do the transformation ? My Code so far:

    Point2D.Double pointonmap = null;
Point2D.Double latlon = null;
MolleweideProjection molproj=new MolleweideProjection();

pointonmap = new Point2D.Double (1400,1000);

latlon=molproj.inverseTransform(pointonmap,new Point2D.Double ());

System.out.println("latlon: " + latlon.getX() + ", " + latlon.getY());

best regards


texts: Extract multiple grids from HDF with subset tool using ArcPy

texts:
This should be fairly simple thing but bugging me at the moment. I have ~5k HDF files in a folder with names:

MOD11A1.A2001001.h24v05.005.2006343042921
and so on.....

I have following which loop through all of them to extract a subset grid from each HDF file. For unknown reason, python script breaks down every time at around file number 3557. I am able to manually extract several grids around that number from HDF files. Maybe, I am hitting RAM limit. But even with all other applications closed, it exhibits same behavior.

rasters = arcpy.ListRasters("*", "All")
for inLST in rasters:
    f = inLST[9:16]
    outLST = "D:/temp/subset/" + str(f) + '_lst'
    arcpy.ExtractSubDataset_management(inLST, outLST, "0")

Is there any alternate to arcpy.ListRasters, which I can use?


texts: How to create a feature class in memory?

texts:
I am trying to create a feature class in memory but I am not sure how to use the syntax.

My code looks like this:

FC = arcpy.CreateFeatureclass_management("", "in_memory/FC", "POLYGON", "", "DISABLED",
                                         "DISABLED", Coordinate_System, "", "0", "0", "0")

Basically I seem to have problem defining the "out_path".


texts: How to make a loop that interates over large raster cells in R go faster?

texts:
I am trying to add +1 to the cell value of certain cells in a raster. I have the cell numbers of the cells on which I'd like to perform said computation saved in a vector. I have a working code but with a large vector and raster it takes ages:

r <- raster(ncol=5,nrow=5)
r[] <- 0
vec <- c(1,1,1,3,4,5,6,7,8,8,8,9)
for(i in 1:length(vec))
      {
      r[vec[i]] <- r[vec[i]] + 1
      }

I have already tried to use r[vec] <- r[vec] + 1 but that doesn't work properly, it just sets the values of the cells stored in the vector to 1.

Does anybody know a faster way to do this ?


texts: When is it better to use 64 bit File GDB access?

texts:
The beta 2 release File Geodatabase API includes 64 bit support.

Esri says arcobjects is still the "preferred" method to interact with the file gdb.  However, arcobjects has no 64 bit support.

Assuming licensing is not a factor, in what situations would it be better to use the 64 bit file gdb api rather than arcobjects?


texts: Adding geometry from .shp file to .csv layer

texts:
Trying to get started with pyQgis and failing hard. My end goal is to take geometry attributes from a .shp file and add it to a .csv file in order to spatially represent details contained in the .csv.

To elaborate: I have a .csv of NYC permits pulled and a .shp of tax maps. I've created matching fields by manipulating Borough Block and Lot data. The trick is, joining is a many to one problem. There is only one tax lot and an unknown number of permits for that tax lot. I see easiest route as adding geometry fields to each permit entry.

I'm using the cookbook and have looked extensively through the questions here but the answers which seem to solve my problems have not.

Below is my code:

from ggis.core import *
import os
from PyQt4.QtCore import *
from PyQt4.QtGui import *

qgis=QgsApplication([], False)

#supply path to where your qgis is installed
qgis.setPrefixPath(os.environ["QGIS"]+"apps/qgis", True)

#supply path to qgis install location--same as above
#QgsApplication.setPrefixPath("C:\PROGRA~1\QGIS2~1.16\apps\qgis", True)

#load providers
qgis.initQgis()
taxlots_uri='file:\\\C:\Users\nd73\Desktop\Data\Clipped_Tax_Lots_11101b.shp'
taxlots_shp=QgsVectorLayer(taxlots_uri,'taxlots','ogr')

permits_uri='C:\Users\nd73\Desktop\Data\DOB_Permit_Issuance_Queens.csv'
permits_csv=QgsVectorLayer(permits_uri, 'permits', 'ogr')

print QgsApplication.showSettings()

print (QgsProviderRegistry.instance().providerList())

if not taxlots_shp.isValid():
  print "Layer1 failed to load!"

if not permits_csv.isValid():
  print "Layer2 failed to load!"

I have been able to run other very similar versions without error, but continues to print 'failed to load'. This version of the code incorporates a few lines which may be wrong. I was trying a few methods to troubleshoot.

I should add that my latest hurdle is that after running the code once qgis stops responding and I must close it and reopen to try again. Not sure what that's about.

From what I've read my problem must be in the paths I'm setting...so below are the paths.

Verbatim according to the properties of the qgis shortcut:
"C:\Program Files\QGIS 2.18\bin\nircmd.exe" exec hide C:\PROGRA~1\QGIS2~1.18\bin\qgis.bat

When I printed the prefix path through the python cmd line I got:
C:/PROGRA~1/QGIS2~1.16/apps/qgis

Any clarification would be great! Obviously I haven't even gotten to my iteration loops as I can't even get past the first step.


texts: Tools/processes that avoid the need for raster to vector converstion?

texts:
A "brute force" method of solving GIS problems that involve acquiring values from raster cells and storing them on spatially overlapping vector features is to first convert the raster to a vector. The main issue with this approach is that converting a raster to vector can lead to very large file sizes - and slow operations. Furthermore, if using ArcGIS without spatial analyst extension, raster to vector converstion is not available.

Could you provide a few of your own examples where you had a project and choose a work-path that purposely avoided a raster-vector conversions? This might have happened either because

i) you had to avoid the conversion because you did not have access to a raster to vector conversion tool, or

ii) you chose to avoid the conversion becuase doing so was a more elegant solution.

EDIT: (Here is an example I used in a previous version of my question before the comments below. It is an example I gave that is referred to in some answers).

Say you had a river network, and you wanted to store the elevation in each line segment in a new column - but elevation values in your GIS database is in raster format (a very large raster). What would be a workaround?


texts:
I am using the QNEAT plugin to calculate catchment are/travel reach over a street network. The network is valid and projected correctly. I have tried changing the cell size and topology tolerance without it making any difference. Running QGIS 3.10 on Mac.
My problem is that the interpolation process simplifies the output to such an extent it is not usable. The picture shows the network and the output polygon. As you can see the polygon fills up gaps in the network in a way that it should not. In this case, the polygon extends far out in the sea where no network is present.
Am I doing something wrong or do I not understand the functionality of the tool?


texts:
I have created a script that splits a shapefile into multiple new shapefiles and adds them to an ArcMap document. When running the program a second time to overwrite those shapefiles, lock-files are created and the program cannot overwrite the files.

What I do is I open the document, clear out all layers from that document, then I try to delete all the old shapefiles and then overwrite them, lastly add the new shapefiles to the document, save and close.

The strange thing is that this happens on some computers and not other. I wonder what arcpy methods could cause such locking and how I can clear the lock?

Update: I have now tried to set env.overwriteOutput = True but without result. The code below will simply print the error message. When I close the IDLE python console the lock files disappear.

def delete_month_shape_file(file_path, *args):

    if arcpy.Exists(file_path):
        try:
            arcpy.Delete_management(file_path)
        except arcpy.ExecuteError:
            message.error("Could not delete month shapefile.")
    else:
        pass


texts: Are there alternatives to the Google Places API, for lat/long retrieval?

texts:
I am looking for an alternative API to Google places both free and paid to extract location data (lat/long) for business locations i.e all of the "coffee shops in London".

I have looked at Factual, OSM and POIPlaza. They do not fit because of different reasons varying from pricing to completeness of data.

Anything suitable out there?


texts: Understanding a complex query in pgrouting

texts:
I am new to pgrouting and basically am familier with sql and postgresql, i came across the following command below:

SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
                SELECT gid AS id,
                         source::integer,
                         target::integer,
                         length::double precision AS cost
                        FROM ways',
                30, 60, false, false);

I am used to seeing sql statements as below:

SELECT param1 , param2

But i am thrown off the hook when i see something like:

SELECT seq, id1 AS node, id2 AS edge

What's with the AS part ?

Also is the below a subquery ?? and can somebody explain it ? I am not very clear with the below query.

(' SELECT gid AS id,
         source::integer,
         target::integer,
         length::double precision AS cost
        FROM ways',
30, 60, false, false);

Got the above query from HERE.


texts: Selecting zip codes without using lasso tool?

texts:
I'm new to ArcGIS Desktop.

I have a layer representing US zip codes and another layer of data points. I need to select all the data points within every zip code and find out the min-max value of the data points within that zip code. Currently, I've been selecting the zip code through the lasso tool and then performing a select by location but this is very time-consuming.

Is there any other way that I could possibly select all the data points zip code by zip code without having to use the lasso tool?


texts: Where can I find a high quality real globe of a decent size without axis?

texts:
I'm trying to find a physical globe of a high map quality that doesn't have any axis or bearings as the axis mounted on standard desk globes completely obscures the arctic regions.

High-quality map globes seem to be available but repeated searches failed me to find a globe model of a decent size (say 12'' and up) with an Atlas quality map and without axis, so that one can view the earth from any angle and the arctic regions are fully mapped and easily visible. Do you know whether/where such a product exists?

Edit: Magnetic globes won't do, as they are too small and the one example from an answer also obscures the polar regions through its magnets.

What I'd really like is a complete sphere without any knobs or visible seams (I don't mind seams if they are well covered) that can be placed in a holder. Inflatable globes seem quite OK, but their map-quality is very low as far as I could find.

Hmmm ... it seems these would go in the right direction (as far as I can trust the images): http://www.worldglobeu.com/movaglobes -- but they only go up to 8,5'' the maps look so-so and they are insanely expensive for the size and map-quality they offer.

This one seems to get the axis/seam thing halfway right, but the mapwork is more artistic that accurate: .

(Note that if that globe offered an Atlas quality map it's price tag might even be acceptable at 30''.)


texts: Adding units to numeric labels in ArcMap?

texts:
How do I add the distance units to labels in ArcMap?

I have the area calculated and would like the label to be displayed as (XXX sq ft).

I was able to write a VBScript of "(" & [AREA] &" )" but am not sure how to include the "sqft" inside the parenthesis.


texts: Neighborhood boundaries and building footprints for Johannesburg?

texts:
Are there any sources for free (or commercial) and distributable neighborhood and building footprint datasets for Johannesburg out there?

I've already looked on:

Google
OpenStreetMap
joburg.org.za
demarcation.org.za
sa-venues.com
SOCIOECONOMIC DATA AND APPLICATIONS CENTER (SEDAC)
geoscience.org.za
gadm.org

and probably some more that I do not remember right now...


texts: Error when assigning a spatial reference to a PointGeometry in arcpy

texts:
I'm creating a PointGeometry using arcpy:
data_path = "E:/Seafile/Seafile/Meine Bibliothek/02_Paper/02_ArcMap/Data/"
point = arcpy.Point(75.9039059275, 35.6928318018)
ptGeometry = arcpy.PointGeometry(point)
arcpy.CopyFeatures_management(ptGeometry, data_path + "point.shp")

So far, so good. However, if I try to add a spatial reference from a prj. file, it gives me this error  RuntimeError: Object: CreateObject cannot create geometry from inputs with this code:
Albers = arcpy.SpatialReference(data_path + "projection.prj")
point = arcpy.Point(75.9039059275, 35.6928318018)
ptGeometry = arcpy.PointGeometry(inputs=point,
                                 spatial_reference=Albers)

I don't get why this would not work. According to the ESRI website I'm following procedure. I don't necessarily need to use the .prj file if there is another way.
So far, I have tried using the Describe function to get the coordinate system from another file, but that didn't work, too. I'm looking for a way to assign a spatial reference to this PointGeometry. It is probably very easy but I'm lost.
EDIT:
The error also occurs if I just add inputs to the code as the parameter name:
ptGeometry = arcpy.PointGeometry(inputs = point)

So maybe the coordinate system is not the problem?
If necessary, this is the content of the prj. file:
PROJCS["Albers_Conic_Equal_Area",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",87.0],PARAMETER["standard_parallel_1",25.0],PARAMETER["standard_parallel_2",50.0],PARAMETER["latitude_of_origin",37.0],UNIT["Meter",1.0]]


texts:
I want to create a Python AddIn button to execute a specific script. I know how to use the python addIn wizard and add the toolbar I created in my ArcMap. The trouble I am having is writing the script itself.

I want my script to pull a specific table into the ArcMap document. Then, I want it to add a field to that table. Then, I want it to run a for loop that updates that new empty field based on the values of another field. I have this script already written for a specific table in a personal geodatabase that I just run with IDLE, but wanted to know how I could implement this as a AddIn button.

Below is the generic script that you edit before adding it as a toolbar.

import arcpy
import pythonaddins

class ButtonClass1(object):
    """Implementation for New folder_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        pass

This is what I have below. Right now I can't even get the button to show up, it just says missing. The code I have written below works as its own script, but trying to get it to do what I want as a button.

import arcpy
import pythonaddins

class Button(object):

"""Implementation for NewAddIn_addin.button (Button)"""
def __init__(self):
    self.enabled = True
    self.checked = False
def onClick(self):
    arcpy.env.workspace = r"C:\Users\sutton\Documents\Python\PythonExcercise\PlanoData.mdb"
    try:
        #create new field to hold new values
        arcpy.AddField_management("SmokeMain","yesorno","TEXT","10")
        print("Field added to SmokeMain")
        with arcpy.da.UpdateCursor("SmokeMain",["Surface_Cover","yesorno"]) as cursor:
            cntr = 1
            for row in cursor:
                #update the indicator field
                if row[0] == "CONCRETE":
                        row[1] = "YES"
                else:
                        row[1] = "NO"
                cursor.updateRow(row)
                print ("Record number {} updated".format(cntr))
                cntr = cntr + 1

    except Exception as e:
        print(e.message)

       pass


texts: Iterating over feature class and performing select by location using each feature as Selecting Features in ArcPy for ArcGIS Pro?

texts:
I have a polygon feature class which contains many (9952) polygons, some of which overlap each other.
I have another point feature class with many points (15998) and one of the fields is called route_number.

What I need to do is for every polygon in the polygon feature class, I need to select the points which intersect it. Then with these selected points I need to find the number of unique values for route_number. For example I can do this in ArcGIS Pro by selecting one polygon feature in the attribute table and then running the "Select Layer By Location" tool. It successfully selects just the points which intersect the polygon.

How do I do this in Python using ArcPy?

I have the following code to do a Select by location however it selects every point in the points feature class. I think this is because the select_features is the whole polygon feature class (polygons). I want this to be just one feature within the polygons feature class.

selection = arcpy.SelectLayerByLocation_management(in_layer='point_fc',overlap_type='INTERSECT',select_features='polygon_fc')
print(selection)
num_points = int(arcpy.GetCount_management(selection)[0])
print(num_points)


texts: Geographic Coordinate System to Projected Coordinate System

texts:
I am trying to convert a shapefile from GCS to PCS to calculate the area of each polygon.  I have tried two different methods:

1) "Data Management Tools->Projections and Transformations->Define Projection" and clearing the previous projection.  I then go back in and define the new PCS projection ("World_Cylindrical_Equal_Area").  This creates a tiny shapefile centered on 0,0 lat/lon with a miniscule amount of area.

2) "Data Management Tools->Projections and Transformations->Feature->Project" and choosing a new Output Coordinate System ("World_Cylindrical_Equal_Area").  When I run this, an error message comes up that says the table name is invalid and cannot be found.

I am trying to convert from "GCS_Bissel_1841" to "World_Cylindrical_Equal_Area" - any help would definitely be appreciated.

Also, here is a link to the shapefile: https://filestogeaux.lsu.edu/public/download.php?FILE=tjoyne1/71691ctWDJg


texts: 'undefined' is not a constructor (evaluating 'new > OpenLayers.Layer.GML'

texts:
My OpenLayers map was working fine before, showing a KML overlay. Suddenly it broke. Now I see this error in my browser error console.

  TypeError: 'undefined' is not a constructor (evaluating 'new
  OpenLayers.Layer.GML')

My code looks like this:

var kmllayer = new OpenLayers.Layer.GML("KML", "./my-kml-file.kml",
   {
    format: OpenLayers.Format.KML,
    formatOptions: {
      extractStyles: true,
      extractAttributes: true,
      maxDepth: 2
    }
   });
map.addLayer(kmllayer);


texts: Use dissolve tool on a polygon layer with attributes from a point layer. ArcGIS

texts:
I'm trying to find a way to automatically merge polygons with others based on a max + minimum value of points in each polygon. So the small ones go together with other small ones and the bigger ones stays the same.
Basically I've made this far with a Buffer around each point Dissolve and Minimum bounding geometry to get the polygons like this (see link in comments for more precise)  Now I want to take this to the next level and make the smaller polygons to merge with other small polygons. I'm curious if it's possible to do all this with tools instead of drawing by hand which could be easy in this case but I've got other much larger maps.
I want the result to look more like this (from Qgis and its drawn by hand)

Note that the lines has nothing to do with this it's all about the polygons around the points. I want about the same amount of points in each polygon.

One thought was to choose all polygons with < 50 points and make the process again in a new layer with larger buffer but it doesnt quite do the job.
(I hope this makes any sense to you, its hard to explain when I'm new to this)


texts: Rasterising shared shapefile using QGIS?

texts:
I am trying to rasterise this shapefile of ancient woodland across England. When I do it I just get a black Tif with out any positive values (just 0 for min and 0 for max). If I change the styling to singleband pseudcolour and manually input a max value of 1, QGIS (2.18.11) crashes.

I've done this in the past for viewsheds which had an integer field called DN and values of 1 which represented land with views to my observation point. I added an integer field in the ancient woodland table, giving each feature a value of 1 but the DN field does not appear in the drop down menu within the rasterisation window.

I need to do this for another 20 shapefiles so can anyone suggest a fix or correction to my method?


texts: Leaflet map url change on geojson click event

texts:
I am setting up a basic Leaflet Dashboard. This combines a Leaflet map view with a live table.

Demo: DAW Dashboard

I have been able to add a WMS service that I have made on Geoserver. This allows me to show project contours by adding de project "Id" in the CQL filter.

I have added the layer in de code:

var geototaal = L.tileLayer.wms("http://h2641061.stratoserver.net/geoserver/DAW_Suport/wms?&CQL_FILTER=Id%3D5", {
layers: 'DAW_Suport:geo_totaal',
attribution: 'Basemap DAW Project'
});

The CQL filter is now set to only get the contour of project with Id = 5

I also have a layer with point data (GEOJSON) that the dashboard uses to build the project marker locations and fill the table view. if I click on a project marker is opens a popup with project info (Id and Name). This also highlights the marker.

Code:

pointToLayer: function (feature, latlng) {
  if (feature.properties && feature.properties["marker-color"]) {
    markerColor = feature.properties["marker-color"];
} else {
  markerColor = "#FF0000";
}
return L.circleMarker(latlng, {
  radius: 4,
  weight: 2,
  fillColor: markerColor,
  color: markerColor,
  opacity: 1,
  fillOpacity: 1
});
},
onEachFeature: function (feature, layer) {
if (feature.properties) {
  layer.on({
    click: function (e) {
      identifyFeature(L.stamp(layer));
      highlightLayer.clearLayers();
      highlightLayer.addData(featureLayer.getLayer(L.stamp(layer)).toGeoJSON());

    },
    mouseover: function (e) {
      if (config.hoverProperty) {
        $(".info-control").html(feature.properties[config.hoverProperty]);
        $(".info-control").show();
      }
    },
    mouseout: function (e) {
      $(".info-control").hide();
    }
  });
}

}

I am hoping that some one can help me find a way that will change the URL of the WMS layer (CQL filter number) with a Id of the selected object so the project contour is shown in the map.

I have also made a zip package for download:

http://85.214.221.88/DAW_Support/development/development.rar


texts: How can I convert an image into map tiles for leafletJS

texts: I'm trying to make it so that people can upload an image onto my site and have it converted into tiles for leaflet that can be used for a custom interactive map. The images are maps from games, not real geographical data. I've been looking around and I cant find a way to slice up the images into the tiles that are usable in leaflet directly through JS. I've found plenty of standalone programs but they won't do the job in this situation.


texts:
Does anyone know where to find disease information in ESRI Shape vector format? I am trying to track zika, dengue, and flu in local communities, but any others are welcome too. My attempts so far are listed below

most of the "GIS" links on the CDC's website are broken
the WHO data is almost entirely pdf
NIH data is pdf and very broad, state/country wide generalizations
USGS social data is "retired"


texts:
I have 50 near tables that I want to append, so I can export them at once.

However, when I use the append tool an error comes up saying that the near tables are 'read only'. Is there anyway to make this append possible?


texts: Converting set of Shapefiles to OpenStreetMap (OSM) file?

texts: How can I convert a set of ShapeFiles to an OSM file?


texts: How to display direction of track with arrows using OpenLayers

texts:
I would like to display track on the map using OpenLayers in the way like this: http://jsfiddle.net/DR/hgCUq/1/light/

How to show direction of the track with arrows here? Any other approach would be great as well.

I googled https://github.com/Sirtrack/PointTrackWithArrows but this project looks buggy and conflicts with OpenLayers.Layer.PointTrack


texts: Cannot convert OneGeology KML to layer file in ArcGIS Desktop?

texts:
I downloaded and saved a KML file of Indonesian geology from OneGeology Portal (http://portal.onegeology.org/OnegeologyGlobal/). However, when I use the KML to Layerfile tool in ArcGIS, the tool does not run and gives me the following error:

  ERROR 000401: No features were found for processing Failed to execute
  (KMLToLayer).

How can I convert the file correctly. I ultimately need it as a shapefile.


texts: How to remove XML tags within a field of attribute table?

texts:
I have this dataset I got it from Wikimapia ,containing bunch of fields . One of fields named "PopupInfo" contains xml tags :

and looks like this when in editing mode :

I want to remove these tags and keep the text between them . is there tool or script can do this .


texts: How to "uncheck" the "No Data Value" box in Transparency on QGIS with Python?

texts:
The title question is overly simplistic, but would accomplish what I am looking for. I have generated a clipped *.tif raster using gdal.Warp and a georeferenced satellite image:
OutTile = gdal.Warp(newRasterLocation, rasterSourceFile, cutlineDSName=shapefileSourceFile, cropToCutline=True, dstNodata = 0)

When I open this raster in a simple image viewer, I get something like this (case A):

When I import this same raster into QGIS as a new layer, I see many regions where what I imagine would be centered around the darker spots, but not pure black values (or any one single color, case B):

If I uncheck the No Data Value indicator in the Transparency menu, the previously Case B image (with transparency holes) reverts to the Case A image.
What does the No Data Value check box control? I am thinking this is related to the dstNodata setting in the generation, but why would this affect more than one color? And why would toggling the No Data Value check box have any affect, if this was a property of the generated image?
(the image is a tree zoomed in, the dark spots are shadows. I am interested in this exclusively graphically, and don't plan on doing anything with the raster data other than showing the image).


texts: Drawing wind rose with QGIS from PostGIS?

texts:
I am working on wind rose data. I have some points on lines (streets) in my Postgres/PostGIS DB. The sample data scenario is:

First, I computed street orientation of this line using PostGIS commands ST_StartPoint(), ST_EndPoint() and ST_Azimuth() like this:

Select
-- line_dir = line direction (orientation) in degrees
(degrees(ST_Azimuth(line_end_pt, line_start_pt)))::int As line_dir,
-- assign the line direction as the lower bound of sector 1 (of 30 degree)
(degrees(ST_Azimuth(line_end_pt, line_start_pt)))::int As   lower_bound_sector1,
-- add 30 degree to lower bound to get the upper bound as each sector is of 30 degree max.
(degrees(ST_Azimuth(line_end_pt, line_start_pt)))::int + 30 As upper_bound

from
(
Select
     pt_on_line.geom as pt,
     line.geom As line,
     ST_StartPoint(ST_LineMerge(line.geom)) As line_start_pt,
     ST_EndPoint(ST_LineMerge(line.geom)) As line_end_pt
from pt_on_line
Left join line on ST_DWithin(pt_on_line.geom, line.geom, 0.1)
)
As foo;

Now starting with this street orientation, I need to generate a wind rose of 12 sectors (each of 30 degree having its lower bound and upper bounds) as geometry around this point on line in a radius of 50 meter. In the above code, I have assigned street direction as lower bound of sector_1 and added 30 degree to the lower bound. This approach might be unrealistic though. The expected output could be something similar to this.

How can I draw these 12 wind sectors (geometry) around this point starting from line's orientation?

Using the answer of @MichalZimmermann, here is the screen shot of wind rose generation using PostGIS and visualized in QGIS.


texts:
Recently, we set up a web app with OpenLayers that allows exporting a map as a jpeg image. The map contains XYZ, WMS and/or WFS layers.
We used the toJpeg function of the html-to-image library, as it was proposed in the OpenLayers "map export" example page a while ago, with a function similar to this:
const printMap = (mapPrintDiv) =>
toJpeg(mapPrintDiv, exportOptions)
   .then(dataURL => {
       try {
          const link = document.getElementById('imageDownload');
          link.href = dataURL;
          link.click();
       } catch (error) {
          console.log(error);
       }
})

Everything was working fine for a while but sometimes, the customer told us that this function was not working (the Promise of toJpeg was not fullfilled but the error message was not informative at all). Recently, the problem was more and more reported. It has been tested by dozens of people by the customer and by us, and it is EXTREMELY hard to see in what configuration this bug could be reproduced (no apparent effect of OS, browser, network, map layers, ...). Sometimes, it works and the next day it does not work.
Then, we saw that the "map export" example of OpenLayers actually changed 7 months ago: the example no longer relies on html-to-image but now use a custom function for printing (with the native canvas.toDataURL).
Why has it changed? What is the best method for exporting an OpenLayers map in 2020?


texts: Random forest classification starting from a JSON file

texts:
I have a lot of images of Earth and I need to find the section of the image where there is water using random forest.

I also have one JSON file (not GeoJSON) that have all the image's information. Inside the JSON file for each image there all the vertex points of the coordinates of the polygon (or polygons) that contains water in the image. The coordinates are related to the image and they are not geospatially referenced.

I have to train and test my network using random forest in order to find in a new image all the sections with water.

Do you have any suggestion how I can implement what I need in Python?
Do you know any website/github code that I can use as reference?


texts: Calculating shadow of trees with DEM in QGIS?

texts:
I have a DEM with a good precision (0,5 x 0,5 m/pix LiDAR) and I´m looking for the way to obtain a raster (or shapefile) with the shadow of the trees located in that place for a whole day (ideally winter solstice)

Is it possible with some tool to calculate it through QGIS?

My intention is to use this raster (or shape) with the maximum shadow to restrict the installation area of a solar power plant.

I have more information, like the height and canopy diameter of each tree.


texts: FME SpatialFilter not outputting desired results

texts:
Have this workflow that pretty much starts with a mapinfo table with either a point feature or a polygon feature that identifies land, licences, and leases, gets filtered based on feature type.  From here all points go to the spatial filter to be compared against a property layer and merge these attributes together.  I have this issue where some points are within the same polygon but the attributes are different (Ie; Local gov leasing an area from another authority and then subleasing to a community group therefore some attributes are different), and I want to replicate this showing two separate property polygon features.  You can see that 715 features go into the spatial filter and 708 come out.

For some background, this is a property register and I use points to get the property information since this comes from a state authority and they increase the accuracy of parcels when they receive more accurate information, so I have gone this method to be more dynamic so the boundaries in our land register change with the States Authority.  After I figure this i want to be able to compare "todays date" (when it gets run) with the expiry of leases and licences and email me so we know when they fall due for renewing.


texts: PostGIS - Intersects feature to split into linestring

texts:
I would like to know if it's possible to intersect feature datasets (multipolygon, polygon, multilinestring and linestring) to get a result as a set of linestrings being split at each intersection point.

Output should give:

4 linestrings for the square
3 linestrings for the line
2 linestrings for the polygon


texts: Is it possible to add layers onto a static image in Google Earth Engine Python API?

texts:
This is my first question on this site, so I apologize if my question's unclear or if it's an oversight on my part. I'm pretty new to Google Earth Engine and am still trying to learn as much as I can.
I'm trying to add an NDVI and Pollution layer onto a static image in GEE Python API but I haven't found any documentation for it. I understand this is possible using geehydro and Folium but the client that I ultimately want to display my results in do not support maps, only static images.
The closest I have come so far managed to produced an image (code below). However, the map is extremely warped and is unusable for analysis. Is there a way to "flatten" this image or create layered images more reliably?
dem = ee.Image('MODIS/006/MOD13A2/2018_01_17').select('NDVI')
AOI = ee.Geometry.Rectangle([-180,-75, 180, 75])
# Display a thumbnail of global elevation.
Image(url = dem.updateMask(dem.gt(0))
  .getThumbURL({'min': 0, 'max': 3000, 'dimensions': 512, 'palette': ['red', 'yellow', 'green'], 'region': AOI,
                }))


texts: What is the minimum / maximum of longitude / latitude?

texts:
I'm currently writing code which checks if coordinates are proper. I'm basically asking which check I need to do:
For longitude:

A: -180 <= LONGITUDE <= 180
B: -180 <= LONGITUDE < 180
C: -180 < LONGITUDE <= 180
D: -180 < LONGITUDE < 180

For LATITUDE:

a: -90 <= LATITUDE <= 90
b: -90 <= LATITUDE < 90
c: -90 < LATITUDE <= 90
d: -90 < LATITUDE < 90

Please give a reference or a good explanation.


texts: Displaying multiple categories on a QGIS map?

texts:
I'm tasked with displaying reserve data on a map but also displaying the category data for each site even though they overlap with other reserves. E.g. I have a category, "Birdwatching" which applies to 50 sites whilst the category, "Grassland" applies to 30 sites but overlaps with the birdwatching category thus a site could have "birdwatching and grasslands".

I created a layer for every category, so a user could select "birdwatching" and view the 50 of 100 reserves that meet this rule. But how do I go about composing a map that includes all ten categories?


texts: Selecting n-random points within polygon as single feature in Google Earth Engine?

texts:
IMPORTS:

var geometry = /* color: #98ff00 */ee.Geometry.Polygon(
    [[[73.85016678395664, 14.304113379426035],
      [73.85016583974448, 14.048667693939866],
      [74.0477686190103, 14.155139129529527]]]);

CODE:

var fg_points = ee.Feature(geometry);

I thought if the geometry is given as a feature it will select all points in that polygon but it isn't.

How to select some n number of Random Points within a polygon as a single feature in Google Earth Engine?


texts: Alternating Horizontal Labels using Maplex

texts:
I am using Maplex label engine to achieve horizontal lines from label to feature as illustrated by below example;

What I want is for the labels to alternate in direction as per the below example;

How can I achieve this?


texts:
I am updating annotations and the annotations feature class reference scale is at 18,000 to 20,000 while my data frame is set at 26,500. When I print the map, the font size seems to be significantly smaller.

Is there a way to increase the font size when it shows on a map?


texts: OpenLayers install npm install --save-dev parcel-bundler error

texts:
I am trying to setup OpenLayers on my computer.

I am following the directions here:
http://openlayers.org/en/latest/doc/tutorials/bundle.html

But when I do:
npm install --save-dev parcel-bundler
I get error:

npm WARN deprecated browserslist@1.7.7: Browserslist 2 could fail on reading Browserslist >3.0 config used in other tools.
npm WARN optional Skipping failed optional dependency /chokidar/fsevents:
npm WARN notsup Not compatible with your operating system or architecture: fsevents@1.2.4
npm WARN ol@1.0.0 No description
npm WARN ol@1.0.0 No repository field.

How do I fix this?

I have tried searching for this on google and duckduckgo but to no avail.

I have found the first error described in more detail here:https://www.npmjs.com/package/browserslist/v/1.5.1

Don't get how to solve it though.


texts: How to write a script that will automatically draw lines from one point to multiple destination points? (ArcMap 9.3)

texts:
I need to create lines from an airport to numerous destination locations. The destinations change often, so drawing them manually becomes tedious. I want to create a script that will draw them automatically. Does anyone know how to do this?

I've written a number of scripts for ArcGIS 10 using arcpy, but nothing for 9.3. If you could point me in the right direction that would be great.


texts: Change style of existing feature in OpenLayers2

texts:
I'm using OpenLayers 2 and have a number of LineStrings with custom attributes's drawn. I need to change the strokeColor on some of the lines styles programatically and have it stay this color even if the user change zoom level or toggles layer visibility.

Currently I'm filtering/selecting features with vectorLayer.getFeaturesByAttribute() and then I'm using vectorLayer.drawFeature(feature, { strokeColor: "red }) to change the line colors. This work well until the user change the zoom level or if I toggle layer visibility back and forth with vectorLayer.setVisibility(), which for some reason changes the strokeColor back to the defalt style.

What is the correct way of changing the feature style and make it "stick"?


texts:
I want to create Attribute-based polygon with SLD but when I use Persian or Arabic filter It's not working.
This is my sld:

    Attribute-based polygon

      SLD Cook Book: Attribute-based polygon

    <Rule>
      <Name>rule1</Name>
      <Title>rule1</Title>

      <ogc:Filter>
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>test</ogc:PropertyName>
          <ogc:Literal>'مقدار ب'</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>

      <PolygonSymbolizer>
        <Fill>
          <CssParameter name="fill">#66FF66</CssParameter>
        </Fill>
      </PolygonSymbolizer>
    </Rule>

    <Rule>
      <Name>rule2</Name>
      <Title>rule2</Title>

      <ogc:Filter>
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>test</ogc:PropertyName>
          <ogc:Literal>'مقدار الف'</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>

      <PolygonSymbolizer>
        <Fill>
          <CssParameter name="fill">#cc5577</CssParameter>
        </Fill>
      </PolygonSymbolizer>
    </Rule>

  </FeatureTypeStyle>
</UserStyle>

Can anyone help me?


texts: Explaining different results from Ordinary Kriging and Universal Kriging?

texts:
When using either methods to create a surface from the same set of observations, these surfaces differ from each other.

Visually it can be seen in the differences of the following contours (overlaid; the Red is Ordinary Kringin, Purple is Universal).

What are the main differences in those methods so that they would create different results?

The surfaces were produces using ArcMap's 3D analyst - Default parameters for both tools. I realize that not many people have access to those tools, so I've put online the produced surfaces here. If you want to recreate the contours, in your software use base heights of 9 and step of 2.


texts: PostgreSQL + PostGIS count with ST_DWithin multiple distance at once

texts:
I am using ST_DWithin to count rows that are within 5KM, What I want is count on multiple distance at once with single query i.e count for 5KM, 10KM, 20M
Current Query
SELECT count(*) AS cnt FROM table_name
WHERE table_name.published IS TRUE AND
ST_DWithin(ST_Transform(table_name.geo_point, 26986),
            ST_Transform(ST_SetSRID(ST_MakePoint(8.398522, 47.000084), 4326), 26986), 5000)

Here, 5000 =  5KM
What I Want
KM_5     KM_10   KM_20
100       215     300

Created Index
CREATE INDEX properties_geo_point_transform_idx ON
properties
USING gist (ST_Transform(geo_point, 26986)) WHERE
geo_point IS NOT NULL;


texts: Returning list of field names in joined table using ArcPy?

texts:
I need to look at all the fields in a joined table performed through arcpy. I am instigating

for fields in arcpy.Listfields(Feature):
  print fields

but it is returning a large list of:

<geoprocessing describe field object object at 0x0EBE4E60>
<geoprocessing describe field object object at 0x0EBE4D40>
<geoprocessing describe field object object at 0x0EBE4B18>
<geoprocessing describe field object object at 0x0EBE4758>
<geoprocessing describe field object object at 0x0EBE4E90>
<geoprocessing describe field object object at 0x0EBE4A28>
<geoprocessing describe field object object at 0x0EBE4428>

I have done this before but forget how its done.


texts: Force Web AppBuilder to use HTTP instead of HTTPS?

texts:
I have an out-of-the-box instance of ArcGIS Web AppBuilder developer edition - the startup.bat file is unmodified and includes the lines:

REM open the url in browser
START http://%FQDN%:%port%/%wabVirtualPath%
EXIT 0

However, when I start a new session on my PC, it runs as https://localhost:3344/webappbuilder - even if I forcibly change the URL to http:// it immediately reverts to https://

This is causing a problem when I try to access ArcGIS Server REST endpoints which are not served via https.

For example, when trying to add a sample layer via the Local Layer Widget I get an error message about insecure content:

I'm working in a corporate environment and don't have the option of serving my local layers via HTTPS, so I'd like to know whether it's possible to force Web AppBuilder to run on HTTP.


texts: Polygons of different layers into one layer. ArcMap

texts:
my problem addresses the union tool. I edited polygons in eight different layers because the area I investigate is quite big. Now I want to have all this edited polygons of these eight layers in one layer. The first problem is that when I use the union tool I have to choose a template which can only be one of the already existing. So I do not get a new one, but yes, this is one thing I could ignore. The more important thing is that when I use the tool the polygons do not exist for themselves anymore. They are ONE and I can not even pick a single one out of it. They do not occure in the attribut table. But I need every single of my created polygons listed with different FIDs because I want to do further calculations. I tried it with Feature to Polygon and this tool solves the problem BUT even worse is now that the main part of my already done work disappears. Most of the polygons are directly neighboured so I had to cut the polygons. And now after the union they are registered as one again so I would have to cut all these polygons again.

Is there someone who can help me handling this problem?

To sum it up again: I need to "union" the polygons of eight layers into one without resulting into neighbouring polygons melting together. They need to stay exactly as I edited them, everyone of them with a single FID.


texts: GEE: How to iterate through a feature collection to calculate annual average indices from MODIS for each feature

texts:
I am trying to calculate the average annual vegetation condition index from MODIS data for a set of raions (districts) for spring months (march, april and may). I am wanting to export this data to a table with this annual data for each individual raion. I have created this script which works for all of the raions as a whole, but I need to adapt the code to loop through the raions feature collection and calculate the average VCI for each raion each year in spring and summer months. I believe I need to use the map() function but my GEE/Javascript experience is minimal and am not sure how to integrate this into the code.
//calculate average seasonal NDVI for a set of raions
//get MODIS NDVI
var collection = ee.ImageCollection('MODIS/006/MOD13A1').select('NDVI');

// reference period for baseline
var reference = ee.ImageCollection(collection.filterDate('2001-01-01', '2020-12-31'))
  .sort('system:time_start');
print("reference",reference);

//filter out images from specific months
//march, april, may (spring)
var spring = reference.filter(ee.Filter.dayOfYear(59,150));

//filter bounds to raions
var spring_soa = spring.filterBounds(raions.geometry());

//calulate min of full series (mam 2001-2020) for each pixel for spring and summer seasons
var spring_min = spring.min();

//calculate max of full series (spr 2001-2020) for each pixel
var spring_max = spring.max();

//mask out agriculture
var spring_min_agrimask = spring_min.multiply(agrimask_ukr);
var spring_max_agrimask = spring_max.multiply(agrimask_ukr);

//calc bottom of VCI equations
var vci_bot_spr = spring_max_agriSoA.subtract(spring_min_agrimask);

//function to iterate through months and calculate average Vegetation condition Index (SPRING) FOR RAIONS
var years = ee.List.sequence(2001, 2020);
print("years",years);

var spring_vci = ee.ImageCollection.fromImages(
      years.map(function (y) {
        return spring.filter(ee.Filter.calendarRange(y, y, 'year'))
                    .select(("NDVI")).mean()
                    .multiply(agrimask_UKR)
                    .subtract(spring_min_agrimask)
                    .divide(vci_bot_spr)
                    .set('year', y);
}));
print("spring_vci",spring_vci);

//get means as table and add to Drive
var spring_vci_means_raions = spring_vci.map(function(image){
  return image.reduceRegions({
    collection:spring_vci,
    reducer:ee.Reducer.mean(),
    scale: 30
  });
});
var ukr_spring_table = spring_vci_means_ukr.flatten();
print(ukr_spring_table);
Export.table.toDrive(ukr_spring_table,"ukr_spring_vci");


texts: How to use default identify tool in ArcMap from C#?

texts:
I made a custom "Find" tool (winform) which is called after specific button pressed in an ArcMap toolbar. In this new "Find" form I need to use identify for one or more selected items in DataGridRow. So far I failed to find a way to call default "Identify" tool from ArcMap. I could make my own tool, but that seams like a waste of time, as I don't need any custom functionality in it.

Is there a way to call ArcMap "Identify" form (with feature properties etc.) from my custom winform tool passing one or more features?

I'm relatevly new to ArcMap as well as C#.


texts:
I am looking to perform raster math on two large rasters (slope ~1.3GB and aspect ~2.4GB). When I use the GIS raster calculator code below, QGIS errors out (out of memory), which is odd considering the computer has 20 cores and 32GB of RAM. I am hoping I could do this using gdal_calc, but am unsure if gdal_calc could perform such an equation, as well as the syntax difference.

("slope@1" <=5)*16 + (("slope@1" > 5) AND ("slope@1" <=38) AND ("aspect@1" >=168) AND ("aspect@1" <=192))*4 + (("slope@1" > 5) AND ("slope@1" <=38) AND ("aspect@1" >=90) AND ("aspect@1" <168))*1 + (("slope@1" > 5) AND ("slope@1" <=38) AND ("aspect@1" >192) AND ("aspect@1" <=270))*1 + ("slope@1" >38)*0 + ("aspect@1" <90)*0 + ("aspect@1" >270)*0


texts: Filter features based on attribute percentage threshold

texts:
I have a layer (vector layer imported from a Shapefile) in QGIS 2.8.3 with a bunch of features with a single numerical attribute, 'level'. I would like to filter the features with a threshold that is a certain percentage of the maximum value. For instance, if the threshold is 20% I would like to perform an expression like:

"level" > .20 * max("level")

or more SQL like:

"level" > .20 * (select max("level") from table)

These obviously doesn't work in the query builder or expression based filter. Is there a way to do this, with a different syntax, or do I need to use a different tool?


texts: How to custom the LayerSwitcher in OpenLayers to toggle all layers?

texts:
I want to have a checkbox in my layer switcher "turn all layers on/off".
how can I do that ?

I thought about adding an empty layer with that name and somehow bind her an event "On Show" and "On Hide" to toggle the rest of the layers.

Any better solution ?


texts: How to preserve ordering of results with the wrapper functions

texts:
When running

SELECT * FROM dijkstra_sp('ways',232, 421);

I get the result ordered by gid, rather than by the order for reaching the target by starting at the source.

What is the reason for that and is there a way to get the results in the natural routing order?

Im using pgrouting 1.05-1~lucid1 from http://ppa.launchpad.net/georepublic/pgrouting


texts:
In QGIS, when I insert a kml contour that was created in Google Earth, it appears misaligned in respect to all the other layers.

Note in this image how the red contour is moved upwards with respect to where it should be (use the house and the road going south as reference).

I am using EPSG:3003 as project CRS, all the different layers I have in my project from WMS services overlay correctly. The kml layer is set to EPSG:4326 (WGS 84) and I have set "automatic on the fly reprojection" in the options of QGIS. The CRS of the raster is set to EPSG:3003.

I have QGIS 2.0.1 but the same problem happens in 2.4.


texts:
I have a shapefile with 1722 rows representing a road network. There is a feature somewhere, somehow that is kind of invisible, specifically the feature with id 1780. When I execute an expression

 "id" = 1780

the bar shows there is a result somewhere. When I try to zoom on the selected result it's nowhere. I'm asking because I try to load the dataset into postgis and this way I'm getting a null geometry.

Any idea what is causing it?


texts: How can someone make this vertical section of a viewshed analysis (e.g. in ArcMap)?


texts: Calculate the distance between points and the start of each lines

texts:
I am trying to calculate the distance between the green points and the red point on their left and add that distance to the green points attributes.
The lines are roads and I made the red dots using the QChainage extension. I then added the roads ID to the red dots. After that I tried to calculate the distance using the Distance Matrix but the results weren't the one I expected.
Is there another way I could get the distance between the points and the point to their left?


texts: Navigator for ArcGIS integration with a laptop via bluetooth

texts:
Does anyone know how to integrate Navigator for ArcGIS with a laptop via bluetooth to display location information/point locations from a web app residing on a laptop?

Or can I just call the Navigator service via the bluetooth connection?


texts:
How can I display my current position in the QGIS canvas, similar to the "blue point" in google mapsa?

I need this to navigate to defined points and I do not want to change the program all the time.


texts: Exporting NDVI CSV data of a rectangular region

texts:
I am using the MODIS MOD13Q1 NDVI in the Google Earth Engine for my thesis project. I created a rectangular region to export CSV data. But I am having trouble exporting the NDVI data associated with the timespan that I am looking at. This is my code so far (the only thing that's not being exported is the NDVI data):
var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
                  .filter(ee.Filter.date('2000-03-01', '2020-06-01'));
var ndvi = dataset.select('NDVI');
var ndviVis = {
  min: 0.0,
  max: 8000.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};
Map.setCenter(30.12884910803246,7.396618943073097);
Map.addLayer(ndvi, ndviVis, 'NDVI');

var rectangle = ee.Geometry.Rectangle([-31.27142723303246, -30.12884910803246,31.27142723303246, 30.12884910803246]);

// Make a collection of points.
var features = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(31.27142723303246,7.396618943073097), {name: 'Ironman'}),
  ee.Feature(ee.Geometry.Point(31.27142723303246,9.591374979848126), {name: 'Cpt. America'}),
  ee.Feature(ee.Geometry.Point(30.12884910803246,9.591374979848126), {name: 'Thor'}),
  ee.Feature(ee.Geometry.Point(30.12884910803246,7.396618943073097), {name: 'Thanos'}),
]);

// Export the FeatureCollection to a CSV file.
Export.table.toDrive({
  collection: ndvi,
  description:'dataToDriveExample',
  fileFormat: 'CSV'
});


texts: Transformation from EPSG:31983 to EPSG:4326 change only one axis (R + ggmap)

texts:
I'm trying to convert my shapefile with spTransform(), but it's not successful. I don't know why it's changing only one axis. I want to change the axis into decimal degrees.

# Libraries ---------------------------------------------------------------
library(rgdal)
library(rgeos)
library(ggplot2)
library(ggthemes)
library(ggsn)
library(maps)
library(maptools)
library(grid)
library(plyr)

URL <- "https://osm2.cartodb.com/api/v2/sql?filename=public.galapagos_islands&q=select+*+from+public.galapagos_islands&format=geojson&bounds=&api_key="
fil <- "gal.json"
if (!file.exists(fil)) download.file(URL, fil)
gal <- readOGR(fil, "OGRGeoJSON")
topo <- readOGR("public.geojson.json", layer ="OGRGeoJSON")

so here if you change CRS("+init=epsg:31983"), it's working, but I don't know why it's not working with only 4326...

gal <- gSimplify(gUnaryUnion(spTransform(gal, CRS("+init=epsg:4326")), id=NULL), tol=0.001)
topo <- spTransform(topo, CRS("+init=epsg:4326")) # epsg:4326 is decimal degree, epsg:31983 is meters only

# double check that they match
identical(proj4string(topo),proj4string(gal))    ## [1] TRUE
gal_map <- fortify(gal)
topo_map <- fortify(topo)
topo_map <- join(topo_map, topo@data, by="id")
names(topo_map)

# Complete Galápagos Island!!!! -------------------------------------------
gg <- ggplot()
gg <- gg + geom_map(map=gal_map,
                    data=gal_map,
                    aes(x=long,
                        y=lat,
                        map_id=id),
                    color="black",
                    fill="#8FBC8F", #comment# #FFFFFF is white, # the colour #8FBC8F is a standard green for map's background
                    alpha = 0.4,
                    size= .1) # was .5
gg <- gg + coord_equal()
gg <- gg + ggsn:::scalebar(gal_map, # scale bar
                           dist=50,
                           location="bottomleft",
                           st.size=5)
gg <- gg + #theme_map() + # veryyyyyyyyyyyyyyyyyy sober
  theme(legend.position = "right",
        legend.text = element_text(size = 12,
                                   colour = "black",
                                   angle = 0),
        axis.text.x = element_text(size=12),
        axis.ticks.length = unit(.2, "cm"),
        panel.grid.major = element_line(colour = "blue",linetype = "dotted",size = .1),
        panel.background = element_rect(fill = "white"), # lightsteelblue2 for water # remove the default grey, you can add blue color for water (# background making illusion that this is water)
        panel.grid.major = element_blank(), # remove grid
        panel.grid.minor = element_blank(), # remove grid
        strip.background = element_blank(), # remove grid
        panel.border = element_rect(colour = "black", # border to the map
                                    fill=NA,
                                    size=1))
gg <- gg + geom_path(data = topo_map, # adding the elevation line
                     aes(long, lat, group=group, colour = elev),
                     size = 0.1 ,
                     alpha=0.4,
                     linetype = 1) +
#         geom_polygon(data = topo_map, # adding the elevation polygon
#                      alpha=0.01,
#                      aes(long,
#                          lat,
#                          fill=elev,
#                          group=group)) +
  coord_equal() +
  scale_color_gradient(low = "aquamarine", # guide = "none" if you don't want
                       high = "darkgreen",
                       guide = "colourbar") +
  labs(title = "Galápagos Islands",
       y = "lat",
       x = "long",
       colour = "Elevation")
gg

In the end, I want that the axis are shown in decimal degrees. Not in meters (31983). I don't want to do this manually (by writing the axis myself), I want to transform the scale.


texts: GeoServer CSS categorized marks with different size

texts:
I am doing categorized marks with GeoServer. The color is fine. I use same template of this question: geoserver css categorized marks  but I can't figure out how to change the size of my points.

This is what I've tried but my points displaying with the same size

* {
  mark: symbol(circle);
}

:mark [total < 500000]{
  mark-size: 2px;
  fill: yellow;
  stroke: yellow;

}
:mark [total >= 500000] [total < 1000000]{
  mark-size: 5px;
  fill: blue;
  stroke: blue;
}


texts: Why my Feature Class created with the tool "create feature from XY table" is outside the Feature Dataset?

texts:
I'm using ArcGIS 10.2.2 and I'm trying to create a Feature Class inside a Feature Dataset from a table opened through OLE-DB connection using the tool "Create Feature Class from XY Table".
The output Feature Dataset Coordinate System is the same than the Feature Class I want to create (WGS 1984), and , when I run the tool, I don't get any error message.
But the output feature class is created OUTSIDE the Feature Dataset!!
(Still inside the File Geodatabase, but outside the desired Feature Dataset).
Does anyone have an idea why this is happening and how I could avoid this problem ?


texts: Identify if the point lies inside the polygon in QGIS

texts:
I created a polygon using Iso-Area as polygons algorithm in QGIS. I have a point's coordinates and added it as a vector layer into QGIS. It is shown on the map that the point lies within the polygon. Then I wanted to create a dummy variable that indicates that the point lies inside the polygon. Using what I found in the following link I created a Virtual Layer by applying this query ('A' is the point layer and 'B' is the polygon layer):
SELECT A.*, st_within(A.geometry, B.geometry) AS dummy
FROM A, B

But the value of the created dummy variable is 0, while the point is inside the polygon. What is the solution?
Generating a dummy variable if point lies inside a polygon using QGIS?


texts: Permanent attribute join (for dissolving)

texts:
I am quite new with QGIS and I have run into a problem that I would love your suggestions to. Here we go:

I want to end up with a municipality map with the sum of perish population.

Here's what I have:

  Layer1: (Shape) Parishes with municipality code and perish code as
  attributes

  Layer2: (Attribute) population number with perish code as attributes

I am joining the layers, no problems so far, but when I want to dissolve and sum up population through the "Dissolve with stats" plugin, I don't see my joined attributes. I only see the original attributes and when I open the .dbf file in excel, I still only have the originals.

I figure that the join is not a "permanent" join, but just a "relational" one. Is there a way to change this or is there another way around it?


texts:
Below is the query I'm using :

SELECT osm_id FROM planet_osm_line where ST_DWITHIN((ST_Transform(way,4326)), st_setsrid(ST_Point(9.17193, 48.76768),4326) ,50000)
and time_cost =0.1 limit 100000;

I did Explain, I see that the index on geom column way is not being used.

I've two indexes, on time_cost and way

Below is the trace on Explain:

  Limit  (cost=0.43..245639.85 rows=2975 width=8) (actual time=0.510..1360.168 rows=100000 loops=1)
   ->  Index Scan using time_cost_idx on planet_osm_line  (cost=0.43..245639.85 rows=2975 width=8) (actual time=0.509..1322.503 rows=100000 loops=1)
         Index Cond: (time_cost = '0.1'::double precision)
         Filter: ((st_transform(way, 4326) && '0103000020E610000001000000050000001AA88C7FDA68E8C0475A2A6FE763E8C01AA88C7FDA68E8C0B9A5D5901870E840E6577380256BE840B9A5D5901870E840E6577380256B
E840475A2A6FE763E8C01AA88C7FDA68E8C0475A2A6FE763E8C0'::geometry) AND ('0101000020E61000000E677E35075822408AE5965643624840'::geometry && st_expand(st_transform(way, 4326), '50000'::double pr
ecision)) AND _st_dwithin(st_transform(way, 4326), '0101000020E61000000E677E35075822408AE5965643624840'::geometry, '50000'::double precision))
 Planning time: 0.284 ms
 Execution time: 1377.343 ms
(6 rows)

Hers is the updated trace after adding index as suggested:

 Limit  (cost=33699.36..183915.21 rows=2975 width=8) (actual time=639.022..2008.438 rows=100000 loops=1)
   ->  Bitmap Heap Scan on planet_osm_line  (cost=33699.36..183915.21 rows=2975 width=8) (actual time=639.020..1971.429 rows=100000 loops=1)
         Recheck Cond: ((time_cost = '0.1'::double precision) AND (st_transform(way, 4326) && '0103000020E610000001000000050000001AA88C7FDA68E8C0475A2A6FE763E8C01AA88C7FDA68E8C0B9A5D59018
70E840E6577380256BE840B9A5D5901870E840E6577380256BE840475A2A6FE763E8C01AA88C7FDA68E8C0475A2A6FE763E8C0'::geometry))
         Rows Removed by Index Recheck: 107717
         Filter: (('0101000020E61000000E677E35075822408AE5965643624840'::geometry && st_expand(st_transform(way, 4326), '50000'::double precision)) AND _st_dwithin(st_transform(way, 4326),
'0101000020E61000000E677E35075822408AE5965643624840'::geometry, '50000'::double precision))
         Heap Blocks: exact=15297 lossy=32836
         ->  BitmapAnd  (cost=33699.36..33699.36 rows=44623 width=0) (actual time=629.404..629.404 rows=0 loops=1)
               ->  Bitmap Index Scan on time_cost_idx  (cost=0.00..4169.81 rows=223117 width=0) (actual time=40.065..40.065 rows=252522 loops=1)
                     Index Cond: (time_cost = '0.1'::double precision)
               ->  Bitmap Index Scan on geom_way_line_idx  (cost=0.00..29527.82 rows=651121 width=0) (actual time=580.924..580.924 rows=3255604 loops=1)
                     Index Cond: (st_transform(way, 4326) && '0103000020E610000001000000050000001AA88C7FDA68E8C0475A2A6FE763E8C01AA88C7FDA68E8C0B9A5D5901870E840E6577380256BE840B9A5D5901870E
840E6577380256BE840475A2A6FE763E8C01AA88C7FDA68E8C0475A2A6FE763E8C0'::geometry)
 Planning time: 0.230 ms
 Execution time: 2025.651 ms
(13 rows)


texts: ArcGIS 10.1 - identify Vehicle Routes from Segment Counts

texts:
We have vehicle counts on every road segment in a small portion of our city, but we are trying to determine the route of individual vehicles on these segments. Is there an algorithm in ArcGIS 10.1 that can let us separate out the individual routes from our segment counts?
For example, if we know the counts in the attached diagram, we are trying to generate the two unique routes from that information.

This is like a reverse traffic-assignment method, from transport modeling.

The volumes are not especially high and segments are small, so we could probably brute-force identify trips, but we'd like to see if anyone has a better suggestion first.


texts:
Is there a way to batch read in ENVI standard image files (they have associated headers) and export them to a GeoTiff ?

I have a folder of ~500 images


texts:
My bigger problem is to convert a raster map to polygons. I've found some questions that give info and I'm trying to use GRASS GIS as part of that. However, I'm having some basic problems and I can't find the solution in the manual or tutorials. If someone can point me somewhere helpful that would be great.

I'm using r.in.gdal to import a tiff file (following this post How to convert image of map into vector format?) of this image

http://images.replacements.com/images/images5/crystal/D/dartington_chalice_clear_water_goblet_P0000019245S0003T2.jpg

It imports the tiff (which I created in Gimp) as red, green and blue bands. Now I'm trying to display them and thin them and here's where the trouble comes in.

When displaying...
Sometimes I get the error that the raster map  can't be found.
Sometimes I get a red box.
One time I did get the goblet, but not sure why!

When I thin, I'm not sure whether to thin the 3 bands separately or whether to composite them into one map and thin that. Neither really works.

As you can see I'm a newbie and I appreciate greatly any help. I've been able to do a fair job of this with WinTopo, but am trying to see if GRASS can do a better job.


texts: QGIS 3.0.0. Coordinate display incorrect

texts:
My CRS is reported by QGIS (project properties) as being ED50/UTM31N (EPSG 23031) and that seems to be correct as all layers are positioned correctly. However, the coordinate display box shows the latitude and longitude (in degrees) as being the same. When the coordinate display is set to "map units(meters)" it displays correctly.

Is this to do with the fact that on-the-fly transformations no longer appear to be available in 3.0.0?


texts: Is there a way to batch resample and mosaic/merge many large raster files at once?

texts:
I have over 400 large raster DEM files at 1 m resolution that I need to resample to 10 m resolution and merge/mosaic them all into one file. I am wondering if there is a way to do this all in one step or as few steps as possible. I assume it would be faster to resample first, then mosaic.
The reason I need to resample is because the other data I am working with is 10 m resolution. I have looked into the problem online, but have not found anything helpful yet.
Does anyone know any tools that can to do this? And which program can do it best? I have ArcGIS and QGIS but may be able to get access to other programs if needed.


texts: QGIS 3.4: error when cutting raster ("Cannot compute bounding box of cutline")

texts:
QGIS 3.4

During the process of trimming raster layers, the command failed.
Looking at the Journal, the root error seems to be :

  ERROR 1: Cannot compute bounding box of cutline. Cannot find source SRS

However, I have defined an CRS for each layout, and it is also shown in the windows selection.


texts: Is it possible to generate a heatmap on the server-side using ArcGIS Server?

texts:
Using ArcGIS Server's JS API it is possible to generate a heatmap at runtime, as in the example at http://tmappsevents.esri.com/website/heatmap/index.html

This requires having client-side access to the individual points in the heatmap, for example by downloading a feature layer to the browser.

The problem is that the browser becomes unusably slow when working with more than a few thousand points.

Has anyone been able to generate the heatmap on the server-side, and send it down as a cached or dynamic layer?

I need the ability to generate the heatmap at runtime, so that the points could be filtered with the change immediately visible in the heatmap.


texts: How to load big geo-data file into web application

texts:
I am building a web app where I show details in map. I load data using geojson file (190 mb).

I was looking in different possibilities to show the data, here there are:

Leaflet GeoJSon: it should make it easy to load huge data, but I am not seeing effects on my work.
Tile Layer http://bl.ocks.org/glenrobertson/6203331
It doesn't support some data viz libs that I am using
Leaflet Vector layer  http://geojason.info/2012/leaflet-vector-layers/
didn't get good results (having probs with data viz)
TopoJson : resized the file from 190mb to 170mb. I still face the same problem.

The simplest way to load the file is using jquery "getJson".

jQuery.getJSON(url, function(data){
..
}

Someone suggested to stream the JSON data rather than bulk-load it. I have tried https://github.com/tmcw/geojson-stream .
I couldn't get great results in all the cases.

Can you please suggest any better approach ?


texts: Finding the major axis of a polygon (the outline of a sinkhole or doline)

texts:
I am using QGIS to do a morphometric analysis of dolines (sinkholes). The dolines have different shapes from circular to elliptical and I need to calculate the major axis of the polygon that represents them. I already use the minimum delimiter box algorithm, which gives me information about the length and width. Now I need to represent the maximum axis inside the polygon and I can not find any method.

Images of minimum bounding box examples taken from Kobal 2015 and Bondensan et al, respectively.


texts: Reconcile Versioned Geodatabase with only DEFAULT version?

texts:
Follow up to this: Spatial Database Connection String

I am trying to reconcile a versioned geodatabase that has only one version, DEFAULT.  I am following the steps from the ESRI documentation: A quick tour of reconciling a version.  It seems like there is no reason to have to do this, but the "OLE DB" table and "spatial" table are out of sync until I do this.  As per comment from @mwalker in linked question above:

  If you're using versioned SDE then you may not see changes right away, since versions are handled with delta tables and views. The source tables won't get updated until the version is reconciled.


texts: Removing overlap within multipart polygon in ArcGIS for Desktop?

texts:
ArcMap 10.1 SP 1
License: Advanced

I am attempting to remove overlapping portions of a multipart without going through the steps exporting or creating a layer from each part and clipping individually. I've looked around at various answers, it seems the answer may reside in topology rules. I admit I'm not very familiar with using topology, so I was unsuccessful. In fact it appears my rule was simply not applied. No doubt due to my incompetence.

Is it possible to accomplish some fashion of "clipping" of a multipart polygon with its own parts?

Here is what I'm doing and what I've done; I will wager someone else has done a similar task. The data is wildfire data taken from modis, essentially I am building a progression map using modis data. The data comes as points, I used minimum bounding geometry to convert the data into a multipart polygon using the Julian date from the point data. This works as I don't need dead accuracy for my purposes, it's more of a quantity over quality project. Because the process of min bounding is a bit coarse, as is the point data, I cannot simply subtract the total acres from the previous day to get acres per day. That is the data I am actually seeking. After the minimum bounding I clipped the polygon to the final fire perimeter as the process created some polygon area outside the fire perimeter. This is really where I am stuck.

I've tried some merging and dissolving based on similar questions. Perhaps I'm missing the vision in those as they have yielded nothing useful. I also tried the topology route by creating a gdb. I created topology with a rule of must not overlap for the feature data set, then exported the data from my map into the gdb as a feature class. This didn't appear to do anything, so I'm rather certain I did not do this correctly. I then fumbled around with the topology tool bar with no success. I could do a union and then pick out the issues; this is tedious. I will likely need do this process again in the future, so avoiding tedium is a plus. Ideally I could make a model with the two inputs of modis data and fire perimeter to get at the end result. My model building is another conversation entirely.

Short version: The picture below may help me explain things. The polygon is multipart, Julian day 223 is sitting atop day 224. I would like 224 to not be under 223, essentially it would be clipped by 223. This would carry on through each day.


texts:
Building from my previous questions:

QGIS Server limitations?

I want to evaluate performance of both MapServer and GeoServer.  Can both servers exist on the same machine (each serving out WMS) or should I set them up on separate servers?


texts: Where to report QGIS and QGIS-Plugin bugs?

texts:
I often find bugs in QGIS and its Plugins.

I would like to ask about this on the internet but I do not seem to find the proper place to address this issue.

Are there any instructions?


texts: For my windows application, I am looking for a method to edit vertexes and geometry of features that are stored in a SQL Server database. I just want a simple function for dividing polygons into 3 parts and so on. Are there any open source library for these functions?