texts: Splitting line at local maxima and minima of raster?

texts:
I'd like to find those points along a polyline where a raster has local minima or maxima (find peaks and valleys where the line slope is zero).

Using ArcMap's slope tool didn't work as desired because it calculates the slope of the raster itself, not the slope of the line. This is a problem when the line goes over a hill, but doesn't touch the top: then the raster slope is non-zero, while the line slope would be zero.

How can I find peaks and valleys along a line?


texts: Extra space in labels on some pages of a map book resulting in invalid searches

texts: I have exported a map book (Arcmap 10.1). All the labels appear fine but on a particular page - if I search for label for a particular layer - I don't receive any result. I copy and pasted this label on notepad and noticed space being added to these labels e.g. this 9239 exists but Adobe Reader page is generating the error and when I copy this label and print to notepad, it is actually being read as '92 39'.


texts: QGIS 2.6.0 Atlas Export to multiple output files fails; RESOLVED

texts:
QGIS version 2.6.0, date 01/05/2014
OS Windows 7 Pro 64-bit

STEP 2 listed below notes resolution - Output should use field name(s) only (not 'attribute') to generate an atlas with multiple output files.

Description

When exporting Atlas features with Output to multiple files (i.e., 'Single file export when possible' checked OFF) and using either the Atlas export menu or the Atlas Preview export menu, Select Folder, the export then creates a single page file of the current atlas feature with the name of the Selected Folder and repeatedly overwrites this file as
each atlas feature is exported. Reproduced for both pdf and images (jpg, png).

Steps to reproduce:

Create a set of features for Atlas generation
Export file name set to be based on one of the attribute fields - THIS WAS THE ISSUE. Simply using the field name(s) for Output file name generates the desired result.
Select Folder
Export

Expected behavior:
This feature worked fine for me in mid-November 2014, creating a folder of multiple output files of one page for each feature. i can report a bug if it appears to be one, but because it worked previously wanted to see if anyone had experienced this and resolved it.


texts: Adding weights to shapefile for points in ArcGIS Pro using ArcPy?

texts:
I have a .csv file with longitude, latitude representing a point on map. And another column that I will use as weights for these maps. I want to calculate weighted mean center for all the points using ArcPy.

This is what I've tried so far:

import pandas as pd
import arcpy

df = pd.read_csv("data.csv")
longs = df['Longitude'].tolist()
lats = df['Latitude'].tolist()
wts = df['Weights'].tolist()
point = arcpy.Point()
pointGeoms = []

for idx in range(len(longs)):
    point.X = float(lats[idx])
    point.Y = float(longs[idx])
    point.M = int(wts[idx])
    pointGeoms.append(arcpy.PointGeometry(point))
arcpy.CopyFeatures_management (pointGeoms,"centroids.shp")

Above code is for generating shapefile. Below code is for calculating weighted means.

import arcpy

workspace = r"C:\Users\sam"
input_FC = "centroids.shp"
weight_field = "M"

try:
    # Set the workspace to avoid having to type out full path names
    arcpy.env.workspace = workspace
    # Process: Mean Center...
    arcpy.MeanCenter_stats(input_FC, MEAN_output, "#", "#", "#")

except:
    # If an error occurred when running the tool, print out the error message.
    print(arcpy.GetMessages())

I am not able to get any positive outcomes from both the codes. The shapefile generated by first code shows empty point values in ArcGIS Pro. And the second code generates a point randomly placed very far away from any point of interest.


texts:
Is TileMill as a service restricted to PNG? Or is JPEG an option?

I setup TIleMill as a service as documented here:

https://www.mapbox.com/tilemill/docs/guides/ubuntu-service/

Tiles are rendering locally e.g:

http://127.0.0.1:20008/tile/q/10/937/596.png

I've tried changing the extension with no luck.

I realise I can export as mbtiles and use TileStache or similar but thought I would as the question.


texts: How to determine elevation of points in shapefile from second shapefile using ArcGIS for Desktop?

texts:
I am a bit new to ArcGIS and have run into a snag.

I have a series of points in a shapefile that only have x and y coordinates. Additionally I have a series of other points in a separate shapefile that have x, y, and z coordinates. I have added both to a map and have generated contours from the latter of the two files.

My question is: Is it possible to determine the elevations of the points in the first shapefile based on the second shapefile and contours and if so how?


texts: Selecting based on distance and change in elevation

texts: I am trying to find areas near a series of roads that are withing 500 ft of the road, and have a change in elevation of +50 ft. I have been provided a 2 ft contour of the area, and a centerline layer for the roads in question. My thought process was to intersect the roads and topo, to create points which I could then create the 500 ft buffer from. From there, I am not sure how to proceed and select from this the areas of elevation im looking for


texts: Create Fishnet with row iterator (in python) - Number of rows and columns not sufficient?

texts:
I have a shapefile that contains over 100 squares of different sizes and I want to divide each of the squares into quarter sections.

I've written a short python script to do a row iteration and for each square, create a fishnet with 2 rows and 2 columns based on some predefined condition. Since the squares have different sizes, I set the numbers of rows and columns to 2 and set the cell height and width to zero. (Find the script below)

However, I get this error when I run the script although it's only the width/height pair that I set to zero.
"ERROR 000124: These two pair parameters cannot be zero at the same time"

Any suggestions how I can work around this?

SCRIPT:

#Loop thru rows of input polygons
with arcpy.da.SearchCursor(outfc, ['SHAPE@', 'OID@', 'SumNet']) as rows:
    for row in rows:
        ext = row[0].extent
        st = '%f %f' % (ext.XMin, ext.YMin)
        orien = '%f %f' % (ext.XMin, ext.YMax)

        #Create fishnet
        out = os.path.join(gdb, 'fish_{0}'.format(row[1]))
        if row[2] > threshold:
            arcpy.CreateFishnet_management(out, st, orien, '0', '0', '2', '2', labels='NO_LABELS', geometry_type='POLYGON')
        else:
            #arcpy.Copy_management(row, out)
            arcpy.CreateFishnet_management(out, st, orien, '0', '0', '1', '1', labels='NO_LABELS', geometry_type='POLYGON')


texts: ST_GeneratePoints returns the wrong number of points

texts:
I'm using PostGIS to generate random points on some bikeway data sets in my region, and I've run into an issue. For a few of the sets, I find that when I run my query I often don't get the right number of points returned. For example, if I ask for 4 points, I'll receive 3 points. If I increase the radius within which I'm querying, there's a higher chance of this happening (versus a smaller radius, where the query works fine most of the time).

Here's an example query, with a "large" radius:

SELECT ST_AsText(st_generatepoints) FROM
(SELECT ST_GeneratePoints(geom, 4)
    FROM (
        SELECT ST_Intersection(
            (SELECT ST_Buffer(
                (
                    SELECT ST_Union(geom)
                    FROM bikeways
                    CROSS JOIN (SELECT ST_GeomFromText('SRID=26910;POINT(-122.98054922394947 49.22650820352046)') AS ref_geom) AS r
                    WHERE ST_DWithin(geom, ref_geom, 0.0462677490234375)
                ),
                0.000010, 'endcap=round join=round'
            )),
            (
                SELECT ST_Buffer(
                ST_GeomFromText('SRID=26910;POINT(-122.98054922394947 49.22650820352046)'),
                0.0462677490234375)
            )
        )
        AS geom
    ) AS s
) as p;

I assume the way that I'm creating the query is the likely culprit, but being new to PostGIS I'm not sure how I could rewrite it. What kinds of situations would result in ST_GeneratePoints returning fewer than the requested number of points?


texts:
I am using the Leaflet List Markers Plugin for showing all features currently visible in the map in a list. Adding a single L.geoJson layer is no problem and works fine.

Now I am using multiple layers stored in a L.layerGroup. In fact, it's still one dataset, but splitted by the filter-option into multiple layers, merged to one layergroup again.

By adding the layerGroup to ListMarkers, there is no output:

var markerlist = new L.Control.ListMarkers({layer: myLayergroup, itemIcon: "blue-marker.svg"});

Output list in my map:

Does anyone know, whether there is a way to insert a layerGroup to ListMarkers or a workaround?


texts:
Is there an editor for Leaflet JavaScript that has intellisense functionality?

I've tried both Visual Studio 2015 and Visual Studio Code ...

in 2015, if I type code like var geojsonlayer = L and hit the CTRL + SPACE combination for intellisense, I get offered  Location (in lib.d.ts) rather than L.
I get no intellisense in VSC


texts: Implementing signals on add feature with PyQGIS?

texts:
I'm trying to implement the code below to run when adding a feature in the layer, but nothing happens even an error message.

However when I call the method directly, the routing is executed normally.

Any suggestions for the problem?

from qgis.core import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.PyQt.QtWidgets import QWidget

def Form(dialog, layer, feature):
    layer.featureAdded.connect(lambda: fillFields(dialog, layer, layer.dataProvider(), 'Id', 'Tag', "CE-C"))

def fillFields(Dialog, Layer, Provider, FieldId, FieldTag, LabelTag):
    Layer.startEditing()
    lineEdit_tag = Dialog.findChild(QLineEdit, FieldTag)
    id = Provider.fieldNameIndex(FieldId)
    tag = Provider.fieldNameIndex(FieldTag)
    num = Provider.featureCount()
    lineEdit_tag.setText(LabelTag + str(num + 1))
    list = range(num)
    i = 0
    j = 1
    for item in list:
        new_values = {id: j, tag: LabelTag + str(j)}
        Provider.changeAttributeValues({i: new_values})
        i += 1
        j += 1
    Layer.commitChanges()


texts: Repairing broken GDAL and PROJ.4 on Ubuntu?

texts:
I solved the library problem, but others remain, only vaguely related to that of PROJ.4 library not loaded. I asked a follow-up question for these.

My GDAL package and PROJ.4 library are damaged, which causes errors while trying to reproject files (see old question below). I don't know whether this problem started while I mistakenly installed another instance of GDAL or before; however, I've uninstalled it and reinstalled QGIS and the problem remains.

After each reinstall, QGIS has the same plugins as before, so the part causing the errors is left untouched by the install/uninstall process.

What have I tried:

sudo apt-get remove build-essential python-all-dev and sudo apt-get remove gdal-bin to purge the manually installed GDAL
(un)install of QGIS with Ubuntu Software Center
sudo apt-get --purge autoremove
dpkg -l | grep '^rc' | awk '{print $2}' | sudo xargs dpkg --purge and dpkg --list | grep "^rc" | cut -d " " -f 3 | xargs sudo dpkg --purge to purge the configuration files

Still no luck. I've found something that seemed related, but that's for Windows, not for my Ubuntu 14.04. What else (and how) should I purge before next install of QGIS?

The original question:

I try to reproject some shapefiles with ogr2ogr. However, I always get the following error:

ERROR 6: Unable to load PROJ.4 library (libproj.so), creation of
OGRCoordinateTransformation failed.
Failed to create coordinate transformation between the
following coordinate systems.  This may be because they
are not transformable, or because projection services
(PROJ.4 DLL/.so) could not be loaded.

I tested this with various transformations, including such basic combinations as from WGS84 (EPSG:4326) to UTM (EPSG:3857), and this error occurs all the time.

The problem might be caused by the fact that although I have GDAL with my QGIS (2.10), I mistakenly installed it once more (from source, doing it as this answer suggested).

After facing this problem few times, I tried to uninstall the new GDAL instance (sudo apt-get remove build-essential python-all-dev and sudo apt-get remove gdal-bin - not sure whether it is cleanly removed or not).

In /USR/LIB/ I have libproj.so.0 and libproj.so.0.7.0, not libproj.so. I can't see it even with hidden files shown. However, when I try ln -s libproj.so.0.7.0 libproj.so (or the same with sudo), it complains that libproj.so already exists.


texts:
I have just started using QGIS 3.4 and want to split a polygon with a line from another layer. I can do this in 2.18 but there is no option in 3.4

Is this correct, is the tool not available


texts: Determining centroid of clusters of points in QGIS?

texts:
I've basically got the same issue as Identifying cluster centroids using ArcGIS Desktop?, but can't find a way to transplant the ArcGIS solution to QGIS.

I have household-level geo data for 10 villages spread out across the country. I need to find the central point of these household clusters, so that I can create a 10km buffer around that central point.

I am new to QGIS.


texts: Supervised classification using endmember spectra (using ENVI 5.1)

texts:
We are trying to map land cover classes over a small watershed. We have selected training sites (during a field campaign in early 2017) and determined their spectral profiles on a Landsat 8 image acquired at the time of field surveying. In order to assess the land cover changes, we wanted to map the same cover classes at a previous year. Since our training sites might not be relevant, we wanted to perform supervised classification using endmembers spectra instead of ROIs. When importing those spectra inside ENVI's Endmember Collection toolbox, it appears that only Spectral Angle Mapper and Spectral Information Divergence classifiers could be used. Common algorithms such as Maximum Likelihood or Mahalanobis distance fail, returning the following error :

  Problem: the selected algorithm requires that the collected endmember
  spectra all contain an associated covariance. ENVI is unable to
  continue because some of the endmembers collected to not have their
  covariance.

Could anyone help here ? Actually is that method relevant ? How can we possibly perform supervised classification using Maximum Likelihood/ Mahalanobis classifiers ?


texts: How to calculate the distance to multiple locations from each cell of my raster?

texts:
I'm working with Arcgis 10.1 and Python 2.7.2
How can I calculate the distance to multiple locations from each cell of my raster?

I tried euclidean distance but it only provides the distance to the
closest location. I need to know the second closest, third closest and so on until the 150th closest.

I tried generating a near table but I'm working with high resolution data. The lowest I could go for this dataset are 1,9 million cells. A near table would generare 150 times 1,9 million rows. Looping through these rows takes up a tremendous amount of time.
Here is an example (scratch of course) of what I'm trying to accomplish.

Red dots are locations, grid represents (49) cells, green lines are to show the distance.
I thought maybe something with an alternated euclidean distance function or an addition of euclidean distances to each single location. I tried combining euclidean rasters, but it didn't work.
Can anyone help me or point me in the right direction? A solution in python would be appreciated but I'm glad for any idea , workaround, anything...


texts:
I've been given the task of assigning address ranges to a centerline feature class.  I have all the data and I am just going street by street.  My question:  When assigning ranges, should I just enter the values I find for present addresses or would it be better to extrapolate the range over of the segment based on present values?  I am attempting to have further geocoding results be somewhere close to where they should be.

From the example in the graphic, where the x's are empty space and not unknown addresses:

Section A would clearly be from 1 to 3.
Section B: from 5 to 6 or from 4 to 6?
Section C: from 12 to 13 or from 7 to 13?

There won't be too many geocoding projects where the results will need to be spot on, but there are long stretches of road with scattered addresses (section C) and areas with clusters (section B)

Another question, should the ranges be continuous or will breaks in the range and null values (stretch or side of the road with no addresses) overly effect the geocoding process?

Thank you

For the need of real example of real data:

The read dots are on the street in blue and range from 3115 - 3141 from left to right.


texts:
I would like to travel to the loneliest place on earth. I define it as being the farthest from any human settlement. The place has to be on the land and not in the sea.

I tried to google, but couldn't find any information on how to calculate it.


texts: Google Earth Engine - counts days of the longest dry period in a time period

texts:
I'm trying to map the max number of consecutive days with rain <1 mm in Google Earth Engine.

This is the link to the full code in GEE, you can see and run it

https://code.earthengine.google.com/22b5c20d2700a2ffb5989f892838ac58

This is the portion of the code that calculates the longest dryperiod in a x time period

// dry period calculation
function drylength(current, previous){
  previous = ee.List(previous);  // cast to establish the client-side type
  var mask = current.remap([0], [1], 0,"precipitation").rename('drydays');
  var sum = mask.add(previous.get(-1)).multiply(mask);

  return previous.add(sum);
}

var list = ee.List(dataset.iterate(drylength, ee.List([
  dataset.max()
    .remap([0], [0], 0, "precipitation")
    .rename('precipitation')
])));
var dryperiod = list.get(-1);

When I run the code for counting the number of days of the longest dry period, the code returns the right number only if the dry period reach the end of the selected time period test.

For instance if I am looking for the longest dry period in 4 days timestep i get the following:

rain days 1 2 3 4
          0,0,1,1 = 0 dry days (expecting 2)
          0,1,0,0 = 2 dry days (OK)
legend:
0 = no rain days
1 = rainy days

How can I solve this?


texts: Arcpy update cursor: stop iteration error

texts:
I have 1 csv with municipalities and their corresponding zip codes and a GDB of houses which contain address, municipality but no zip codes

Objective: add the zip codes to the GDB table

import csv
import arcpy
import arcgisscripting
arcpy.env.overwriteOutput = True
gp = arcgisscripting.create()
file = "path\\munzip.csv"
houses = "path\\geocoding.gdb\\houses"
csvfile = gp.SearchCursor(file) #personally prefer using the object search cursor than arcpy.da
munzip = {}
for x in csvfile:
    mun = x.MUN.split()
    zip = "0"+str(x.ZIP)
    mun = str(mun[0])
    munzip[mun] = zip
with arcpy.da.UpdateCursor(houses, ("City","ZIP")) as cur:
    for row in cur:
        for k,v in munzip.items():
            if k == str(row[0]):
                row[1]=v
                cur.updateRow(row)


texts: Which layer creation option (lco) to pick in v.in.dxf?

texts:
I'm trying to properly upload a DXF via the GRASS option into QGIS. It is a building structure of LoD2 quality. Unfortunatelly, while the importer runs I always get the error:

WARNING: Vector map  is 3D.
Use format specific layer creation options (parameter 'lco') to export
<in 3D rather than 2D (default).

I couldn't really find a specification what lco one has to pick to properly display 3D structures. From the documentation I tried

lco="SHPT=POLYGONZ"
lco="SHPT=POINTZ"
lco="SHPT=ARCZ"

I'm using QGIS 3.16 on a MAC OS X.


texts:
When I start qgis with the terminal I get the following error:

QGIS died on signal 11[New LWP 15278] [New LWP 15279] [New LWP 15280] [New LWP 15283] [New LWP 15284] [New LWP 15286] [New LWP 15292] [Thread debugging using libthread_db enabled] Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1". 0x00007fa1f82b96c2 in __GI___waitpid (pid=15293, stat_loc=0x7fff74e3b6ac, options=0) at ../sysdeps/unix/sysv/linux/waitpid.c:30 30   ../sysdeps/unix/sysv/linux/waitpid.c: No such file or directory. [Current thread is 1 (Thread 0x7fa1fc1a5980 (LWP 15275))]
#0  0x00007fa1f82b96c2 in __GI___waitpid (pid=15293, stat_loc=0x7fff74e3b6ac, options=0) at ../sysdeps/unix/sysv/linux/waitpid.c:30
        resultvar = 18446744073709551104
        sc_cancel_oldtype = 0
#1  0x000055c5eb4217ce in  ()
#2  0x00007fa1f8213f20 in <signal handler called> () at /lib/x86_64-linux-gnu/libc.so.6
#3  0x00007fa1f835a5ab in __strncmp_sse42 () at ../sysdeps/x86_64/multiarch/strcmp-sse42.S:270
#4  0x00007fa1ec4b5d77 in pj_param_exists () at /usr/lib/x86_64-linux-gnu/libproj.so.13
#5  0x00007fa1ec4c99b4 in proj_create () at /usr/lib/x86_64-linux-gnu/libproj.so.13
#6  0x00007fa1f6dd8a67 in  () at /usr/lib/libgdal.so.20
#7  0x00007fa1f6dc8520 in OGRSpatialReference::importFromProj4(char const*) () at /usr/lib/libgdal.so.20
#8  0x00007fa1fa957244 in QgsCoordinateReferenceSystem::setProj4String(QString const&) () at /usr/lib/libqgis_core.so.3.4.3
#9  0x00007fa1fa95a2e5 in QgsCoordinateReferenceSystem::createFromProj4(QString const&) () at /usr/lib/libqgis_core.so.3.4.3
#10 0x00007fa1fa95b01c in QgsCoordinateReferenceSystem::fromProj4(QString const&) () at /usr/lib/libqgis_core.so.3.4.3
#11 0x00007fa1fa9a7f40 in QgsEllipsoidUtils::ellipsoidParameters(QString const&) () at /usr/lib/libqgis_core.so.3.4.3
#12 0x00007fa1faa6cfbb in QgsMapSettings::setEllipsoid(QString const&) () at /usr/lib/libqgis_core.so.3.4.3
#13 0x00007fa1f649071a in  () at /usr/lib/libqgis_gui.so.3.4.3
#14 0x00007fa1f8e206cf in QMetaObject::activate(QObject*, int, int, void**) () at /usr/lib/x86_64-linux-gnu/libQt5Core.so.5
#15 0x00007fa1faf44845 in QgsProject::ellipsoidChanged(QString const&) () at /usr/lib/libqgis_core.so.3.4.3
#16 0x00007fa1faaf48db in QgsProject::setEllipsoid(QString const&) () at /usr/lib/libqgis_core.so.3.4.3
#17 0x00007fa1fb601172 in QgisApp::fileNew(bool, bool) () at /usr/lib/libqgis_app.so.3.4.3
#18 0x00007fa1fb615e92 in QgisApp::QgisApp(QSplashScreen*, bool, bool, QString const&, QString const&, QWidget*, QFlags<Qt::WindowType>) () at /usr/lib/libqgis_app.so.3.4.3
#19 0x000055c5eb41d5d5 in  ()
#20 0x00007fa1f81f6b97 in __libc_start_main (main=0x55c5eb41b760, argc=1, argv=0x7fff74e3d998, init=<optimized out>, fini=<optimized out>, rtld_fini=<optimized out>, stack_end=0x7fff74e3d988) at ../csu/libc-start.c:310
        result = <optimized out>
        unwind_buf = {cancel_jmp_buf = {{jmp_buf = {0, -8979781846726237085, 94308543893936, 140735154477456, 0, 0, -2948634281101908893, -2906720311133100957}, mask_was_saved = 0}}, priv = {pad = {0x0, 0x0, 0x7fa1fc013733 <_dl_init+259>, 0x7fa1fbfd1538}, data = {prev = 0x0, cleanup = 0x0, canceltype =
-67029197}}}
        not_first_call = <optimized out>
#21 0x000055c5eb4211da in _start () gdb returned 0 Aborted`

No idea what is happening.

I had some problems in the installation, I had to create a fake gdal-abi-2-2-3 package and had to install grass-core_7.4.0-1_amd64 instead of grass-core_7.4.3 for dependency problems with Virtual package grass740.


texts:
I am trying to find out if the v.net.connect GRASS tool is on QGIS 3.4 because I'm not finding it anywhere in the toolbox.

I have to be reverting to QGIS 2.18 to run this tool.


texts: How do I separate .txt data into multiple columns when exporting to .CSV?

texts: I have .txt data that is described as Lat Lon Elevation (meters) and when I export to .CSV in Excel it puts all data into 1 column. So when I import the data into QGIS it all dumps into 1 column, and I cannot change the X and Y value column names when importing as Delimited Text File. Help!


texts: Plotting point based on angle and distance from another given point in QGIS?

texts:
I am using QGIS 2.18.10.

If I have a given point with an easting and a northing, and know its distance and angle from a transit which took the initial point, how do I plot the location of the transit if I don't already have its easting and northing?

As a secondary question; the transit took a further 900 points (an archaeological scatter which I'm plotting). Once I have the location of the transit, how would I calculate and plot those other points as a batch, based on their distances and angles from the transit, instead of inputting them individually?


texts:
I need to build an app to analyze some raster layers in QGIS.
I thought it was possible to just write code in python (2.7.6) (IDLE) and let it automatically open QGIS, load some layers, do a distance matrix, export that and load this file with python and use that information in the app.
Is something like this possible?
I have a mac running OSX10.9.1, using python 2.7.6 because I constantly got errors when using python 3.


texts: d3.js rectangle feature plotted as arc in Mollweide projection

texts:
I am trying to plot a rectangle using d3.js under the Mollweide projection. However, when I run the following script the rectangle appears as an arc or curved rectangle. I find the bahaviour somewhat strange since the parallels in Mollweide projection are just straight lines. Could someone explain the effect? The same effect persists in other latitudes and other projections (even the Equirectangular one) as well.

<!DOCTYPE html>

<head>
    <script src="https://d3js.org/d3.v4.min.js"></script>
    <script src="https://d3js.org/d3-geo-projection.v1.min.js"></script>
</head>

<body>

</body>

<script type="text/javascript">

var mainSVG = d3.select('body').append('svg').attr('width', 500).attr('height', 500);

var projection = d3.geoMollweide().scale(400).translate([250,750])
var path = d3.geoPath().projection(projection);

d3.json('test.json', function(error, vData) {

    var features = vData.features;

    mainSVG.selectAll('path')
        .data(features)
      .enter().append('path')
        .attr('d', path)
        .style('fill', 'red')
        .style('stroke', 'black');

})
</script>

the geojson file:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "DN": 0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 7.6, 65.0 ], [ 60.4, 65.0 ], [ 60.4, 63.4 ], [ 7.6, 63.4 ], [ 7.6, 65.0 ] ] ] } }
]
}

and the resulting figure:

An even more weird result when I try to plot several neighbouring rectangles at different latitudes:


texts: Seeking open source solution to access rasters in File Geodatabase?

texts:
I'm looking for a way to access raster layers in an Esri File-GeoDatabase without using ArcObjects or ArcPy, specifically I would like to access the raster layers with GDAL. I understand that the GDAL FileGDB and OpenFileGDB give access to vector datasets within File-GeoDatabases, and that ESRI's own FileGDB C++ API does not support raster access.

Is anyone out there aware of any open source solutions?


texts: Loaded OSM Data into PostGIS Database. Cannot view in ArcMap "The coordinate system identifier is invalid"

texts: I loaded OSM File for Pennsylvania into a PostGIS Database in PostgreSQL using osm2pgsql. Everything seemed to go without a hitch until I finally connected to the Database in arcCatalog and tried to view the layers. My plan was to use these for a basemap. When I tried to view the layers I was getting the following error and I'm not sure if there is a resolution to this, I have never seen it before.


texts:
I have layer containing a area that will be flooded if the water level rises 6 meters. I also have a layer containing roads for the whole area. I want to export all the roads that will be under water in case of a 6 meter water level rise.

What I have tried is to use Select by location with roads as target layer and the 6meter flooding layer as source layer. I have used ALL the spatial selection methods but none of them selects only roads that inside the flooding zone. I guess it's because it selects the roads, but the roads is composed by different line segments.

How should I do to select the roads correctly?

I am using ArcMap 10.2.2 with all the extensions.


texts: Counting upstream dams along river network in QGIS

texts:
I am working on QGIS. I have a river network where each segment is a line with ID from the upstream one (if any) and the downstream one (if any).
I also have information on:

location of dams
location of fisheries communities

I want to assess the impact of dams on quantities fished by each community.
Is it possible to count the number of upstream dams along the river for each community in QGIS?
Here is an example where dams are represented in red and fisheries communities in green. I would like, for each community, to know how many dams are located upstream along the river. For example, there is 0 upstream dam for communities A. For point C, there is one dam, but no dam for point B since it is located along a tributary of the main river and there is no dam along this tributary.The same applies for point E. Finally, for points D, there are 2 upstream dams.


texts: Is preproccessing a Landsat scene necessary for image classification?

texts:
I am classifying Landsat TM scenes for the entire Chesapeake Bay as part of my Master's research. Scenes have been collected from 1986 to 2016 in 5 year increments for my area of interest. I recently came across various papers that mentioned how to preprocess Landsat scenes for various analyses but none of them mentioned image classification (be it unsupervised or supervised). Is it necessary to preprocess these scenes before classifying them? As stated in my previous question (Separating spectral similar training sets to improve accuracy of classified Landsat scenes in ArcMap?), I have had issues with distinguishing beach and urban pixels in my supervised classifications. The mosaic of Landsat scenes I am working with is below:

The color distortion in the right scenes concerns me. Will this negatively affect my classifications? How could I go about fixing this? Will preprocessing help with this at all?

If so, are there specific tools or preprocessing procedures I should be using for my scenes?


texts:
I have a WMS that I would like to consume with ArcGIS for Server.  The WMS is secured by way of an IP whitelist, and network connectivity from the server to the WMS has been verified.

Is there a way to consume this WMS and then re-publish as a service to the ArcGIS API for Javascript?

Calling the WMS directly from the Javascript code is not an option as each client will have a different IP address, and this will be unknown until run-time.  I therefore won't be able to pre-populate the IP whitelist with a list of all known IP addresses.

I know I can do this with MapServer, but hoping not to have to move down the path of having multiple software components installed on the infrastructure.


texts: Create Strictly Simple Polygon in PostGIS/Shapely

texts:
The clipper lib has an amazing function to create a strictly simple polygon (definition below).

Is it possible to do this in simple features (pref postgis or shapely)?

A strictly simple polygon is a simple polygon that does not contain 'touching' vertices, or 'touching' edges.

http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Classes/Clipper/Properties/StrictlySimple.htm


texts: Updating firstPoint and lastPoint using ArcPy cursor?

texts:
I want to update a Line feature's start and end.
I got Point of first and last.
BUt, I can't update them.

What is wrong ?
(I'm using ArcGIS 10.2)

   rows = arcpy.UpdateCursor(srcFile)
   feat = row.getValue(shapeName)

   bufFeat1=feat.firstPoint.X
   bufFeat2=feat.firstPoint.Y
   bufFeat3=feat.lastPoint.X
   bufFeat4=feat.lastPoint.Y

   row.firstPoint.X=bufFeat3
   row.firstPoint.Y=bufFeat4
   row.lastPoint.X=bufFeat1
   row.lastPointY=bufFeat2

   rows.update(row)


texts: Python script for raster calculation using GDAL

texts:
I am trying to loop over a raster calculation using GDAL. However, the loop appears not to work and QGIS crashes (without the loop the code which I found somewhere here works). The problem appears to be in the lines where I define input or output file. Any idea what could be the problem?

s = ['F101992','F101993','F101994','F121994','F121995','F121996','F121997','F121998','F121999','F141997','F141998','F141999','F142000','F142001','F142002','F142003','F152000','F152001','F152002','F152003','F152004','F152005','F152006','F152007','F162004','F162005','F162006','F162007','F162008','F162009','F182010','F182011','F182012']

for i in range(len(s)):

nighttime = "{mypath}//"+s(i)+".tif"

GDALnighttime = gdal.Open(nighttime, GA_ReadOnly )
band1 = GDALnighttime.GetRasterBand(1)
Rdata = BandReadAsArray(band1)

dataOut = Rdata-((Rdata>63)*(Rdata-63))-((Rdata<=5)*Rdata)

outFile = "{mypath}//"+s(i)+"_ic_1.tif"

driver = gdal.GetDriverByName("GTiff")
dsOut = driver.Create(outFile, GDALnighttime.RasterXSize, GDALnighttime.RasterYSize, 1, gdal.GDT_Float32)
CopyDatasetInfo(GDALnighttime,dsOut)
bandOut=dsOut.GetRasterBand(1)
BandWriteArray(bandOut, dataOut)

dsOut =None
bandOut = None
driver = None


texts:
Does anybody know is CRS:84 in WMS v.1.3.0 getCapabalities response is same projection as WGS84 (EPSG:2180)(EPSG:4326)?

And is [westBoundLongitude],[eastBoundLongitude],[southBoundLatitude],[northBoundLatitude] in EX_GeographicBoundingBox respectively equivalent with [minx],[miny],[maxx],[maxy] in earlier version of WMS response (1.1.1)?


texts: Executing functions after multiple data sources have been loaded in Leaflet AJAX

texts:
I have a leaflet map that takes in data from multiple sources via the LeafLet Ajax plugin.

I want to execute certain functions only after ALL the data has been loaded, not just one layer.

The data loading part of my code is fairly simple and looks like this:

data_layer1=L.geoJSON.ajax(data1.geojson,{options});
data_layer2=L.geoJSON.ajax(data2.geojson,{options});
data_layer3=L.geoJSON.ajax(data3.geojson,{options});

Now if I want to execute a function after data_layer1 is loaded, I will have to

data_layer1.on("data:loaded",function(){
      //function body
});

But what I want to do is execute a function after data_layer1, data_layer2, and data_layer3 is loaded. How would I do that?

I am fairly new to AJAX, so the Asynchronous part of this is a bit disorientating.


texts: Zoom after apply CQL_FILTER in OpenLayers

texts:
I apply a filter to one layer and works nice. I want to maximize the bounds where I applied the filter, like this:

// mapbounds = array(minx,miny,maxx,maxy)

activeLayer = new OpenLayers.Layer.WMS(
      layername, layerurl,
      {
         width: map.size.w,
         srs: map.projection,
         layers: layername,
         height: map.size.h,
         styles: '',
         format: 'image/png',
         transparent: false
      }
      );

 // for example: conditions = 'igh > 5000 AND BBOX('the_geom,1,2,3,4)'
 activeLayer.mergeNewParams({
         'cql_filter': conditions
      });

bounds = new OpenLayers.Bounds(mapbounds[0],mapbounds[1],mapbounds[2],mapbounds[3]);

map.zoomToExtent(bounds);

The filter works fine, but the zoom is restored too quickly after apply it. It's applied, but the layer restore to the initial view instantly.

What am I doing wrong?


texts:
I am working with MODIS Terra (and simultaneous ASTER) Level 1B swath data using ENVI and am having one of these "this should be easy" questions about the most direct way to get to geo-referenced datasets of radiances (and reflectances). There are two ways I'm familiar in reading in this data:

If I download the HDF L1B swath datasets they come as files named MOD021KM.Ayyyyddd.[whatever].hdf containing scaled integers, with the scaling factors in the header. The geolocation information is in separate files MOD03.Ayyyydd.[whatever].hdf, which I have as well. ENVI is able to convert the scaled integers into radiances, but not to apply the geolocation (to effectively turn this into a Lat/Long projecton). So I'm currently stumped with this approach.
Alternatively, I can download my swaths band-by-band as GeoTIFF from the NASA LAADSWEB site, which combines the relevant MOD021KM and MOD03 files. I read them into ENVI and stack the bands I need. But then I still need to convert the data into radiances/reflectances, which requires me to download the HDF file as well and search around in the header for the correct coefficients.

Right now I'm doing 2. But is there a way of going down the path of 1? ENVI's "Geolocate MODIS" function stumps me -- it seems there's something for L3 data and L1A data, but nothing for L1B.


texts: How to install GDAL with TIFF support on Linux?

texts:
I am trying to install GDAL on a linux machine and tested some sample commands like gdalinfo --version This worked fine, but when I started working with images like .tif I got the following error .

Warning 1: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order
ERROR 1: MissingRequired:TIFF directory is missing required "ImageLength" field
gdalinfo failed - unable to open '1.tif'.

I reinstalled libtiff but I am still getting the same error.

Iam working on CentOS release 5.7 (Final) .Is there any one who worked on this.
Any help is appreciated.


texts: What is the best way to use Leaflet and Geoserver to highligh a feature after a WMS request?

texts:
I have a web map, and I'm hoping to make it interactive.

The layer I have is reasonably large, at 15 mb in GeoJSON format 4 mb compressed, so I don't want to load the whole layer as one file. But, I do want to be able to interrogate each feature for around 6 attributes and the geometry.

What I have is:

Leaflet 0.5
GeoServer 2.2
Access to the server and data, happy to customise the data.

What I want is:

user can click on a map, and is returned a nice dialog box with feature information in a table
the map centers on the feature, possibly zooms, and highlights the feature.

I'm pretty sure this can be achieved using a GeoJSON request (like the example here), but I can't find examples of this anywhere and I'm struggling to implement it myself.

Does someone have a nice example of this?


texts: Using updated wikimapia API to export data

texts:
My goal is to get a json export from Wikimapia. I'm trying to use their updated API but I always get no objects returned no matter what coordintes or categories I specify.

Here is my example query (?key=example is a valid key, but I also tried my own API key with the same result):

http://api.wikimapia.org/?key=example&function=place.getbyarea&coordsby=bbox&bbox=46.7674%2C56.1952%2C48.5939%2C56.8782&category=17&count=100&format=json

P.S. this answer is no longer valid as it uses deprecated methods


texts: Why does Field Calculator say failure during processing, check geoprocessing results window for details?

texts:
I am using the Field Calculator and keep receiving the following error: 'There was a failure during processing, check the geoprocessing results window for details'. I am not very tech-savvy but I have tried a few simple fixes and still have had no luck. It's only a simple division equation so i'm confused.

My equation is [Sum_NO_OF_] / [KM2AREA]


texts: Keep only the last marker when geolocating

texts:
I'm working with OpenLayers 2 and I managed to create a javascript to locate the users and it works fine. I added a marker representing the current location and when the position change, it creates a new marker. The problem is that it also keeps all the other markers. I'd like to find a way to keep only the last one.

Here is my code:

    function startWatch(){

      if (navigator.geolocation)
        watchId = navigator.geolocation.watchPosition(successCallback,
          errorCallback,
          {enableHighAccuracy:true,
            timeout:Infinity,
            maximumAge:0});
      else
        alert("Your browser can't get your location");
    }

    function stopWatch(){
      navigator.geolocation.clearWatch(watchId);
      Emplacement.clearMarkers();
      NearestHospital.setVisibility(false);
    }

    function successCallback(position){

      document.getElementById("lat").innerHTML = position.coords.latitude;
      document.getElementById("long").innerHTML = position.coords.longitude;

      latitude = position.coords.latitude;
      longitude = position.coords.longitude;
      var coord = new OpenLayers.LonLat(longitude,latitude).transform(new
      OpenLayers.Projection("EPSG:4326"),("EPSG:900913"));
      map.setCenter(coord, 16);

      Myposition = new OpenLayers.Layer.Markers("Myposition");
      map.addLayer(Myposition);
      var size = new OpenLayers.Size(15,15);
      var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);
      var icon = new OpenLayers.Icon('img/loc_marker.png', size, offset);
      var mylocation = new OpenLayers.Marker((coord),icon);
      Myposition.addMarker(mylocation)

    }

    function errorCallback(error){
      switch(error.code){
        case error.PERMISSION_DENIED:
          alert("User did not share geolocation data");
          break;
        case error.POSITION_UNAVAILABLE:
          alert("Could not detect current position");
          break;
      }
    }


texts:
I have been asked to add an interactive map to my company's annual report.

The requirements are that it should have multiple layers that can be turned on and off, zoom in/out and pan functionality, and clickable points on the map that will display images.

It must be entirely standalone as there's no guarantee that the reader will be connected to a network when they're viewing the document.

Any suggestions?


texts: Histogram chart of categorical attributes in QGIS

texts:
Using mmqgis I calculated the hub distance from source points to destinations but I would like to add histograms or similar diagrams to show the hub characteristics of the hub members (count of Child by type, count of transport by type). I've tried adding histogram and piecharts but I'm not sure if I'm starting from the correct point or whether its the configuration of my diagram which is at fault.

Layer Properties: Line (WKB type: "LineString")
Feature Count:100

Attribute Structure

RowID     HubName     Child     Transport
1         School1     Boy       Bike
2         School1     Girl      Car
3         School1     Boy       Car
4         School2     Boy       Walk
.         .           .         .
.         .           .         .
100       School6     Girl      Bike

Desired Output

Beside School1 I would have diagram showing the composition (CountOfBoy =x, CountOfGirl =y, Transport-Bike=z, Transport-Car=a). in chart form bdeside each Hub. I expect I'll be repeating the analysis with different destination points so I'd love to get an efficient way of accomplishing this without too much back and forth importing and exporting attributes and creating summary tables


texts: Convert values from square meter to square kilometer in ArcMap

texts:
I am using ArcMap 10.3.1 on a Windows x64 system.

I use the Universal Soil Loss Equation (USLE) to calculate the average annual rate of erosion on a field slope based on rainfall pattern, soil type, topography, vegetation system and management practices.

                         A= R*K*LS*C*P

A represents the average annual soil loss in tonnes per square kilometer per year (t km−2 yr−1);
R stands for rainfall erosivity factor;
K is the soil erodibility factor;
LS is the slope length-gradient factor
C is dimensionless vegetation cover factor;
P refers to the soil conservation practice.

In my case, all the data are 30m spatial resolution, and I calculate each parameter in ArcMap. So the spatial resolution of my soil erosion raster map is 30m, which means the unit is in tonnes per 900m2 per year.

I’d like to convert unit from ‘tonnes per 900m2 per year’ to ‘tonnes per square kilometer per year’. I know that the Resample tool, block statistics, Focal statistics tool, and Zonal statistics tool in ArcMap are useful for changing the resolution of a raster to a coarser cell size, but it will mass values (e.g.Figure 1). I want to change the units used for area value, but without loss value. But I do not know how to convert this in ArcMap.

                         Figure 1


texts: Generating statistics from line and hex grid in ArcGIS Pro?

texts:
I have a process I use QGIS/mmqgis for.  I have a hex grid that I use for binning large data, and I have that hex grid overlaying a road network.  I use the "Join Attributes by Location (Summary)" tool to generate a new road network that has a series of statistics derived from all the cells a line traverses.  Maybe more helpful to see the principle:

QGIS will give me a new road network and each feature has a min, max, mean, sum, range, etc derived from all the hex cells that feature intersects.  In the example graphic, the processed top line would have a sum of 500 and the tributary line would have a sum of 495.

I have tried to reproduce this analysis in ArcGIS Pro for a series of users who only have Esri tools, and have not been able to do so.  I have tried using Summarize Within, and this "sort of" works, but it creates a duplicate of each feature (line) for each hex cell it traverses and thus gives me no statistics.

How would I go about accomplishing this?

Pythonic solutions are fine as long as there is an explanation of which method to use in particular.


texts: Applying same color ramp to two line symbols in ArcMap

texts:
I am working in ArcMap 10.2.

I want to show two different line feature classes, one with just plain lines and the other with arrow-headed lines. The feature classes have the same categorical field I use to symbolize. I want to maintain the different line symbols, but use the same color assignments.

If I save one FC's symbology out as a layer, this (obviously) resets the actual line symbol. I've then messed w/changing Properties for All Symbols for  back to the arrow at end line but can't seem to propagate that back to all symbols in the FeatureLayer.

How do I do this?


texts:
I'm not able to find the "Diagram Overlay" plugin in QGIS 2.6 (OS X) from the Manage plugin list.

Could anyone help me to download it?


texts: OpenLayers map with custom projection not showing WMS

texts:
Trying to add to my map an WMS layer with custom projection

When executing, the map does not show any layer.

Custom projection:

proj4.defs("EPSG:25830","+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");
register(proj4);
const customproj = getProjection('EPSG:25830');

WMS Layer:

var layer =
  new TileLayer({
    source: new TileWMS({
      url: 'https://idena.navarra.es/ogc/wms',
      params: {'LAYERS': 'catastro'},
      projection: customproj
    })
  });

Map:

var map = new Map({
  target: 'map',
  layers: [layer],
  view: new View({
    projection: customproj,
    center: fromLonLat([-1.629950,42.63]),
    zoom: 9
  })
});

What's wrong here?


texts:
I am new to oracle spatial. I want to load shapefile into oracle spatial. I am using oracle 11g. I have already look at the following examples and answers

http://docs.oracle.com/cd/E11882_01/appdev.112/e11830/sdo_shapefile_converter.htm#SPATL1430

How do use the Shapefile Converter in Oracle 12c to load a Shapefile?

Where I found the following example to load a shapefile

  java -classpath
  %ORACLE_HOME%\jdbc\lib\ojdbc5.jar;%ORACLE_HOME%\md\jlib\sdoutl.jar;%OR
  ACLE_HOME%\md\jlib\sdoapi.jar
  oracle.spatial.util.SampleShapefileToJGeomFeature - h gis01 -p 1521 -s
  orcl -u scott -d tiger -t states -f states -r 8307 -g geom

I have the files on the following directory

ojdbc14.jar - D:\app\user\product\11.2.0\dbhome_1\jdbc\lib

ojdbc5.jar - D:\app\user\product\11.2.0\dbhome_1\jdbc\lib

sdoutl.jar -D:\app\user\product\11.2.0\dbhome_1\jdbc\lib

sdoapi.jar - D:\app\user\product\11.2.0\dbhome_1\jdbc\lib

Java - C:\Program Files\Java\jdk1.7.0_21\bin

shapefile (bd.shp) - E:\data\bd.shp

So I used the following codes

java -classpath D:\app\user\product\11.2.0\dbhome_1\jdbc\lib\ojdbc5.jar;D:\app\user\product\11.2.0\dbhome_1\jdbc\lib\sdoutl.jar;D:\app\user\product\11.2.0\dbhome_1\jdbc\lib\sdoapi.jar oracle.spatial.util.SampleShapefileToJGeomFeature -h localhost -p 1521 -s hr -u hr -d password -t bd -f E:\data\bd -r 8307 -g geom

Would anyone help me to run this code with those directory?


texts: Creating collection of polygons from collection of points?

texts:
I have a set of points, which are spread over a wide geographical area.  I need to be able to create a set of polygons which best represent the clustered points.

In parts, some points may simply lie on their own, but in the main they are grouped, but not linked.

I have effectively created a grid of, lets say, the north western hemisphere and in each grid, have created a point with an arbitrary value.  For all values of zero, I have deleted the points, leaving the points with non zero values in them.  I now need to turn these points into polygons and am struggling a little...

Using ArcGIS and python, I thought about clustering the points, then 'convex hulling' them, but am struggling through not really knowing ArcGIS that well; how do I cluster?

thanks in advance for any help.


texts:
I have a LZW-compressed GeoTIFF file and would like to split it into several smaller regional files and change the format to png. I tried to do the spliting  with gdal_translate -srcwin. I only get a black image. I guess that I first have to deal with the LZW-compression, but I cannot find a way to decompress files.

To compress the file I have the option -co COMPRESS=LZW. But how do I go the other way around. I am very much looking forward to your answers! :)

Here is the command I used:

$ gdal_translate -ot UInt16 -srcwin 6500 3300 1000 1000 input.tif output.tif

$ gdal_translate -of PNG -scale output.tif output.png

This is what gdalinfo says:

$ gdalinfo input.tif
Driver: GTiff/GeoTIFF
Files: EUD_CP-DEMS_2500035000-AA.tif
Size is 40000, 40000
Coordinate System is:
PROJCS["ETRS89 / ETRS-LAEA",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.2572221000027,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",52],
    PARAMETER["longitude_of_center",10],
    PARAMETER["false_easting",4321000],
    PARAMETER["false_northing",3210000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3035"]]
Origin = (2000000.000000000000000,4000000.000000000000000)
Pixel Size = (25.000000000000000,-25.000000000000000)
Metadata:
  TIFFTAG_DOCUMENTNAME=EUD_CP-DEMS_2500035000-AA.tif
  TIFFTAG_IMAGEDESCRIPTION=File written by egcs_wrgtif 2.1
  TIFFTAG_SOFTWARE=IDL 8.2, Exelis Visual Information Solutions, Inc.
  TIFFTAG_DATETIME=2013:09:17 18:46:37
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 2000000.000, 4000000.000) ( 26d54'5.12"W, 53d35'57.62"N)
Lower Left  ( 2000000.000, 3000000.000) ( 20d45'24.21"W, 45d41'42.74"N)
Upper Right ( 3000000.000, 4000000.000) ( 12d15'58.58"W, 57d13'49.19"N)
Lower Right ( 3000000.000, 3000000.000) (  8d 7'39.52"W, 48d38'23.47"N)
Center      ( 2500000.000, 3500000.000) ( 16d56'50.88"W, 51d31'45.31"N)
Band 1 Block=40000x1 Type=Float32, ColorInterp=Gray


texts: How to keep initial zeros when importing csv to dbase in ArcMap

texts:
I have a CSV table with millions of records. When I imported in into dbase via ArcMap, those inital zeros in a certain field just went away.
E.g. one field show be "01001", but it displayed as "1001".

How can I keep the inital zeros?


texts: Raster comes in black square outside my shape layer while performing unsupervised classification in QGIS

texts:
Any idea why my raster comes in black square outside my shape layer while performing unsupervised classification in QGIS? How can I correct this?

Tried correcting it by adding 100% transparency for pixel values 0, but, post processing report includes these no data values in calculation of area and pixel sum.


texts:
I've got an idea that I'd like to try in order to extend the QGIS 3 "wordwrap" label engine function.  Ultimately, I'd like to enable users to set min and max length for a single line, and then wrap the text onto as many lines as necessary.

The idea is described here, and I have been able to emulate the concept:
I can emulate the result I want using a method similar to this post:
Specifying Minimum and Maximum Word Wrap length in QGIS

However, I'd like to but, to make it easier to use,  so I was hoping to simply create/call a brand new function that looks something like:

dynamic_wordwrap ( < string_or_field > , < min >, < max >)

It seems like I'd need to start by creating a new python function in the "Function Editor" (on the "Expression Dialog").  I thought if I could find the current code behind wordwrap(), that would be the easiest way to extend it.
 However, for the life of me I cannot find the actual Python code behind the wordwrap() function.

The intellisense in the Function Editor seems to be pointing me towards "worWrap" in PyQt5.QtWidgets.QLabel, or "wordWrapMode" in PyQt5.QtWidgets.PlainTextEdit, or PyQt5.QtWidgets.TextEdit, but I can't figure out how to implement these.  Can someone help point me in the right direction?


texts: Drive time analysis with ArcGIS Network Analyst?

texts:
What is the most common SETTING for the Layer Properties when creating time-travel polygons?

I am in the process of creating 4, 5, 6 minute polygons using Network Analyst.  Every time I generate these polygons they have different shapes.

I would like to be consistent from one year to another based on the number of streets that are being added to our Center Line Street Network.

Using ArcInfo – 9.1 SP #1


texts:
I have a map application where:

I load polygons as GeoJSON objects
I draw all polygons in a single layer
I highlight a polygon on mouse hover

What I need to do is the following: on mouse hover, I need to highlight multiple polygons (I get the list of polygons somehow - irrelevant -). I have the IDs of all the polygons I need to highlight, but I can't see how can I access them though the API.


texts: Batch processing multiple files with MCC-LIDAR?

texts:
I want to classify ground points from multiple LiDAR tiles (.las files) using MCC-LIDAR software.

I already made some research and I have found this sourceforge thread; however, it is not that clear (I am not a programmer).

Any ideas how to batch process multiple files with MCC-LIDAR?


texts: Getting GDAL Color distortion around nodata value stitching?

texts:
I am trying to extract JPEGs from SID images using GDAL. I ran into a an issue where black lines were present and figured out that it was a nodata issue, so I added -src/dst/vrt nodata arguments and it stitched together properly but there are color distortions where the black lines used to be.

Here is what I am doing, using GDAL 10.1.1 packaged with QGIS:

Building the VRT from SID files

gdalbuildvrt -overwrite -srcnodata 0 -vrtnodata 0 -input_file_list C:\imaging_test\pa18.txt pa18.vrt

Warping into latlong projection (from NAD83 utm18)

gdalwarp -t_srs "+proj=latlong +datum=WGS84" -srcnodata 0 -dstnodata 0 -of vrt -overwrite pa18.vrt pa18_ll.vrt

Extracting the JPEG (subsection) from the VRT

gdal_translate -of jpeg -projwin -74.957684 40.176146 -74.880596 40.147794 pa18_ll.vrt pa_test_discol.jpeg

Here is what I am observing:

I am a programmer, and have limited experience with GIS issues.


texts: Buffering a polygon keeping the same angle/shape using GeoPandas

texts:
I have an 'area.shp' (a polygon) file that I would like to buffer with a distance of 1000m. The polygon needs to be read from a folder, and the output needs to be saved with a new name at the same folder.
@BERA's suggestion creates a circle polygon around the polygon using GeoPandas - see script below.
Can GeoPandas expand/buffering keeping the same angle/shape?
import geopandas as gpd
df = gpd.read_file("C:\GIS\data\testdata\Grid_sample.shp")
df2 = df.buffer(1000)
df2.to_file("C:\GIS\data\testdata\Grid_sample_buffered.shp")


texts: Computing volume of pits with a variable surface edge in a DEM? Multiple TINS per surface?

texts:
I have a DEM of the seafloor. The grid is highly variable, with scour pits throughout the lower right quadrant. I want to compute volume for these pits. It seems all the standard tools such as Cut/Fill using the DEM or Polygon Volume using a TIN (converted from DEM) use a set plane to represent the depth of the second surface. However, the edges of my scour pits have variable depths.

Is the best way to get a true volume to convert each individual scour polygon to a TIN, and then do a cut/fill or some TIN analysis/difference?

Is there a way to do this and have all polygons converted to TINs within a single layer? There are quite a few polygons, so separate layers and analyses would be very time intensive.

I have ArcGIS 10.5 on my Windows VM.

Related question: How to calculate volume of bund around pit from TIN and external feature edge

Looks like they opted to convert polygons to TINs, but unknown If they had to do each one separate or if they found another solution.


texts: EPSG of "Albers equal area conic projection (WGS 1972)"

texts:
I am trying to work with a raster I found here:

http://geoplaza.vu.nl/data/dataset/continental-scale-quantification-of-landscape-values-using-social-media-data

The authors report the spatial reference system is  "Albers equal area conic projection (WGS 1972)". When I load it into to QGIS the CRS is not recognized. The proj4string is:

+proj=aea +lat_1=32.5 +lat_2=54.5 +lat_0=51.4 +lon_0=22.65 +x_0=0 +y_0=0 +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +units=m +no_defs

Can anyone tell me what the EPSG of this CRS is?


texts: Why time differs between gps and ntp exactly 90 seconds

texts:
EDIT: This question can be closed due to figured out the problem was in the NTP configuration.

While the gps device is no gps data received, gps device is send gps time to server as zeros ("00000000000000") or starting of the unix time (1970-01-01 00:00:00) depending on the model of the gps device. If I get time data like this, I have to override it before record it to database as in current UTC time.

After this datas when the gps is starts to receiving gsp signals and begin to send GPS time, a problem apearing because of the 90 seconds difference between gps time and ntp time.

Say, If I get no time data from device I record data as:

2016-01-01 12:53:00
2016-01-01 12:53:10
2016-01-01 12:53:20

When device start to send GPS time records are looks like this:

2016-01-01 12:53:00
2016-01-01 12:53:10
2016-01-01 12:53:20
...
2016-01-01 12:51:30
2016-01-01 12:51:40
2016-01-01 12:51:50
...

If the device is change it's location before start sending the correct gps time, device is displayed as moved back and forth on the map.

Does this 90 seconds difference between with gps and ntp time normal and I have to make adjustment for this in my codes, or is something wrong with my system.

Additional Information:

This problem occurs once if device is newly mounted and powered until faced with power failure.
We have many different tracking devices but unfortunately I'm not allowed to give more information.
Devices are sent various information like vehicle status, location information etc. as in HEX based based data packages. I am pretty sure that I process data accurately.
I process data by splitting each information piece in data package and process each part as described in devices document.


texts: Applying categorized symbololgy to many layers using QGIS?

texts:
I have many layers that I want to symbol with same rules

Currently, I right click on each layer-->properties-->style-->
categorized-->column-->use(my case)-->classify-->
-switchboard - star symbol
-other - regular circle symbol

now this rules apply to all of the layers that I have.
switchboards get a star and other points get a circles.

it's just that I have many layers that I need to apply it to.

Is there a way to do it faster?

For each layer I use an individual color, but I use the same color for the star and the circle
e.g
if in switchboard 3 I have 24 regular points and 1 switchboard point
for the regular I use blue circle and for the switchboard I use blue star
I do that by changing the color of each symbol
and that's only for switchboard 3. for switchboard 4 I'll use the color red


texts: Automatically marking all endpoints as exception in ArcGIS topology?

texts:
I have inherited a dataset of lines representing roads, which were inconsistently connected. I have applied the Must Not Have Dangles rule to fix the topology, but the visualized errors also include endpoints. In my situation, endpoints are not errors.

I know I can mark each endpoint as an exception using the Error Inspector, but I have over 5,000 points which are as such. Is there a way I can mark them systematically, so to focus on the undershoot/overshoot inconsistencies?

A Python approach would also do for me.


texts:
My story is like this about geowebcache. I am useing geoserver 2.5.x and GeoWebCache version 1.5.2, build 1.5.x. OS is ubuntu.

I can not found geowebcache.xml file in server. But there is a file named gwc-gs.xml in .../geoserver/data/

I do not use EPSG:4326 and EPSG:900913, I use custom projection named EPSG:987654 and I set it in epsg.properties file.

And I created a Grid Set  for custom projection.

1- Should I set all layers Tile Caching settings to my Gridset?

2- When I open demo page of layer geoweb page (http://localhost:8080/geoserver/gwc/demo/pg:cimy?gridSet=EPSG:987654&format=image/png) my CPU usage is going %104 and system is locking. After 30-40 second it is normal.

What settings I should do?


texts: Error when adding single image layer from WMS server

texts:
No one question from stackExchange with similliar title works for me. Cause I face another situation. I added WMS layer, and chrome's debuger throws this Group.js?d2e6:207 Uncaught TypeError: layer.getLayerStatesArray is not a function. So map havn't displayed.
Code is simple:

import { Map, View } from 'ol';
import { Tile as TileLayer, Vector as VectorLayer } from 'ol/layer.js';
import OSM from 'ol/source/OSM';
import TileWMS from 'ol/source/TileWMS.js';
import ImageWMS from 'ol/source/ImageWMS.js';
import ImageLayer  from 'ol/Image.js';
import VectorSource from 'ol/source/Vector.js';
import { bbox as bboxStrategy } from 'ol/loadingstrategy.js';
import GeoJSON from 'ol/format/GeoJSON.js';
import { Circle as CircleStyle, Fill, Stroke, Style } from 'ol/style.js';

var imageGeoServer = new ImageLayer({
    source: new ImageWMS({

        url: 'http://mydomen.local:8080/geoserver/test_osm/wms',
        params: {
        'LAYERS': 'wspace:layer_set',
        'TILED' : true,
        },
        projection: 'EPSG:3857',
        serverType: 'geoserver'

    }),
    zIndex: 0
})

var vectorSource = new VectorSource({
    format: new GeoJSON(),
    url: function (extent) {
        return 'http://localhost:/collections/1/items/.....'+extent.join(',');

    },
    strategy: bboxStrategy
});

var vector = new VectorLayer({
    source: vectorSource,

    style: styleFunction,
    zIndex: 1
});

var map = new Map({
    layers: [vector,imageGeoServer],
    target: document.getElementById('map'),
    view: new View({
        center: [4237635.289265558, 6645635.289265558],
        maxZoom: 19,
        zoom: 6
    })
});

What is wrong?


texts: Mapping imageCollection to list of values?

texts:
I have an imageCollection where I want to create a new property called "DayOfRange" which will contain values similar to "Day of Year" except with day 1 in the middle of the year. Without going too much into details, it has to be this way so that I can make a reduction later based on this property in the metadata.

Below is the code I have tried, but it returns a list of "unknown elements" with the entire "longList" as property:

var chirpsColl = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
                                .filterDate('2004-01-01','2007-01-01')
                                .filter(ee.Filter.dayOfYear(3,28))
print(chirpsColl)
var shortList = ee.List.sequence(1,26)
var longList = (ee.List.repeat(shortList,3)).flatten()
print(longList)

var mapped = longList.map(function(i){
  return chirpsColl.set('DayOfRange',longList)
})

I was thinking perhaps I had to iterate over it? The reduceColumns also looks useful, but I'm not entirely sure how to use that with properties.


texts: AttributeError when counting entity number using PyQGIS

texts:
I'm trying to count the number of entities but each time I have this error:

`AttributeError: 'dict' object has no attribute 'featureCount'

while I have only one selected layer
uri = r"C:\impression_plan\FicheCana.shp"
layer = QgsVectorLayer(uri, "", "ogr")

layer_split = processing.runAndLoadResults("native:splitfeaturesbycharacter",
                        {'INPUT': layer,
                         'FIELD':'NOM_CONCAT',
                         'CHAR':', ',
                         'REGEX': False,
                         'OUTPUT':'TEMPORARY_OUTPUT'})["OUTPUT"]

#dissolve features based on values in an attribute table field
layer_circuit=processing.runAndLoadResults("native:dissolve",
                        {'INPUT':layer_split,
                        'FIELD':'NOM_CONCAT',
                        'OUTPUT':'C:\impression_plan\circuit.shp'})

iter = layer_circuit.featureCount()
print(iter)


texts: Georeference and warp a drone image based on sensor orientation and location

texts:
I have been flying some drone surveys over the ocean and need to properly project and georeference my images. I have all the information I think I need: lat, lon, altitude, yaw, pitch, and roll along with sensor optics params. I have to imagine this can be done with some existing package, but I can't find any. So I've been using the cameratransform package to get the GPS positions of the corners of my image:
import cameratransform as ct

# camera parameters

cam = ct.Camera(ct.RectilinearProjection(focallength_mm=f,
                                         sensor=sensor_size,
                                         image=image_size),
               ct.SpatialOrientation(elevation_m=img.altitude,
                                     tilt_deg=pitch+sensor_offset,
                                     roll_deg=roll,
                                    heading_deg=yaw))

# gps pts are lat lon
cam.setGPSpos(img.latitude, img.longitude, img.altitude)

# these are the coordinates of the image corners
coords = np.array([cam.gpsFromImage([0 , 0]), \
    cam.gpsFromImage([image_size[0]-1 , 0]), \
    cam.gpsFromImage([image_size[0]-1, image_size[1]-1]), \
    cam.gpsFromImage([0 , image_size[1]-1])])

From there I'm not certain what to do. I've tried basically saying these corners are GCPs and trying to warp the image in rasterio:
import rasterio

gcp1 = rasterio.control.GroundControlPoint(row=0, col=0, x=coords[0,1], y=coords[0,0], z=coords[0,2], id=None, info=None)
gcp2 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=0, x=coords[1,1], y=coords[1,0], z=coords[1,2], id=None, info=None)
gcp3 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=image_size[1]-1, x=coords[2,1], y=coords[2,0], z=coords[2,2], id=None, info=None)
gcp4 = rasterio.control.GroundControlPoint(row=0, col=image_size[1]-1, x=coords[3,1], y=coords[3,0], z=coords[3,2], id=None, info=None)

# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():

    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile

    # create rasterio transform
    tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])

    # I also tried this function but to no avail
    #tsfm = rasterio.warp.calculate_default_transform(rasterio.crs.CRS({"init": "epsg:4326"}), rasterio.crs.CRS({"init": "epsg:4326"}), img.size()[0], img.size()[1], gcps=[gcp1,gcp2,gcp3,gcp4])

    crs = rasterio.crs.CRS({"init": "epsg:4326"})

    profile.update(
        dtype=rasterio.uint8,
        transform = tsfm,
        crs=crs)

    with rasterio.open('example.tif', 'w', **profile) as dst:
        dst.write(src.read().astype(rasterio.uint8), 1)

This does produce an image but it is not properly warped. Since my sensor is a 40 deg off nadir the warping should be reasonably significant, when I warp it with cameratransform's cam.getTopViewOfImage() function I get this:

Though I'm not sure how to go from that to a georeferenced version or I would just use that function.


texts: Creating virtual layers using the processing toolbox

texts:
A few days ago Cy T posted an interesting question about creating virtual layers using the Python console:

Creating new Virtual Layer programmatically in QGIS?

Is it also possible to create virtual layers through the processing toolbox?

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# ----- define interface -----
##Virtual Layers=group
##Point on surface (ST_PointOnSurface)=name
##Input_polygons=vector

from qgis.core import QgsVectorLayer, QgsMapLayerRegistry

inptlyr = processing.getObject(Input_polygons)

vlayer = QgsVectorLayer( "?query=SELECT rowid AS gid,ST_PointOnSurface(geometry)AS geom FROM inptlyr.name()", "vlayer", "virtual" )
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

The log says Algorithm Point on surface (ST_PointOnSurface) finished, but nothing happens.

Edit: ?query= added


texts:
At this link https://sisteme-ig.com/questions/10129/creacion-automatizada-de-lineas-perpendiculares-entre-una-ca#232721 @Xunilk shared a Script:
mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

p_lyr = layers[0]
l_lyr = layers[1]

epsg = p_lyr.crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&field=distance:double(20,2)&index=yes"

dist = QgsVectorLayer(uri,
                      'dist',
                      'memory')

QgsMapLayerRegistry.instance().addMapLayer(dist)

prov = dist.dataProvider()

lines_features = [ line_feature for line_feature in l_lyr.getFeatures() ]
points_features = [ point_feature for point_feature in p_lyr.getFeatures() ]

feats = []

for p in points_features:

    minDistPoint = min([l.geometry().closestSegmentWithContext( p.geometry().asPoint() ) for l in lines_features])[1]
    feat = QgsFeature()
    feat.setGeometry(QgsGeometry.fromPolyline([p.geometry().asPoint(), minDistPoint]))
    feat.setAttributes([points_features.index(p),feat.geometry().length()])
    feats.append(feat)

prov.addFeatures(feats)

But from that QGIS throws me the following error:
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.12\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 5, in <module>
TypeError: index 0 has type 'QgsPointXY' but 'QgsPoint' is expected

How can I change this Python code so it doesn't give this error?


texts: Extract band values to points from mosaic (ArcGIS) in R

texts:
I am classifying satellite images in ArcGIS. I would like to partially evaluate my results in R, e.g. through Crossvalidation. I want with mosaic (raster with 6 bands) extract band values to points but from my raster it only gets values from 1 band. How to get values from each bands?
Below is the code I used:
df <- mosaic %>%
extract(y = points) %>%
as.data.frame %>%
mutate(id_cls = points@data$id_cls) %>%
left_join(y = unique(poly@data), by = c("id_cls" = "id")) %>%
mutate(id_cls = NULL)

I would like to receive data.frame with the values of all bands at selected points.


texts:
Working with point data from Texas and am unfamiliar with what coordinate system it's using. Unfortunately, I am on a work computer which doesn't allow me to do much without admin privileges, so downloading and running Projectionuesser isn't a possibility at this moment.

Below is some examples of the point data:

2127838.753,    13709652.5

2061615.337,    13718121.64

2061615.337,    13718121.64

These points should all be in the San Antonio area of Texas.


texts:
I'm trying to get detailed geojson data for various African countries. I think I found a good source at http://download.geofabrik.de/. The problem is that the downloaded file is in osm format. I know of the osmtogeojson tool, but there is no way this would work on the 32.2GB (!) file I got from Geofabrik. I also have the pbf file which is 1.2 GB.

My question is:

Is it a way to split an OSM file (in my example, I could extract data country by country and hopefully feed them to osmtogeojson then)
Is it any other way to get detailed geojson info at the subnational level (provinces)?


texts: Converting GeoTIFF file to numpy array to QPixmap with PyQt4?

texts:
I'm a student in CS doing an internship and I'm not familiar at all with GIS yet.

I have to display a GeoTIFF image (with a single band) into a QPixmap with PyQt4. So, I found some code that could help me in the GDAL library, but I think my conversion from my GeoTIFF file to a numpy array isn't correct. Here's why:

When I display my QPixmap from a QImage (see code below), the only thing I have is a transparent band at the top and the rest of the image is just a huge black rectangle.
When I type gdalinfo <FILENAME> in my console, it tells me that the maximum value is 165.146. However, when I get the max value from my numpy array, I get 544.843.

Here's the code:

import numpy as np
from osgeo import gdal
ds = gdal.Open("path/to/file.tiff")
mr_array = np.array(ds.ReadAsArray())
maxes = []
for sub in mr_array:
    maxes.append(np.amax(sub))
print(max(maxes)) # => 544.843
mr_img = gray2qimage(mr_array) # comes from this script: https://github.com/sjara/brainmix/blob/master/brainmix/gui/numpy2qimage.py
self.pxm_mr = QtGui.QPixmap.fromImage(mr_img).scaled(371,301,QtCore.Qt.KeepAspectRatio)
self.lbl_pxmap_mr.setPixmap(self.pxm_mr)

The idea was to search where the problem was coming from: the GeoTIFF -> numpy array conversion, the numpy array -> QImage or the QImage -> QPixmap conversion. My guess is that it comes from the GeoTIFF -> numpy array because of the different maximums, but every code I find about that kind of conversion is basically the same as mine.

Am I right in thinking it comes from here?

If so, how could I correct this?

If not, where is my problem most likely to come from?


texts: Use Windows function GetProductInfo via MapBasic

texts:
I develop a MapInfo-based GIS and I need to be able to determine the Windows OS of a given user. So far we have used a Windows function GetVersionEx (kernel32 library) to do this but the function doesn't behave properly from Windows 8.1 onwards...unless you do Application Manifesting which I'd rather not mess around with unless absolutely necessary.

As an alternative I found a post suggesting to use GetProductInfo instead but I don't know how to set this up in MapBasic.

The GetVersionEx function is first defined like this:

Type OSVERSIONINFO
    dwOSVersionInfoSize As integer
    dwMajorVersion As integer
    dwMinorVersion As integer
    dwBuildNumber As integer
    dwPlatformId As integer
    szCSDVersion As String * 128
End Type

Declare Function GetVersionEx Lib "kernel32" Alias "GetVersionExA"
    (lpVersionInformation As OSVERSIONINFO) As integer

(OSVersionInfo documentation here)

And then called in MapBasic:

Dim OSInfo as OSVersionInfo
Dim nRet,iMaj_Version As Integer

'Set the structure size
OSInfo.dwOSVersionInfoSize = 148    'Len(OSInfo)

'Get the Windows version
iRet = GetVersionEx(OSInfo)

'get values e.g.:
iMaj_Version = OSInfo.dwMajorVersion

so it will be similar to this I expect...but I don't know how the new function needs to be defined or how I then get the output in to a MapBasic variable.

Finally, we are using the 32-bit version of MapInfo so for the time being I'm not interested in potential help/solutions for the 64-bit version.


texts: PyQgs tool that replaces layout text labels with attributes from an input table

texts:
In the Qgs Project is a csv table file with two fields and two features/Objects.

The script I'm writing is supposed to replace two textlabels in my only open layout with the attributes in the second field of the table.

My script has the following procedure:

1) take the table from the project and store it in a variable

2) take the text labels from the layout and store it in a variable

3) compare the text labels with the attributes from the first field and if they are equal overwrite the text label with the content from the second field.

My script doesn't give me an error message but it is not doing anything. apparently I missed something but I don't see what it is. Hoping for help.

Here is the code:

    def getTable(self):
    """gets table from the QGIS Project"""
    layer = None
    for lyr in QgsProject.instance().mapLayers().values():
        if lyr.name() == "table":
            layer = lyr
            break
    return layer

def setVariable(self):
    """adds the table to a variable called self.inTable"""
    self.inTable = self.getTable()

def changeLabel(self):

    """saves the Layout label texts in a variable"""
    labels = None
    for l in QgsProject.instance().layoutManager().layouts()[0].items():
        if isinstance(l,QgsLayoutItemLabel):
            labels = l

    """gets the features of the table"""
    features = self.inTable.getFeatures()

    """compares the list of labels with the attributes in the first
    column of the features and replaces the label with the attribute
    from the second column"""
    for label in labels:
        for feat in features:
            if label == feat.attributes()[0]:
                labels[label] = feat.attributes()[1]
    return labels

def run(self):
    """Run method that performs all the real work"""
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        self.setVariable()
        self.changeLabel()

I guess there is something wrong in the changeLabel function but I don't see what it is. Or I completely forgot somethin.


texts: Is there a way to retrieve shapefile attribute table information for a given coordinate in R?

texts: Does anyone have an idea how could one get information out of shapefile attribute table (class SpatialPolygonsDataFrame in R), in my case data from column that contains species names, based on a list of coordinates that I have. Each coordinate that is within a certain polygon should be given the species name, and as a result I would like to have a dataframe with three columns: X,Y,species name. Are there any packages and their specific functions for this kind of task? Maybe sp package, aggregate function?


texts: Sending coordinates of features of any polygon featureclass to a table

texts:
I am trying to store x-coordinates and y-coordinates of a polygon feature in a table. But, i am getting an error message: Item not found in this collection. The code is like this:

Do
 pFeature = pFeatureCursor.NextFeature
 If Not (pFeature Is Nothing) Then
   pPointCollection = pFeature.Shape
   For lPnt = 0 To (pPointCollection.PointCount - 1)
     ppolyrow = ptable.CreateRow()
     ppolyrow.Value(coorxind) = pPointCollection.Point(lPnt).X()
     ppolyrow.Value(cooryind) = pPointCollection.Point(lPnt).Y()
     ppolyrow.Store()

The error comes in this following line of code:

ppolyrow.Value(coorxind) = pPointCollection.Point(lPnt).X()

Any help....


texts: Python : How to plot a raster .tiff image as lines instead of colours on a map?

texts:
On python 3.8.2, Windows 10.

I'm currently plotting a monthly mean temperature map from a .tif image on a basemap with the gdal package (see code). I'm following a method that I've seen online.

I'd like to add, on top of this mean map, a map of variance that can be seen as some sort of shaded grey lines. Bigger grey lines would indicate higher variance and lower gray lines would indicate smaller variance. Although, perhaps it would be better to, instead of having the strength of the variance show as line size, show the variance as line transparency. I'd be curious to try those out, but I don't even have any idea on how to show a .tif image as lines!

Here is an example of what I have in mind. The areas with no bars basically have little to no variance.

Here is the code used for plotting the map (considers the mean data only).

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
from mpl_toolkits.axes_grid1 import make_axes_locatable
from osgeo import gdal

## Read the data and metadata
ds = gdal.Open(Tiff_path + Tiff_Filename)
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()

xres = gt[1]
yres = gt[5]
# get the edge coordinates and add half the resolution
# to go to center coordinates
xmin = gt[0] + xres * 0.5
xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5
ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5
ymax = gt[3] - yres * 0.5
ds = None

# create a grid of xy coordinates in the original projection
xy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
xy_source = np.delete(xy_source,(-1),axis=1)

#subtitute -999 to nan
data[data == -999] = np.nan #Set all -999 to nans
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
# set Basemap with slightly larger extents
# set resolution at intermediate level "i"
m = Basemap(epsg=3857,
            llcrnrlon=xmin, \
            llcrnrlat=ymin, \
            urcrnrlon=xmax, \
            urcrnrlat=ymax, \
            resolution='i')
m.drawcoastlines(color="gray")
# m.fillcontinents(color='beige')
im1 = m.pcolormesh(xy_source[0,::], xy_source[1,::],data.T,alpha=0.6,cmap="jet",vmin=Min_val,vmax=Max_val,latlon=True)


texts: Placing a WMS layer above a WFS layer makes WFS layer not selectable

texts:
I have two different layers which I call from Geoserver. The first one is a WMS layer with points and the second one is WFS layer with lines.
I need to display the points above the lines and I also need to be able to select each line and get some data.
I use the Openlayers2 library to display the layers on the map.
I set the index of the WMS layer like this:

wms_layer.setZIndex( 1001 );

So I am sure it will appear on top when I display it from the layer switcher.
But the problem is that when I make visible the WMS layer then the WFS layer is not clickable.
Is there a way to fix this?


texts: Seeking a method to print all visible layers in QGIS (3.10), using Atlas (and a map theme)?

texts:
Hope to obtain: sample expressions for printing each visible layer of a project using Print Layout's Atlas functionality
Based on @Erik's comments on another question (about changing QGIS behavior), I'm trying to follow through on setting up printing an "atlas."
The dialog (below) is imposing for a neophyte. The first reference returned from a search on 'atlas' in the QGIS docs is to an Example, which nevertheless is enough different from my scenario as to not shed light.
It focuses on using a vector layer composed of an array of rectangles to print a set of smaller maps covering a single layer.
In my application, I wish to use a layer with a single rectangle to print a set of maps for the same region from each of a set of 4 layers. (Erik had also suggested perhaps using a Map Theme to reference all the visual layers, but I do not find any reference to a map theme among the items available in "Expressions.")
Beyond opening Atlas->Settings in Print Layout, selecting "Generate an atlas", choosing my "Coverage layer" (with the vector layer containing my polygon), and hiding that coverage layer (so the polygon won't print on the outputs), I'm lost.
In particular, I'd like to know how to fill out:

Configuration->PageName, ->FilterWith, and ->SortBy fields
Output->OutputFilenameExpression field.

I'd also like to know how to get a list of the visible layers and how to step through them in Expressions Dialog (which comes up when clicking on the side-button of the PageName field). If using a map theme is the way to go, how is it referenced in the Dialog?
Thank you.


texts:
I am currently using QGIS and it's database tool with custom plugins. I have to update a specific Table using an SQL query but I cannot use Database specific programs like PostGIS. Here is the test query I used:

UPDATE  "BAL"
SET
"BAL"."BAL_ADRESSE" = "BAL2"."bal_adresse"
FROM "BAL", "BAL2"
WHERE "BAL"."BAL_ADRESSE" != "BAL2"."bal_adresse";

It is quite a simple query, I am trying to replace values in the first BAL table using a csv I created myself and imported in the same database afterwards, calling it BAL2. Yet everytime I try to do anything than basic selections, like the update I need to do, it returns an error message and never the same.

So here are my questions: Is my query wrong? And is QGIS unable to use the UPDATE query?


texts: Numpy Logical Syntax in QGIS GDAL Raster Calculator

texts:
I need to set a logical expression in the box of the GDAL's Raster Calculator using Numpy. My knowledge in Numpy and python is extremely limited (I'm not a Developer just a simple Geographer!). I've seen on this post (Using variable in graphical modeler for raster calculation? (QGIS 3.8)) an example of how to append 2 parts in that box, the part using the letters (A, B, C...) and the other from other inputs like numbers. In that example is easy ('a + ' || @numberinput). My question is about logical statements. I tried 'a >=' || @numberinput and the result is wrong. I checked the numpy documentation and all is about arrays which is notice is kind of complex, I don't need to do too much, just get those A values less or greater than the result from a formula expressed as string getting one variable number from the @numberinput.   I'm attaching the images of what I'm doing for you to help me. Threshold set using String concatenation, the output will be used in GDAL Raster Calculator

Then inside GDAL Raster Calculator you can see my issue:

Finally, you can see what I'm trying to do inside the Expression window.

I think this is an easy question for python experts.


texts: How to get Openlayers cluster strategy work for large data?

texts:
The map Layer has around 1500 points. I am able to display all the points in the map without any problem but adding the cluster strategy to the layer does not work. The Cluster strategy functions when I limit the number of points to 100.

Is there a way to cluster large amount of points in Openlayers?


texts: Remove Duplicate Records IDENT field -Shapefile using GDAL OGRINFO or OGR2OGR

texts: I have two shapefiles that I am merging/combining. However, this produces more than 4,000 duplicate records.  I need to have only unique records.  Is this possible with OGRINFO or OGR2OGR?   I've tried Select DISTINCT but it only lets me do select statement on one column/field.  I can get a list of unique records but I need to create new shapefil.


texts: Error in extracting Raster Values Greater than 0.5 in QGIS

texts:
I want to extract raster values which is greater than 0.5 .I Used Raster Calculator to create a new raster which have values greater than 0.5 .

Formula for Raster Calculator:

("raster" > 0.5) * "raster"

but it create a tiff having values from 0 .
Now I want to change 0 to No Data

For this i tried Raster--> Conversion --> Translate(Convert Format)

gdal_translate -a_nodata 0 -of GTiff input.tif output.tif

But output does'nt have correct values.

Is there any other way to do this by Python Console or by GDAL .


texts:
I have 4 feature types and writing to a single xls file using FME. Each feature type has its own worksheet.

How can I sort the worksheets with FME before creating the xls-file?

I use FME 2012 SP3


texts:
I am installing Mapnik v2.2.0 on my 64 bit Windows with (64 bit) Python 2.7. After following the usual install guide, when trying to import mapnik I get the following error:

Traceback (most recent call last):

File "<ipython-input-9-1290ece2bb4b>", line 1, in <module>
import mapnik

File "x\mapnik-v2.2.0\python\2.7\site-packages\mapnik\__init__.py", line 69, in <module>
from _mapnik import *

ImportError: DLL load failed: %1 is not a valid Win32 application.

This definitely sounds to me like a 32 bit-64 bit issue, and Mapnik does say that its installation is for 32 bit. However, Googling around seems to suggest it should also work on 64 bit installation (as e.g. none of the guides explicitly say "this will absolutely not work if you're on 64 bit).

What would you think is the issue? Should I install 32 bit Python on my 64 bit OS? I would rather treat that as a last resort.


texts:
I want to symbolize a line layer like this:

With three dots - empty space (not a white covering line, since i have other data underneath) - three dots.

The closest I come is using "hashed line" and rotating the hashes 90 degrees, but it does not look good, the dots are short lines (rectangles, not circles)

How can I create the symbol? Maybe using geometry generator?


texts: Need steps to print a Map in ArcGIS Silveright

texts: I want to know how to print a map in ArcGIS Silverlight application.


texts: How is NokiaHere Geocoding Score calculated?

texts:
When geocoding an address in NokiaHere, a Component Score is calculated according to the information provided in the original address.

Example:

  {"componentsScore": {"CITY": 1.0, "COUNTRY": 1.0, "DISTRICT": 1.0,
  "HOUSE_NUMBER": 1.0, "POSTAL_CODE": 0.0, "STATE": 0.0, "STREET": 1.0},
  "totalScore": 0.92}

If you could provide any information on the way this estimation is made, it would very helpful.


texts: Can collector auto populate related feature ID field?

texts:
I have staff conducting hydrants inspections using ESRI ArcGIS collector. The hydrants have a unique ID. If I make a relationship class between hydrants and hydrant inspections, does Collector know to auto populate the hydrant inspection facility id with the hydrants unique ID?

The hydrant inspection layer is a feature class, not a table.

I know collector is able to auto populate related tables; however, what about related feature classes.


texts: Add shapefile from a local file to a JS API map

texts:
Using ArcGIS.com's Map Viewer it's possible for the user to add a local shapefile to the map, using an Add Data approach.

Is it possible to add this functionality to an ArcGIS Server map built on the JavaScript API?

(I posted this on the ESRI forum but haven't had a reply yet. I'm hoping the GIS.SE experts will know the answer)


texts: How to tell if a device supports Assisted GPS, and if how to get the Assisted GPS cooridinates?

texts:
1)  Is Assisted GPS implemented at the hardware and OS level, or is it a hardware and secondary software thing.  I've never used this so wondering what is required to 'access' this.  I noticed a wiki article so is seems it is very device specific, but wondering if this is full and complete list.

2) How are the aGPS value accessed?  Are they streamed through a diff port (v. true GPS)?  Is that configurable?

I tried to ask appropriate questions to understand more about this, but I'm kind of lacking context so my questions may need some edits from some who knows more about this.

thanks
df


texts: Google Street View: Getting photos from different years

texts:
Given a pano ID or longitude and latitude of a location, I want to obtain all of the photos of that location in different years.
All I could do now is to get all of the photos that are linked to that pano ID using geo2.ggpht.com.


texts:
I know there is Google API for distance where we can calculate the driving distance but similar to that I want to pass two city names and I want it to give me just the duration like 4 hours 30 minutes.

Is there anyway I can achieve it. I want to pass my values into the API and get the travel times.


texts: Creating .lyr file to import symbology using ModelBuilder?

texts:
I am trying to create a .lyr file in order to import symbology to a .shp file. The ESRI link says to use ModelBuilder to do this, but I am having trouble understanding the directions bc I haven't used ModelBuilder in forever....

http://desktop.arcgis.com/en/arcmap/10.3/analyze/modelbuilder/creating-a-layer-symbology-file.htm

I got as far as the attached photo, but am unsure of actually how to save the file. It mentions 'adding to display' but when I do nothing comes up.


texts:
What is the best way to validate SLD?

I'm using geoserver-manager but i cannot find any method to validate SLD. Even if i publish invalid SLD over geoserver-manager the response is 200 OK but it is not saved.

Of course I'm asking for soultion over geoserver rest api.


texts: Removing mid-division label on subdivisions before zero on ArcGIS Pro scale bar?

texts:
Using ArcGIS Pro 1.3.1 I have created a scale bar using the Scale Line 1 Metric template with a few modifications like:

Show one division before zero
Divisons: 11
Subdivisions: 10
Label Position: Below center
Label: KILOMETRES

It appears as below:

How can I remove the mid-division (0.5) label (circled above) on the subdivisions before 0?


texts: Clustering method based on building attribute data

texts: I have a building footprint and each building has a quality metric attached to it from a scale 1 to 6. I want to know if there is a certain pattern or cluster associated with the building spatial distribution. I am using ArcGIS Pro and there is a lot of clustering methods with different parameters there, but I don't know what would be the best method.


texts: Count the number of points in each Voronoi polygon

texts:
I have a Data.frame consisting of Longitude and Latitude called Points1, with which I use Voronoi polygon.
I also have another Data.frame of Longitude and Latitude called Points2.

I then want to count the number of occurences of Points2 in each Voronoi Polygon. Regarding the plot, I want to count the occurences of triangles in each Polygon.

Sample code:

library("deldir")
Points1 <- data.frame(x=runif(10), y=runif(10))
Points2 <- data.frame(x=runif(100), y=runif(100))
x=Points1$x
y=Points1$y

vtess <- deldir(x, y)

plot(x, y, type="n", asp=1)
points(x, y, pch=20, col="red", cex=0.5)
points(Points2$x,Points2$y,pch=2)
plot(vtess, wlines="tess", wpoints="none", number=FALSE, add=TRUE, lty=1)


texts:
I am using from gps for get latitude and longitude of user when is walking, but when I see points(latitude and longitude) it get me points with 2 or 3 meter distance difference.How can I have the exact coordinates .

I am writing my app with android studio.

Please see too:

Link


texts:
How can I adapt the following expression in the ArcGIS field calculator? Below is the python code:

import arcpy
fc = "2014_2"
keyField = "ID_segment"
affectedField = "KKN2"
maxVal = 2500
cnt = 1
Value = -1
oldValue = str(maxVal)
queryString = '"' + keyField + '" < ' + "" + oldValue + ""
rows2 = arcpy.UpdateCursor(fc, queryString)
row2 = rows2.next()
while row2:
        iValue = row2.getValue("ID_segment")
        if iValue != Value :
                cnt = 1
                Value = iValue
        else :
                cnt = cnt + 1
        kVal = str(row2.getValue("KKN")) + "/" + str(cnt)
        print(Value,kVal)
        row2.setValue(affectedField, kVal)
        rows2.updateRow(row2)
        row2 = rows2.next()


texts: What does "Information [18]: Select at least one signature" in the semi automatic classification plugin mean?

texts: I am using the semi automatic classification plug in to classify vegetation and water in a satellite image. I create training shp and signature list. When i hit "classification preview" it gives me the error "Information [18]: Select at least one signature". I tried selecting one but it did not help. What is meant be "Information [18]: Select at least one signature"?


texts:
I'm starting a project in which I intend to use GeoServer/GeoWebCache as a provider of WFS, WMS, WCS, and WPS services. I've configured the GeoServer to use PostGIS for storage. The frontend is using OpenLayers.

I, however, would like to use C# to program my business logic and to send requests for WMS, WFS, WCS and WPS services. I then make the response available to my view layer using View Model. I've been looking for materials that explain how this can be achieved, but no success so far. Any advice will be highly appreciated.


texts:
What is more correct, km2 or sqkm?

I'd prefer to use km with a superscript 2, but it's not always possible to format text that way so I was just wondering if there is a convention that should be used when you can't use superscript formatting.


texts: Changing Layer Definitions / Refresh Not Behaving in ArcGIS API for JavaScript?

texts:
In ArcGIS Javascript API 3.3, I have an ArcGISDynamicMapService layer in my application that points to a fairly simple (6 layers) map service.  In the application, there is a GUI/widget that allows the user to pick various options - which then dynamically changes the layer's layerDefinitions property to re-render the layer based on their selection.

After I define and set the layerDefinitions of the layer, I call the layer.refresh method to re-render the layer.  However, the re-rendered layer is always wrong - either the wrong layers are rendered, or wrong features are rendered, and there doesn't appear to be any pattern in these errors that I can determine.  However, once the map is panned or zoomed, the re-rendered layer is drawn exactly as expected.  Below is the relevant javascript code.  The "format" function is a .NET-like string formatting function I've added to the script.

var layerDefinitions = [];
layerDefinitions[1] = "TransmitterId IN ({0})".format(visibleTransmitterIds.join(","));
layerDefinitions[2] = "TransmitterId IN ({0})".format(visibleTransmitterIds.join(","));
analysisResultLyr.setLayerDefinitions(layerDefinitions);
analysisResultLyr.refresh();

Using various http sniffing tools (Google Developer tools is preferred choice), I can see that the outgoing http request for my layer of interest is exactly the same between the refresh call - and the pan/zoom call.  I have also verified the sql contained in the layer definitions, and am positive the layers/features returned in the refresh call are definitely wrong.

Anyone have any ideas as what may be causing this?  Better yet, does anyone know an alternate way to refresh a layer without panning/zooming the map?

A couple of other notes:

I have disabled cacheing at the client via setDisableClientCaching
My sql where clauses make use of IN statement.  Is this a bad idea with file geodatabases?