texts: How can I compute confusion matrix using Python

texts:
My project is about "Extracting water and urban features from satellite images". I have extracted urban and water features using Python code.
Now, I want to do accuracy assessment. I know two primary requirements for accuracy are reference data and predicted data. So, as reference data I have GeoTIFF file and as predicted data I have my extracted image (GeoTIFF format).
I tried with accuracy assessment using Python but I am getting accuracy as 99% which is not possible.
Can somebody suggest how I can achieve this using Python?

Please find my code:

import numpy as np
import matplotlib.pyplot as plt
import cv2
from PIL import Image
import os
import gdal
import osr

gtif = gdal.Open('E:\\tanvi_iirs_work\\Task 3\\Images\\geosample.tif')
srcband = gtif.GetRasterBand(1)
trans = gtif.GetGeoTransform()
proj = gtif.GetProjection()
#print(proj)
#print("1")
arr1=srcband.ReadAsArray()
arr = srcband.ReadAsArray()
print("urban area")                       #Extracting urban area
for i in range(arr.shape[0]):
    for j in range(arr.shape[1]):
        if 0.6 <= arr[i,j] <= 1:
            arr[i,j] = 1
        else:
            arr[i,j] = 0
plt.imshow(arr,cmap = 'gray')

cols = arr.shape[1]
rows = arr.shape[0]

driver = gdal.GetDriverByName('GTiff')

if arr.ndim == 2:
    band_num = 1
else:
    band_num = arr.shape[2]

outRaster = driver.Create('E:\\tanvi_iirs_work\\Task 3\\Images\\outputtrial1.tif',cols,rows,1,gdal.GDT_Float32)
outRaster.SetGeoTransform(trans)
#srs = osr.SpatialReference()
#s1 = srs.ImportFromEPSG(4326)
#print(srs)

outRaster.SetProjection(proj)
#print(outRaster.GetProjection())

outband = outRaster.GetRasterBand(1)

outband.WriteArray(arr)
#data2 = outband.ReadAsArray()
#print(data2[345][456])
outband.FlushCache()
outRaster=None

#Accuracy Assessment

tp=0
fp=0
tn=0
fn=0

for i in range(arr.shape[0]):
    for j in range(arr.shape[1]):
       if arr[i][j]==1 and 0.5=<arr1[i][j]<=1:
           tp+=1
       if arr[i][j]==0 and 0.0=<arr1[i][j]<=0.4:
            fp+=1
       if arr[i][j]==1 and 0.5=<arr1[i][j]<=1:
            tn+=1
       if arr[i][j]==0 and 0.0=<arr1[i][j]<=0.4:
            fn+=1
print(tp,tn,fp,fn)
acc=(tp+tn)/(tp+tn+fp+fn)
print(acc*100)


texts: Extending Identify widget and connect it with Legend widget

texts:
I am using latest Flex Identify widget and I have created my own legend widget.

What I am trying to do is, when a user select a point using Identify tool it selects some feature on the map. That feature refers to some common value in our legend. What I want is to highlight that value in legend widget.

For further explanation:

If within a feature we have 6 attributes coming from our layer and we have the same 6 features present in Legend widget, I want to connect/enhance the Identify widget so on point-click I can see the widget on the map as well as its properties on the Legend Widget.


texts: Why the footprints of the identical Sentinel-2 image from GEE and SciHub are different?

texts: The footprint of one Sentinel-2 tile (e.g., the image whose ID is "COPERNICUS/S2/20200702T023551_20200702T024443_T51RTQ") from Scihub is rectangle, while that from GEE is polygon. How to explain this difference?


texts: Degrees Minutes delimited Seconds conversion

texts:
How do you handle delimited seconds?

I can find a bunch of stuff on converting degrees minutes seconds. That is easy enough. My question is what about when there is a 4th set of numbers for example, 41 39'51.91N, 71 15'00.47W.

What do I do with the 91 or 47?


texts:
I am trying to download NDVI data from the Google Earth Engine MODIS satellite. I used the following dataset:

  img = ee.ImageCollection("MODIS/006/MOD13A2")

To actually download it, I ran img.getDownloadUrl() to produce a URL and retrieve the data.

According to the website https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD13A2, it states that NDVI has a minimum value of -2000, a maximum of 10000, and a scale factor of 0.0001, inorder to scale the NDVI between the standard values of -0.2 and 1.0.

However, the images downloaded from ee.ImageCollection are tif grayscale images, and using cv2.read() on those images produces a matrix of values between 0 and 255 inclusive. Can anyone explain what the true NDVI values of the image are? Furthermore, most of the values are between 0 and 35, while water is returning values of 250+. Im confused as to what is going on? How is the data that I am downloading being scaled?

Update: I have the code below

print(img.projection().getInfo())
path = img.getDownloadUrl({
    'crs': 'EPSG:4326',
    'crsTransform': [926.625433056, 0.0, -20015109.354, 0.0, -926.625433055, 10007554.677],
    #'scale':1000
    'bands': [{'id':'NDVI'}],
    'region': '[[' + str(currentMinLng) + ',' + str(currentMaxLat) + '],['+str(currentMaxLng) + ',' + str(currentMaxLat) + '],['+str(currentMaxLng) + ',' + str(currentMinLat)  + '],[' + str(currentMinLng) + ',' +str(currentMinLat) + '],[' + str(currentMinLng) + ',' + str(currentMaxLat)+ ']]'
    #'region': '[[16.499050771335575,-29.07796275949132],[18.584367695114103,-29.07796275949132],[18.584367695114103,-30.882],[16.499050771335575,-30.882],[16.499050771335575,-29.07796275949132]]'
    #'region': '[[20.2431,-28.882],[22.2431,-28.882],[22.2431,-30.882],[20.2431,-30.88]]'
})

In the above code, the "img" variable is a Google Earth Image from the "MODIS/006/MYD13A2" image collection.The variable "path" produces a download link to the image. When I open the image using opencv and read the .TIF file, I get a bunch of pixels that are oddly scaled. What could cause this to happen?

I realized that it may not be the .TIF thats causing the problem but rather the way I am reading it. What would be the correct way to read in this .TIF using python?


texts:
I am trying to style a new single-band raster with a platted style, but when I choose the ‘platted’ no colors or values appear, and I can’t find a way to add them.

QGIS Essen, Win 7.


texts: Adding a delimited text layer in QGIS 1.8 only imports a fraction of the total rows

texts:
I am trying to add a delimited text layer, and I think I've formatted it correctly as a CSV file. There are no errors when I import the data. But even though there are 1300 rows total, only 189 appear in the table when I right click on the layer and open the attribute table. Is there any reason why it would import so few rows?

Thanks!


texts: Using Collect Values in ArcGIS ModelBuilder with submodel?

texts:
I have a large model with submodels.

Basically the Model iterates through all Districts, selects in each district an attribute above a certain threshold, creates a layerfile of it. The submodel iterates through all the Attributes of the created layerfile and then selects the neighbor of this attribute below a certain threshold and creates another layerfile of it.

In the end I want to merge the layerfiles. I though of using a collect value tool to collect all the layerfiles (red circle) which got created in the iterations in the submodel and then add a Merge tool to it in the main model ( or should it be in the submodel)? However, I get no result if I add it to display.

The whole Model works if I just run it and see the result for the last run in the iterator.

Can anyone help me with this? I attached the model but the jpeg probably doesn't make any sense as the text is too small.

I also asked the question on Geonet: https://geonet.esri.com/message/609161


texts: Is there a way to test HERE Maps offline?

texts: I am developing web apps that implement HERE Maps. I would like to ask if there is a way to test it out offline (ie. in localhost or offline development machine)?


texts:
I am trying to create groups based on a primary object as its centroid.

Groups are defined as all objects within a certain distance from a primary object. I am currently using spatial join in a ModelBuilder to create the groups (one to many) with the defined distance, but this causes objects to join more than one primary object's group.

Many objects can be joined to one primary object but objects cannot be in more than one primary object group. If an object falls within the buffer of two primary objects, I would like the object that is closest to be added in that group.

I can use ArcPy if necessary.


texts: I have 3 shapefiles(shp) having containing the same attribute names. How do I merge these 3 shp files into one shp file in QGIS 3.0.0 desktop. I don't have the option of Merge Shapefiles to One under Data Management tools. Can someone help me?


texts:
i have a problem with "tile layers". i want to display my layers in tile layers but i always have a blank page. i'm already displaying them in the "Preview layer" with openlayers and it worked. i don't know where is exactly the problem. i'm just a beginner in geoserver so if you have some useful tutoriel it will be great
thanks guys


texts: Mosaicing many large ECW files for display in ArcGIS Server apps

texts: I have about 200 ECW files of recently acquired 1-foot resolution aerial photos. Its about 15GB worth, covering over 250 sq. miles. Unfortunately, you need a license from ER Mapper that runs about $10K to publish ECWs in ArcGIS Server. I brought them all into GlobalMapper, mosaiced them, and pushed out a jpeg2000 that is about 31GB in size. It actually draws quite well, looks gorgeous, and is pretty snappy. Anyone have any other suggestions on how to deal with imagery datasets of this size (and by size I mean both areal extent and filesize)?


texts:
I need to know which points from Table 1 are contained in Table 2.

In Table 1 I have longitude and latitude and I made points with ST_SetSRID(ST_MakePoint(d.longitude, d.latitude),4326) as geom

Table 2 has a collection of multipolygons: geom geometry(MultiPolygon,4326)

I know it is a collection because when I run
select ST_IsCollection(geom)from table2 it is TRUE, but when I want to know how many collections with SELECT ST_NumGeometries(geom)from table 2, it returns 1.

I know that since I have a collection I have to use st_INTERSECTS.

select TABLE1.*, TABLE2.*
from TABLE1, TABLE2
where ST_Intersects(TABLE2.geom, TABLE1.geom)

But it is empty. And I know from plotting the data that they intersect.

Do you know what might be wrong?


texts: Using Field Calculator to assign consecutive letters for each feature?

texts:
I have a feature class generated in ModelBuilder using a transect script that contains multiple transects.  Before moving on to further build the survey I would like to assign ID's to each feature in the class using field calculator.  It needs to be string, i.e. A,B,C....and so on.  After 26 features I would like to name it AA, BB, CC and so on.

Does anyone have a script that could accomplish this?


texts: Displaying row/column of point in raster using QGIS?

texts:
Is it possible in QGIS to display column and row index of pixel in raster?

I mean something similar to value tool.


texts:
I will try explain this as best I can.

I have a point on a map and I want to create buffers surrounding the site; for example, 10Km, 20Km 50Km.  This is easy.  But I also need to show buffers for distance from site, based on road network. lets says 30mins, 60mins 180mins etc.

I have a detailed road network map showing the various road classes, and I could assign average speeds to them but from here... I am lost.

Is this possible in QGIS 2.6.1.

Also if I could use Googles route calculator it would be AMAZING!

Oh and I have tried the route calculator tool, but I want to put in point a, then identify max travel time, then show polygon of area.


texts:
i have used gdal2tiles.py for tiling about 100 tiff tiles which i have generated virtual dataset with gdalbuildvrt for my 1:250.000 scaled maps...

everything goes well from 15 to 11, but when it comes to 10, the following parts was not generated from top to downs.(I gave a picture of the sample below) and it finishes at 70 like this ...50...60...70...

i think gdal dont get extent of map for making it small(resolution) for zoom levels.

is there anyone meet this problem or somebody help me about this trouble?

Example Result:

Black areas indicate regions which wasnt rendered.


texts: Getting elevation readings at any lat,lon from shapefile with scattered irregular points

texts:
I have a set of x,y coordinates, each of which I'm trying to get an elevation reading for.

The elevation data source, however, is not a regular grid, but consists of a shapefile containing scattered, irregular points, as this QGIS screenshot shows:

I'm used to using grid-based files such as SRTM or ASTER, which have a reading on a regular grid, and using a TIF reader to extract the locations. However, the above data is obviously quite different.

What is the best way to deal with this? I'm not a GIS expert - much more a coder.

Approaches that occur to me are:

Interpolate the points into a regular grid. Would gdal_grid as noted at https://gis.stackexchange.com/a/24610/58752 be workable?
Convert the points to triangles somehow and then interpolate the x,y from the triangle that it is in.

Crucially, I want to do this using command-line processing rather than point-and-click in GIS if possible, so that it is easily reproducible from the original data download.

Presumably the above screenshot indicates that the elevation field is called 'ALTURA'.

Any pointers would be useful.

(I'm aware that SRTM covers the area that this second data source does, but the data is a ground-based survey and avoids the problems of incorrect readings in SRTM due to building heights and trees, which are proving problematic.)

PS Is there a simple way to get the points out of the shapefile as a CSV? I can't seem to get this using either ogr2ogr or gdallocationinfo.


texts: How to set datum in GPS connection setup window in ArcFM Viewer for ArcGIS Engine Application

texts:
I have a custom application based on ArcFM Viewer for ArcGIS Engine. I want to use an external GPS receiver with that application and the setup window does not let me choose a different geographic datum. In comparison the corresponding drop down field in ArcGIS Desktop has a variety of datums to choose from. Where is this list defined? How can I get this in my ArcFM Viewer application? It might also be helpful to know where ArcGIS actually gets that list from. Here are some screenshots. Sorry for the German UI.

ArcFM Viewer:

ArcGIS Desktop:


texts:
I have a map of a city and a .csv-file containing 700 geographic coordinates pointing to locations within that city.

How can I plot these 700 locations on the map?


texts: Show images related to features in QGIS?

texts:
Say for example, if there's a point data with old trees, then when you use the identify tool, you can see the photograph of the tree

Or even one level further, the labels displaying these images right away without any clicking needed

I believe the term is file streaming?


texts:
I'm using PostGIS 2.5.2 on PostgreSQL 10, and I'm trying to build a function which would split a LineString into portions, depending on points positionned close to it. I'm getting a Geometry type (Point) does not match column type (LineString) error, and I can't figure out why. I've stumbled upon this old  answer but it's not helping me. Here's where the error comes from :

INSERT INTO portions_lignes(id_parcours, ordre_portion_parcours, id_arret, geom)
VALUES(line_id, prior_point.ordre_arret_parcours, prior_point.id_arret, ST_LineSubstring(linestring, 0, ST_LineLocatePoint(linestring, current_point.geom) / 2));

Where portions_lignes.geom is a LineString geometry type, linestring is also a LineString, and current_point.geom is a Point geometry type.

All geometry types are in SRID 4326, my goal here is just to insert into portions_lignes a particular line_substring (using ST_LineSubstring) from the linestring parameter, which is passed to the function and declared like this :

CREATE OR REPLACE FUNCTION splitLine(
    line_id character varying,
    linestring geometry(LineString,4326)
)


texts:
I'm trying to use the saga algorithm 'quadtree structure to polygon' (with points to get polygons), in the python console.
But I think I'm not able to set the parameters correctly.

When I call the algorithm help, I get :

----------------
Input parameters
----------------

SHAPES:  <QgsProcessingParameterFeatureSource>
    Shapes

ATTRIBUTE:  <QgsProcessingParameterField>
    Attribute

POLYGONS:  <QgsProcessingParameterVectorDestination>
    Polygons

LINES:  <QgsProcessingParameterVectorDestination>
    Lines

POINTS:  <QgsProcessingParameterVectorDestination>
    Duplicated Points

----------------
Outputs
----------------

POLYGONS:  <QgsProcessingOutputVectorLayer>
    Polygons

LINES:  <QgsProcessingOutputVectorLayer>
    Lines

POINTS:  <QgsProcessingOutputVectorLayer>
    Duplicated Points

Moreover, I get the error :

NameError: name 'QgsProcessingFeedback' is not defined

So I tried this...

from qgis.core import (QgsMessageLog,
                   QgsApplication,
                   QgsMapLayer,
                   QgsProcessingProvider,
                   QgsProcessingAlgorithm,
                   QgsProcessingException,
                   QgsProcessingParameterDefinition,
                   QgsProcessingOutputVectorLayer,
                   QgsProcessingOutputRasterLayer,
                   QgsProcessingOutputMapLayer,
                   QgsProcessingOutputMultipleLayers,
                   QgsProcessingFeedback)

... but still the same error.

How can I use this algorithm?


texts: Calculating area of polygon within square mile in ArcGIS?

texts:
I would like to create a raster, where in every cell it returns the area of land in a polygon feature.  So for example, every cell in the raster would give park land per square mile, or vacant land per square mile.

Edit: I'm having a hard time communicating what I want to do.  I have a polygon layer with park land.  I want to find, for each unit in a raster grid, the amount of park land within a square mile of each raster cell.


texts: GeoServer WFS fails to load in OpenLayers

texts:
I'm working with PostGIS->GeoServer->OpenLayers on localhost.  I followed a GeoServer tutorial to import some nyc sample data into my postgres/PostGIS database and link it to GeoServer.  The layer seems to display correctly in GeoServer using the OpenLayers layer preview.  I then attempted to use OpenLayers to display a WFS-T from GeoServer on localhost, routing the LH:8081 (where my GeoServer is located) through a proxy.cgi script to avoid Cross-Scripting. When I visit: proxy.cgi I do get the OpenLayers home page (without images).

But when I try to display the WFS vector layer in my OpenLayers HTML/JavaScript page in either Apache LH or Jetty LH:8081/geoserver/www/ the POST call to LH:8081/geoserver/wfs doesn't complete.  I am dumbfounded as to why this is happening.


texts: Matching feature classes by similar name for map automation

texts:
I was working on this a bit ago and had something which I thought worked however, it does not seem to be working at the moment. What I'm trying to do is match feature classes having similar files. One will be a buffered version of the original and have 'Buffer_' in the name the other will simply be the original feature class. I'd like to then bring the buffered version and original feature class into a .mxd, symbolize both, create a legend, save the .mxd and export the map to a pdf. Everything is working as it should with the exception of getting the buffered feature class to match it's original.

I am able to bring the buffered feature class into the .mxd, rename the .mxd to the name of the buffered feature class. I am also able to bring in an original feature class to the map document however, it is not the original feature class which corresponds to the buffered feature class.

The portion of my code which attempts to match the 2 feature classes looks like:

buffered = arcpy.ListFeatureClasses('Buffer_*')

original = [buffer.split("Buffer_")[-1] for buffer in buffered]

fccount = -1
buffered.sort()
for fc in buffered:
    fccount += 1
    print fc + " " + str(fccount)
    #REPLACE sandbox WITH tip_2014
    fcclean = fc.replace('tBuffer_','')
    print "fcclean = " + fcclean

    copymxd = templatemxd.saveACopy(r"map output directory\\" + fcclean + ".mxd")
    print "Saved Copy"
    mxd = arcpy.mapping.MapDocument(r"map output directory\\" + fcclean + ".mxd")
    print "MXD set"
    df = arcpy.mapping.ListDataFrames(mxd,"*")[0]

    buffLayer = arcpy.mapping.Layer(buffered[fccount])
    orgLayer = arcpy.mapping.Layer(original[fccount])

I have a feeling I am missing something with the count but am puzzled as to what. Thanks for any help you can provide with this issue.


texts: Joining Excel tables to data attribute in ArcMap?

texts:
I am trying to join an Excel table to my parcel attribute table in ArcMap.

I have run a validation and all checks out except it says there are no matches.

What would cause this?


texts: No ogr module despite having gdal installed?

texts:
I am trying to use the mapping component of the wonderful PySAL package, but for some reason I do not have the ogr module.  I was under the impression that it came packaged with gdal.  Am I missing something?

import gdal
import ogr

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-29-ef9eb3418665> in <module>()
----> 1 import gdal
      2 import ogr

ImportError: No module named gdal

UPDATE:  I am afraid I was in a bit of a hurry to catch a plane, so I failed to include some critical information.  I already installed python-gdal using the apt-get method described here, because I wanted to leave little to chance.  Is it possible that the package manager did not install it correctly?

choct155@choct155-Q550LF:~$ sudo apt-get install python-gdal
Reading package lists... Done
Building dependency tree
Reading state information... Done
python-gdal is already the newest version.
The following packages were automatically installed and are no longer required:
  firefox-globalmenu libfftw3-3 libgeos-3.3.3 libogdi3.2 libxerces-c28
Use 'apt-get autoremove' to remove them.
0 upgraded, 0 newly installed, 0 to remove and 33 not upgraded.

I am afraid I am not finding osgeo either...

from osgeo import gdal

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-28-a04b80bf4e65> in <module>()
----> 1 from osgeo import gdal

ImportError: No module named osgeo

Additionally, I should add that I am using Ubuntu 13.04.


texts: Qgis 2.10-kml is not a valid or recognized data source

texts: After upgrade from Qgis 2.8, importing the kml format doesn't work. The Kml file is not problem, because it works well on other computers. I tried to reinstall Qgis, but the problem is still here...Any ideas?


texts: What does it mean to access a WMS using a proxy?

texts: I am trying to access a WMS service and bring it in through ArcMap, and was given instructions on how to access it using a proxy. Using the proxy script I'm able to download a layer directly to my snapshot. However I don't really understand what it works, and how one gets from that step to using a proxy with ArcGIS so that ArcMap can bring in WMS layers. I have tried searching but I am not exactly understanding what the proxy does and how it is incorporated into a map service/online map.


texts:
I have a column with several data separated by a comma, I wanted to separate them into several columns and to rotate the X column created to have only one
N ° | circuit | Cx1 | Cx2
1   | C1, C2  | C1  | C2
2   | C1      | C1  |
3   | C1      | C1  |

Result hope
N ° | circuit | CP1 |
1   | C1, C2  | C1  |
1   | C1, C2  | C2  |
2   | C1      | C1  |
3   | C1      | C1  |


texts: What do the parameters for Google Streetview mean?

texts:
I've noticed that the new Google Streetview uses a url of the following structure:

https://www.google.com/maps/place/London/@51.508506,-0.125532,3a,90y,72h,90t/data=!3m4!1e1!3m2!1s-KZZeMFtZx0aq3L0MLqE-A!2e0!4m2!3m1!1s0x47d8a00baf21de75:0x52963a5addd52a99

Specifically, there is a string like this @51.508506,-0.125532,3a,90y,72h,90t

I understand lat/lon/heading (I think) but I'm unclear on zoom and tilt. I think it could be structured like this:

@51.508506 -> lon
-0.125532 -> lat
3a -> zoom?
90y -> tilt
72h -> heading (measured from North)
90t -> ?

Are there any references to these parameters, anywhere? I can't quite figure out the zoom parameter, as well as the 90t


texts: I'm working on QGIS 2.2, and I would like to make a vector layer extent bigger (thus the clip option is not suitable) to match exactly my mask layer extent. I would like the new added area to be empty. I can't seem to find a way to edit the extent.


texts:
I would like to have search control for Leaflet, with suggestions to street addresses, like on this image from Google Maps:

Some of the plugins I tried (https://labs.easyblog.it/maps/leaflet-search/) don't give such results, when I typed some street name in them they offered either less or no suggestions.
Here Google Maps responds properly for some random street name in my city.


texts: Easiest way to update a large vector dataset

texts:
About 6 months ago, I downloaded a county dataset of all of the parcels contained in that particular county, which was a zip file with about 78,000+ records. I added it to my QGIS as a vector file. I have a map with map views and print views using this dataset layer and all of the symbology set up exactly the way I like it. I called this layer "Land Use".
Today, I downloaded an updated file from the county, the new file has all the same attributes obviously, just with updated information. I added it to my map called, by default "parcels".
I was just going to delete the existing "Land Use" layer, and add the new "parcel" layer, but that would require me building all my definition queries, rules, and symbology all over again. I know there has to be a better way.
I tried opening the Attribute Table in "Land Use" and deleting all of the content, then opening the "parcels" table and copying all the data and then pasting it in the "Land Use" table since all attributes are the same. That worked for about the first 25,000 records (out of 78,000) which took about 15 minutes but then QGIS bogs down and stops responding, and my laptop gets hot and slow.
Surely there is a better way to update the data without adding a new layer, I just am not sure how to go about it. I am fairly inexperienced with QGIS/GIS in general.


texts: In QGIS 2.14.0 The measuring tool is set to Meters even after changing it to feet in settings

texts:
I need to measure some lengths in feet, in previous versions I could change the units in the drop down window and it would stay on feet now it reverts back to meters, if you take focus off the measure tool by for example selecting the next point to measure to.

I have gone to Settings >Options >Map Tools and change the distance units to feet but the tool still opens in meters


texts:
I have aerial imagery as RGB jpgs and I would like to be able to clip them to a specific extent but retain a world file so I can use them in other software.
I could do this by splitting them and then clipping and then joining them back together but that is not exactly a smart solution if doing it regularly.
I have considered using a custom sized composer but that is a bit rough.

This functionality is supported by by the clipper tool but for some reason the XYZ driver is chosen instead of JPEG or GTIFF - this can be corrected by editing the command manually.


texts: PHP-MBTiles Server implementation with Leaflet and Wax

texts:
Using PHP tile server from infostreams with Leaflet and Wax Connector for Leaflet to visualize my .mbtiles. I upped the source and the .mbtiles on a free hosting service.

I am getting a blank Leaflet container on my site right now. There are the zoom buttons and Leaflet attribution but no map.

When I point the browser to domain.com/tiles/tileserver.php it gives the open_basedir restriction error.
When I replace the .json URL (in java.js) with another map's .json that is uploaded on MapBox Hosting, that map shows (interactivity doesn't work though).

EDIT: Since I'm on a free hosting service, I cannot access the httpd.conf or php.ini file. Therefore I have no way of suppressing the open_basedir error. I could use another free hosting service of your choice.

Server directory looks like:

DIR tiles
DIR wax
.htaccess
index.html
java.js
css.css

index.html goes:

<!DOCTYPE html>

<link href='wax/ext/leaflet.css' rel='stylesheet' type='text/css' />
<link href='css.css' rel='stylesheet' type='text/css' />

java.js goes:

var url = "/tiles/Century.json"
wax.tilejson(url, function(tilejson) {

// Create map and add image tiles
var map = new L.Map('map-div')
    .addLayer(new wax.leaf.connector(tilejson))
    .setView(new L.LatLng(51, 0), 1);

// Create map legend
wax.leaf.legend(map, tilejson).appendTo(map._container);

// Add map interaction (tooltips)
wax.leaf.interaction()
    .map(map)
    .tilejson(tilejson)
    .on(wax.tooltip().animate(true).parent(map._container).events());

});

"tiles" directory:

.htaccess
Century.mbtiles   //This is just what I downloaded from infostreams' tile server.
tileserver.php


texts: Layer switcher triggered by HTML button in OpenLayers?

texts:
in this code you have an example of layer switcher use in OpenLayers

https://rawgit.com/walkermatt/ol-layerswitcher/master/examples/layerswitcher.html

The button which triggers the action is written in .js using ol.control.LayerSwitcher, though.

var layerSwitcher = new ol.control.LayerSwitcher({
    tipLabel: 'Légende' // Optional label for button
});
map.addControl(layerSwitcher);

What do I have to change if I want the same pop-up to be displayed, but by clicking the HTML icon in the "navbar" div of the site, not this button?


texts: Creating CoverageStore using REST API from GeoTIFF already in GeoServer DATA_DIR

texts:
I have a GeoTIFF already available in GeoServer DATA_DIR, specifically in a subdir named "rice" and I woluld like to publish it as a CoverageStore.

What I could understand is that REST API allows to create the coverage with a PUT request passing the file as the payload.
How can I skip this step?
The only documented option seem to allow for just the 'file[.extension]' option that requires passing in a file.
On the other side, datastores allow for an 'external' option that allows picking from DATA_DIR.


texts: Seeking introductory texts about coordinate system projections?

texts:
Even after some time working in this GIS field, I still don't have a solid understanding of spatial references, how they are projected, and most importantly, if you still need to worry about the fact that coordinates are geographic once they have been projected. In which cases should one pay attention when handling coordinates in a projected coordinate system?

For example, assume I am working in a projected coordinate system. I have a polygon whose vertices are now specified with X / Y coordinates. Let's say I want to manually calculate the polygon's area — could I now take any normal algorithm (for Cartesian/Euclidean coordinate spaces) for computing the polygon's area, or do I still have to account for the fact that these are geographic coordinates? What would I have to do to get an accurate result for my area computation?

It's questions such as the above example that I still don't understand. Any cues to introductory texts explaining such details are welcome.

(I was already considering looking at the PROJ.4 catrographic projection library and how it is used to deduce from that how to deal with projected coordinates, but that would've probably been too early.)


texts: How to check if a geometry of type Point is inside another geometry of type MultiPolygonZM using Postgis functions?

texts:
I want to know if there's a way to check if a geometry of type Point is inside a geometry of type MultiPolygonZM, ignoring the Z coordinate, I just want to know if my Point(lat, lng) is inside a MultiPolygonZM.

This is an example of a MultiPolygonZM geometry in my database; I used ST_AsText to get a text representation of the geom :

"MULTIPOLYGON ZM (((
570277.5312 3592073.8077 460.900200000004 -1.79769313486232e+308,
570262.0283 3592073.7806 460.712700000004 -1.79769313486232e+308,
570262.4922 3592086.4664 460.705199999997 -1.79769313486232e+308,
570278.0406 3592086.5792 460.697700000004  (...)"

I tried ST_CONTAINS, ST_INTERSECTS like this :

SELECT * FROM CITIES
WHERE
ST_Intersects(ST_SetSRID(geom, 4326), ST_SetSRID(ST_MakePoint(3.67622580, 32.48549230), 4326))

SELECT * FROM CITIES
WHERE
ST_Contains(ST_SetSRID(geom, 4326), ST_SetSRID(ST_MakePoint(3.67622580, 32.48549230), 4326)).

My table uses 4326 as SRID, but i wanted to set it using ST_SetSRID to be sure that all the geoms have the same SRID

But it doesn't work! All the results are false!

Any solution ?


texts: How to convert EPSG:2163 coordinates to wgs84 in javascript

texts:
I'm currently testing a web app that calculates EPSG:2163 coordinates based on user interaction with a map created through d3.js. I wish to convert these coordinates to WGS84 in order to use various APIs whose parameters are in the form of WGS84 coordinates. Is there an simple formula or javascript library that can do this?

To further clarify, the map seen http://spatialreference.org/ref/epsg/2163/ takes WGS84 coordinates as input and converts them to EPSG:2163. All I want to do is reverse that calculation.

EDIT

Here's an example, however, I'm not convinced I'm getting correct results.  I started with WGS84 coordinates and projected them to epsg:2163 using cs2cs.  In the Javascript I transform the epsg:2163 coordinates back to WGS84, but as you will see, they don't end up the same as my starting coordinates.  I'm a little puzzled, myself, and I sure could have botched a step:

<input type="button" value="Project" onclick="javascript:project();" />


texts: Clipping Raster Object with Polygon Object in ArcPy?

texts:
I have about 400 rasters that I want to clip each of them based on a polygon object.

The code I tried to use is the following:

import arcpy, os
origpath = ur'd:\test'
jpgs = [f for f in os.listdir(origpath) if f.endswith('.jpg') and os.path.isfile(os.path.join(origpath,f))]
for jpg in jpgs:
    r = arcpy.sa.Raster(os.path.join(origpath,jpg))
    p = arcpy.Polygon(arcpy.Array([arcpy.Point(562879.802,4538631.518),arcpy.Point(583906.505,4538842.192),arcpy.Point(584139.129,4518489.715),arcpy.Point(563054.256,4518279.244)]))
    arcpy.Clip_management(r,p,os.path.join(origpath,'clipped_'+jpg),'#','#','ClippingGeometry')

When it runs it produces an error:

ERROR 000628: Cannot set input into parameter rectangle.

I assume that it has to do with the Polygon Object, I have inserted, for the Extent parameter.

The last solution would be to put this one polygon in a feature class and then add this feature class as the extent parameter.

Is there any direct way to use the Polygon Object to clip my rasters?


texts: Quantum GIS Reproject a projected coordinate system with rotation

texts:
I have two shapefiles:- one in WGS84 Z30S and one in a custom, system based, projection. I need to convert the custom projection into WGS84 Z30S.

I have a set of parameters including a shift in X and Y rotation, scale, and origins for rotation and scale.

How do I format a custom coordinate system in QGIS to achieve this? Is there a better way of doing it?

The syntax is foxing me at the moment.


texts: What's the correct transformation for Lambert Conformal Conic for Wyoming?

texts:
I am a newb in QGIS. I have loaded QGIS 1.8 (Lisboa) and practicing with a number of maps for portions of the state of Wyoming. My maps are mostly in UTM 12N NAD 27 & NAD 83 some are in decimal degrees and one is Lambert Conformal Conic NAD 83 with false Easting and False Northing offsets.

There is no problem handling the maps in UTM 12N NAD 83 to overlay on those in UTM 12 NAD 27. But I cannot find a CSR Transformation for the Lambert Conformal Conic for Wyoming. The list contains practically all the states but not Wyoming.

Am I blind, or am looking in the wrong place, or is it simply not there? Help would be greatly appreciated. Thank you.


texts: Spatial join not assigning IDs correctly in QGIS

texts:
I'm using some of the spatial join functions (join attributes by location) to assign points to a polygon grid (that was created in QGIS using spatial tools --> create grid), but it's assigning different random IDs to points clearly in the same grid. I've tried within, contains, etc and all have the same issue. Yes I've double checked projections.

Any idea what's going on?

In QGIS v 3.8

Solution (from @Leo M): this occurred when working with the grid layer as a temporary layer. I went in and saved the grid layer as a .shp and it's now joining by location accurately.


texts: World map with clear difference between land, sea and inland waterways

texts:
I'm looking for a free high resolution map that clearly indicates the land and sea in the world. I need it to calculate distances over sea between ports.
Additionally it should also include inland waterways and canals.

Does anybody know whether this data is freely available?


texts:
I have file with coordinates X, Y of point objects in EPSG3301 coordinate system (that means it is in meters).

I want to spatially cluster points. Here is my R code (dbscan package):

db2 <-dbscan(cbind(point2$X, point$Y), eps = 500, MinPts = 25)

In what units is eps in this case? How can I achieve eps to be axactly 500 meters?


texts:
I am working on Geo-server 2.16.x on Cent-OS & I have configured "auth-key" extension.
While configuring extension I used /**  "ANT patterns".
I don't know the meaning of ANT patterns, I simply copied default and re-named as 'api-key'
Could anyone tell me what is 'ANT patterns' ? And why they are required in Geo-server ?

My primary research says that it's related to Java Spring framework and related to security but it's great if anyone explains in simple language.


texts: Display final desired color with RGB + transparent grey hillshade

texts:
In QGIS 2.18.16, I am trying to end up with a final desired color specified in RGB (from a style branding book of a client).  The map also contains a hillshade with some transparent grey over top of the polygon colors.  When I specify the polygon colors in the RGB I need, predictably the final colors are actually the sum of the RGB value plus the effect of the partially transparent grey hillshade.

Is there a way to correct for this in a predictable way?  All I can come up with so far is to trial-and-error adjust the colors via RGB in the presence of the transparent hillshade to get closer to my desired RGB from the style book.


texts: Creating virtual raster produces gaps along edges in QGIS?

texts:
I've got almost 1500 .adf raster files of LiDAR data that I need to merge/mosaic into one raster. The most efficient way that I found was to create a virtual raster of all of the files:

gdalbuildvrt -input_file_list out.txt LiDAR_virtual.vrt

and then convert the virtual raster into a GeoTiff:

gdal_translate -of GTiff -co "TILED=YES" LiDAR_virtual.vrt LiDAR_mosaic.tif

However, after the first command, the virtual raster that is created has some gap-lines that become more apparent at finer scales and also get translated to the .tif image.


texts: I am new at QGIS and GIS in general. Basically, I have a shapefile of continents. However, I need to include further refinements. For example, within "Oceania", I need to identify New Zealand and New Guinea. How can I create a new feature that identifies "New Zeland" while it is no longer one of the polygons identified as "Oceania"?


texts: Sharing a Folium Interactive Map w/ Embedded Markers

texts:
I have been experimenting with a few Python packages that have great potential for my dissertation research.  Since I use pandas for pretty much everything, geopandas was a natural choice for working with spatial data (within the range of capability of the package), and the interactive properties of folium are ... sweet.

In any case, I want to be able to share the maps that I have made with other folks who are not using or even familiar with IPython Notebook.  My test case is plotting the location of traffic cameras in DC.  (The shapefile is provided to the public by DCGIS via opendatadc.org, a new and growing source for District info.)  The Notebook Viewer seemed like an obvious choice, but as can be seen from my rendered gist, the maps are 1) static, and 2) displayed only in their final form.  (Markers were not added until the last map, even though they appear in the first one.)

I then thought I could leverage Mike Bostock's awesome viewer for interactive graphics, bl.ocks.org.  I am missing something, however, because while my map is now interactive, the markers have disappeared.

Has anyone run into this problem?  Any help would certainly be appreciated.


texts: Creating SqlView Layer with GeoServer Manager?

texts:
I would like to register a new SQL view layer to GeoServer through the REST API, I am using geoserver-manager.

The registration is OK, the layer is working fine, except the WMS GetFeatureInfo, which returns:

"Operation on mixed SRID geometries".

If I manually set the SQL view's geometry column's SRID in the GeoServer administration page, it works fine.

But I cannot find any way to set it in the XML request. I'm sending the geometry info:

<geometry>
    <name>geom</name>
    type>Point</type
    <srid>4326</srid>
</geometry>

And the geometry attribute:

<attribute>
    <name>geom</name>
    <minOccurs>0</minOccurs>
    <maxOccurs>1</maxOccurs>
    <nillable>true</nillable>
    <binding>com.vividsolutions.jts.geom.Point</binding>
</attribute>

But I cannot set the attributes SRID and the geometry info is not enough.
The geometry attribute will be registered as SRID -1.

Is there an API call to refresh the layer fields?


texts: Working with ESRI WKT Spatial Reference in ArcGIS API for JavaScript?

texts:
I have been practising ArcGIS in JavaScript (for WMS layers ) and have always faced numerical values of SpatialReference like wkid: 102100, but now I faced in this format:

'PROJCS["NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601",' +
'GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",' +
'SPHEROID["GRS_1980",6378137.0,298.257222101]]' +
'PRIMEM["Greenwich",0.0],' +
'UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],' +
'PARAMETER["False_Easting",8202099.737532808],PARAMETER["False_Northing",0.0],' +
'PARAMETER["Central_Meridian",-120.5],' +
'PARAMETER["Standard_Parallel_1",44.33333333333334],' +
'PARAMETER["Standard_Parallel_2",46.0],' +
'PARAMETER["Latitude_Of_Origin",43.66666666666666],UNIT["Foot",0.3048]]'

I came across  the explanation of the wkt usage but do not know how to use it together with extent. Is there any correct syntax / solution / tutorial I could learn from?


texts: Using Python or ArcPy to extract usable shapefile data from several hundred MXD files

texts:
I have several hundred (potentially more) old MXD files which were made to show grass areas to be cut along city roads. The employee before me opted to do these one by one and rather than create actual data, used the ArcMap drawing tool to draw the shapes over ortho imagery in the ArcMap layout view.
I have a current request to provide this data on City-wide maps to management and am entirely baffled as to a solution being possible beyond having to manually redraw these and loading each file to do so.
How can I tackle this problem?


texts: I created an image in QGIS using the Google Maps layer and I would now like to publish this as part of an academic journal article. However, the journal editor has asked if I have permission to use the Google Maps base layer. Does anyone know if I need to obtain permission from Google to use this image?


texts: issues in table reading from MS Access using ODBC

texts:
I want to add two tables from my MS Access db.
So I follow the procedure to add ODBC (add vector layer, database, etc.) and it works.
I choose the tables that I want to add and hit ok.
Qgis adds both the tables I chose, BUT are showed in a different way.
In the first table I’m not able to see any data, just the name of the fields.
In the second one I see everything.

Somebody can explain to me if I’m doing something wrong?


texts: Is there a limit to have many layers you can have in a leaflet map?


texts: US Census roads - names of overlapping street segments

texts:
I'm working with recent Census (TIGER) road files in ArcGIS 10.5 / ArcPy (Advanced license).  I want to get the road name(s) of the closest road segment (between intersections) to a point. Roads in this dataset frequently have different names which are represented by overlapping lines. The example below is output from Intersect with line output; this type of overlap occurs for a good portion of my dataset.

To get the roads into segments broken by the intersections, I used Unsplit Line then Feature to Line.

Using Near analysis would get me the closest road to specific points, but doesn't reflect the other road names / aliases. I'm looking for something similar to Spatial Join where you can have a combined output field, but Spatial Join requires two input layers. What I've thought to try next are either a) copy the layer and use Spatial Join or b) Use the Near Tool with a maximum of 10 nearest objects and somehow query the table (cursor?) based on NEAR_FID and NEAR_RANK and write the results to a new field. Am I missing any better approaches here?

This question is similar to how to identify cross streets for points that are not at intersections? but I'm looking for all of the names of the road the point is on, not the names of the cross streets.


texts: Adding graduated symbology to created arrows using PyQGIS

texts:
I would like to add graduated Symbology to my custom created arrows. Unfortunately it does not seem to work. When manually creating graduated symbology the steps are rather easy.
Could someone help me with this last step?

And this is what I would like to see eventually:

And my code I created so far:
from qgis.core import (
    QgsVectorLayer,
    QgsSingleSymbolRenderer,
    QgsLineSymbol,
    QgsFillSymbol,
    QgsProject,
    QgsRectangle,
    QgsArrowSymbolLayer,
    QgsMultiRenderChecker,
    QgsProperty,
    QgsSymbolLayer,
    QgsMapSettings,
    QgsSymbol
)

layer = QgsVectorLayer( "?query=SELECT from_country, ubc_processing_plant, delivered_tonnes_2019, make_line(Centroids_countries.geometry, UBC_processing_2030.geometry) FROM ubc_flows JOIN Centroids_countries ON ubc_flows.from_country = Centroids_countries.NUTS_NAME JOIN UBC_processing_2030 ON ubc_flows.ubc_processing_plant = UBC_processing_2030.name", "vlayer", "virtual" )
QgsProject.instance().addMapLayer(layer)

def test_3(self):
        sym = self.renderer().symbol()

        shape_sym = QgsGeometryGeneratorSymbolLayer.create({'geometryModifier': 'difference(difference(make_line(start_point($geometry), centroid(offset_curve($geometry, length($geometry)/-10.0)), end_point($geometry)), buffer(start_point($geometry), 10000)), buffer(end_point( $geometry), 10000))', 'outline_color': 'black', 'symbolType': '1'})
        shape_sym.setSymbolType(QgsSymbol.Line)

        arrow = QgsArrowSymbolLayer.create(
            {
                "arrow_width": "3",
                "head_length": "3",
                "head_thickness": "3",
                "head_type": "0",
                "arrow_type": "0",
                "is_curved": "1",
                "arrow_start_width": "3"
            }
        )

        shape_sym.subSymbol().changeSymbolLayer(0, arrow)
        sym.changeSymbolLayer(0, shape_sym)

        vals = []
        fld = 'delivered_tonnes_2019'
        for f in self.getFeatures():
            vals.append(int(f[fld]))

        colors = ['#fffff4', '#fcfdce', '#ffffbf', '#abdda4', '#2b8300']

        vals = list(filter(None, vals))
        vals = [float(s) for s in vals]
        lower = float((sorted(vals))[0])
        upper = float((sorted(vals))[-1])
        step = (upper-lower)/len(colors)
        range_list = []
        for c in colors:
            cat = [lower, lower+step, c]
            sym = QgsSymbol.defaultSymbol(self.geometryType())
            sym.setColor(QColor(cat[2]))
            rng = QgsRendererRange(cat[0], cat[1], sym, '{0:.1f}-{1:.1f}'.format(cat[0], cat[1]))
            range_list.append(rng)
            lower = (lower+step)

        renderer = QgsGraduatedSymbolRenderer(fld, range_list)
        # Convert to rule based renderer ("Single Symbol" to "Rule-based")
        # If not required, just comment below two lines and keep the last line
        newRenderer = QgsGraduatedSymbolRenderer.convertFromRenderer(self.renderer())
        self.setRenderer(newRenderer)
        self.triggerRepaint()

test_3(layer)


texts:
I'm trying to create what I call a "many to one" join. I don't know if that would be the correct term or not. I have a table that has unique account numbers for mobile homes (i.e. - M1007970) per parcel account number (R0003285). (Many mobile homes per parcel - many to one.) I need to join this table to our parcel geometry - and still only have one polygon per parcel.

So, for example, the table may have three rows that have mobile home account numbers M1007370 on one row, M1007371 on another, and another one with M1059370, but all have the same parcel number R0032585. Our parcel geometry would only have the same field of R0032585.

When joining I have 12,088 mobile home records and 44,103 parcels. If I "keep all records" I have 44,103 records with only 7,947 mobile home account numbers (of the original 12,088). If I join based on "keep only matching records" I end up with just 7,947 records total.

I have done it successfully in the past and created a model. In this model I use the table for the mobile homes to join to the parcel layer (.lyr - the only way you can/could join in a model) based on the parcel account number. I copy the features keeping only matching records to a file geodatabase. From the file geodatabase I then append it into our SDE system. This currently quit working for reasons that I can't fathom, as nothing has changed.

Perhaps someone can convey better than I can what I'm trying to do, and if it's called something other than a many to one relationship (I don't believe it's a one to many...).


texts:
When I try to create chart by bands and image spectra i get "Error generating chart: Data column(s) for axis #0 cannot be of type string"
// Plot band values at points in an image.
var L8_SR = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

var COLOR = {
  built_up: 'ff0000',
  grassland: '00ff00',
  forest: '228B22',
  water: '0000ff',
  snow: 'add8e6',
  rocky: '808080',
  cultivation:'ffff00'
};

// 7 known locations.
var built_up = ee.Feature(
    ee.Geometry.Point(83.978710, 28.3202722), {'label': 'built_up'});
var grassland = ee.Feature(
    ee.Geometry.Point(84.008199, 28.218142), {'label': 'grassland'});
var forest = ee.Feature(
    ee.Geometry.Point(84.2250, 28.371391), {'label': 'forest'});
var water = ee.Feature(
    ee.Geometry.Point(83.947601, 28.215001), {'label': 'water'});
var snow = ee.Feature(
    ee.Geometry.Point(83.895411, 28.588669), {'label': 'snow'});
var rocky = ee.Feature(
    ee.Geometry.Point(84.042354, 28.523533), {'label': 'rocky/Barren'});
var cultivation = ee.Feature(
    ee.Geometry.Point(84.085761,28.170022), {'label': 'cultivatiom'});

var kaskiPoints = ee.FeatureCollection([built_up, grassland, forest,water,snow,rocky,cultivation]);
var landsat8SR = L8_SR.filterBounds(kaskiPoints);

var KaskiImage = ee.Image(landsat8SR.first());

// Select bands B1 to B7.
KaskiImage = KaskiImage.select(['B[1-7]']);

var bandChart = ui.Chart.image.regions({
  image: KaskiImage,
  regions: kaskiPoints,
  scale: 30,
  seriesProperty: 'label'
});
bandChart.setChartType('LineChart');
bandChart.setOptions({
  title: 'Landsat 8 band values at 7 classes in kaski district',
  hAxis: {
    title: 'Band'
  },
  vAxis: {
    title: 'Reflectance'
  },
  lineWidth: 1,
  pointSize: 4,
  series: {
    0: {color: COLOR.built_up},
    1: {color: COLOR.grassland},
    2: {color: COLOR.forest},
    3: {color: COLOR.water},
    4: {color: COLOR.snow},
    5: {color: COLOR.rocky},
    6: {color: COLOR.cultivation}
  }
});

// From: https://landsat.usgs.gov/what-are-best-spectral-bands-use-my-study
var wavelengths = [.44, .48, .56, .65, .86, 1.61, 2.2];

var spectraChart = ui.Chart.image.regions({
  image: KaskiImage,
  regions: kaskiPoints,
  scale: 30,
  seriesProperty: 'label',
  xLabels: wavelengths
});
spectraChart.setChartType('LineChart');
spectraChart.setOptions({
  title: 'Landsat 8 Surfacereflectance spectra at 7 classes in kaski district',
  hAxis: {
    title: 'Wavelength (micrometers)'
  },
  vAxis: {
    title: 'Reflectance'
  },
  lineWidth: 1,
  pointSize: 4,
  series: {
    0: {color: COLOR.built_up},
    1: {color: COLOR.grassland},
    2: {color: COLOR.forest},
    3: {color: COLOR.water},
    4: {color: COLOR.snow},
    5: {color: COLOR.rocky},
    6: {color: COLOR.cultivation}
  }
});

print(bandChart);
print(spectraChart);

code link : https://code.earthengine.google.com/c24c20cbf4b0b9a2e820e541f82acc6b


texts: Why are some features missing from Raster menu in QGIS for linux?

texts:
I am running QGIS 1.8 in Ubuntu/Linux (13.04), with GDAL 1.10 installed.  There seem to be features missing from GDAL tools that are present in Windows version. For example, I can't find Rasterize or Contour options in the menue, and Clipper only has "extent" option - no "mask to layer", which is what I need!

Am I missing something?  How do I get these additional functions?

Thanks


texts: Removing/excluding items from legend in an ArcFM Silverlight viewer

texts:
I'm customizing an ArcFM Silverlight viewer (which is essentially an ArcGIS Silverlight viewer in a different package and with some additional tools). I managed to make a legend via the Legend Class and Sample: Legend with Templates.

The last three layers (InsetMapMarker, TextTemporary, and BufferLayer) are things automatically generated by the ArcFM viewer tools. It will confuse to end users to see them ("what's a buffer layer, I'm not making any buffers"), and potentially mess with functionality if they turn them on or off. So I really ought to exclude them from the legend.

I was very happy to find this Q&A How do I exclude raster layers from legend? (Silverlight, ArcGIS Server) and adapted the code snippet to exclude the last three layers (see below). However, this has no effect -- they're still visible in the legend.

private void Legend_Refreshed(object sender, Legend.RefreshedEventArgs e)
{
    LayerItemViewModel removeLayerItemVM = null;

    if (e.LayerItem.LayerItems != null)
    {
        foreach (LayerItemViewModel layerItemVM in e.LayerItem.LayerItems)
        {
            if (layerItemVM.IsExpanded)
                layerItemVM.IsExpanded = false;
            if (layerItemVM.Label == "InsetMapMarker")
                removeLayerItemVM = layerItemVM;
            if (layerItemVM.Label == "TextTemporary")
                removeLayerItemVM = layerItemVM;
            if (layerItemVM.Label == "BufferLayer")
                removeLayerItemVM = layerItemVM;
        }

        if (removeLayerItemVM != null)
            e.LayerItem.LayerItems.Remove(removeLayerItemVM);
    }
    else
    {
        e.LayerItem.IsExpanded = false;
    }
}

I know very little about C# and the Silverlight API, and guessing isn't getting me very far. What am I missing?


texts:
How to update an .img DEM in ArcView 9.3?

I do have DEM which need to be updated only if few parts. The area to be updated has been produced in separate DEM with Spline Interpolation with Barriers.

The area to be updated itself is S-shape, but I received larger- rectangular DEM (btw. I was hoping to receive -9999 value outside barriers, but it seems that outside area was also interpolated).

My aim is to to update only area within the polyline/barriers. How to do that?


texts: How to get 4 values for different Indexes and names of rasters for set of points?

texts:
I am trying to  get 4 values for different Indexes and names of rasters  for set of points, but get problem with defining function    def extractIndex(raster, pathName).

#List of Folders  4 Folders     for   rasters with  different Indexes:
path1 = r"C:\Users\1_fire_dNBR/"
path2 = r"C:\Users\1_fire_RdNBR/"
path3 = r"C:\Users\2_fire_dNBR/"
path4 = r"C:\Users\2_fire_RdNBR/"

#Points:
inputPoints = r"C:\Users\select3.shp"

# Rasters in specific folder:
RExtents = [arcpy.Describe(R).extent for R in arcpy.ListRasters()]

#Names of rasters:
 RNames = [arcpy.Describe(R).name[:-9] for R in arcpy.ListRasters()]  # extract  names of rasters

def extractIndex(raster, pnt, pathName):
    cr = arcpy.ListRasters(RNames[counter]+'*')  # names of fires
    rast = arcpy.Raster(cr[0])  #Raster Object   (path and raster   here)
    V = arcpy.RasterToNumPyArray[rast, pnt, 1, 1]   # specific value of Point from   raster
    return V[0]
#Code:
with arcpy.da.UpdateCursor(inputPoints, ["SHAPE@","dNBR_1", "RdNBR_1", "dNBR_2", "RdNBR_2", "name"]) as iCur:
    for pnt in iCur:
       counter = 0
       for iR in RExtents:
            if iR.contains(pnt[0]):

                # Names of path  with  Raster's Data Sets:
                dNBR_1= extractIndex(RNames[counter], pnt, path1)
                RdNBR_1= extractIndex(RNames[counter], pnt, path2)
                dNBR_2= extractIndex(RNames[counter], pnt, path3)
                RdNBR_2= extractIndex(RNames[counter], pnt, path4)

                iCur.UpdateRow[pnt[0],dNBR_1, RdNBR_1, dNBR_2, RdNBR_2, RNames[counter]]
              counter =+ 1

Running the above code gives an error:

  File "", line 33, in
      dNBR_1 = extractIndex(RNames[counter], pnt[0], path1)

  File "", line 20, in extractIndex
      V = arcpy.RasterToNumPyArray[rast, pnt, 1, 1]   # specific value of Point from specific  raster

  TypeError: 'function' object has no attribute 'getitem'


texts:
Is there a limit of how many chars that can be used in a GeoNetwork search?

We have a card called "Grönområden inom 200m för förskolor 300m för skolor ÖP 2030 rev 2018", 69 characters long. It is possible to write that in the search field but it won't display as a result.

Is there a limit of 50 characters or so?


texts: Arcpad Message Box and Keyboard Button Issue

texts:
I'm using ArcPad 8 along with ArcPad Studio to create a custom application for use on Windows Mobile 6 devices. There are data entry forms and on those forms are buttons that, when clicked, present the user some information in a message box. For a button's onclick event, I specify an event script (using JScript) such as:

Application.MessageBox("Some information for the user...");

That all works fine. However, when the user dismisses the message box, the button for toggling the input keyboard/transcriber/recognizer that is located in the lower right corner of the application window immediately disappears, which prevents the user from showing and hiding the keyboard as they attempt to enter data in other various controls on the form. Has anyone else ever had this problem and know of any solutions? Thanks!


texts: GRASS 7.8.3 problem with r.watershed and r.water.outlet

texts:
I need to find the basin of a small river and this is what I did:

I found the drainage direction with r.watershed
I used r.water.outlet to find the basin for the coordinates of the water outlet point.

The problem is that this works for some bigger rivers but not for the one I need.
I saw that in the output "stream" of r.watershed my river doesn't appear, could be it the problem? What can I do? Here the images of "accumulation" form r.watershe in which you can see the river and the images of "stream" where there river doesn't appear.


texts: Shapefile for same region but they are not overlapping

texts:
i have two shapefile for same region but they are not overlapping with each other. one data set is admin boundaries which i have purchased form the survey of India and other data set is ward boundaries of the cities which i have downloaded form the government website. i have projected the ward boundaries in same projection system as the admin boundaries but still its not overlapping. what i should do.??

Admin boundaries details:
Projected Coordinate System:    LCC_WGS84
Projection: Lambert_Conformal_Conic
False_Easting:  4000000.00000000
False_Northing: 4000000.00000000
Central_Meridian:   80.00000000
Standard_Parallel_1:    12.47294400
Standard_Parallel_2:    35.17280600
Scale_Factor:   1.00000000
Latitude_Of_Origin: 24.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_WGS_1984
Datum:  D_WGS_1984
Prime Meridian:     Greenwich
Angular Unit:   Degree


texts: Making New Feature Class with Midpoint of Line Feature using ArcGIS Desktop?

texts:
I have a problem with creating a new shapefile ( for Points) with a midpoints from Line shapefile.

I cannot use FeatureVerticesTo Points_management because I only have access to Basic and Standard level licenses of ArcGIS Desktop 10.5.1.

Is there another way to do this using ArcGIS Desktop tools or ArcPy functions?


texts: Convert NAD_1983_Georgia_Statewide_Lambert to NAD83

texts:
All I have a ESRI generated shapefile that uses: NAD_1983_Georgia_Statewide_Lambert

Here is the projection information:

PROJCS["NAD_1983_Georgia_Statewide_Lambert",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",656166.67],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-83.5],PARAMETER["Standard_Parallel_1",31.41666666666667],PARAMETER["Standard_Parallel_2",34.28333333333333],PARAMETER["Latitude_Of_Origin",30.25],UNIT["Foot_US",0.3048006096012192]]

I believe the EPSG code is 102604

Can I use GDAL To re-project this to either a NAD83 or WGS84 projection?
is there a better way to re-project this using a free tool? I'm having trouble verifying that ogr2ogr support this since it isn't a EPSG assigned range.


texts: How to use a clipped Raster with Raster Calculator in a script

texts:
I want to make a script which will asking me a raster layer, a vector layer and a double.
The point is to clip the raster and then process to filter data, like (raster > filter)*raster.
My script without the init algorithm:
layer = self.parameterAsRasterLayer(parameters, self.INPUT, context)
clip = self.parameterAsVectorLayer(parameters, self.MASK, context)
filter = self.parameterAsDouble(parameters, self.FILTER,context)
cliplayer = self.parameterAsOutputLayer(parameters, self.OUTPUT_CLIP, context)
calclayer = self.parameterAsOutputLayer(parameters, self.OUTPUT_CALC, context)

cliplayer = processing.run("saga:cliprasterwithpolygon",
            {'INPUT':layer,
             'POLYGONS':clip,
             'OUTPUT': cliplayer})

processing.run("gdal:rastercalculator",
            {'INPUT_A' : cliplayer,
             'BAND_A' : 1,
             'FORMULA' : '(A >' + str(filter) + ') * A',
             'OUTPUT' : calclayer})

I have an error saying the gdal:rastercalculator is invalid. But I don't know how to code it correctly, if someone could help me ! :)

EDIT
I find it ! I may call it back with ['OUTPUT'] as it is temporaly layers, here is my code if someone have the same problem :

cliplayer = processing.run("saga:cliprasterwithpolygon",
            {'INPUT':layer,
             'POLYGONS':clip,
             'OUTPUT': 'memory:' })

processing.run("gdal:rastercalculator",
            {'INPUT_A' : cliplayer['OUTPUT'],
             'BAND_A' : 1,
             'FORMULA' : '(A >' + str(filter) + ') * A',
             'OUTPUT' : calclayer})


texts: How do I not specify a CustomerWeightField when scripting CustomerDerivedTA_ba?

texts:
When you run the Business Analysis toolbox wizard for creating a Customer Derived Trade area, you have the option of weighting each customer by count or for specifying a field (such as sales or number of visits) to use as a weight. You select this with radio buttons and if you select "count", then you don't have to specify a field.

However, when scripting this tool, the CustomerWeightField is a required field. Leaving it blank yields a "This field was not found in the customer layer" type error message. In my case, I just want to use the count as the customers I'm analyzing do not have any other valid fields for weighting the results.

Is there a magic value to enter for CustomerWeightField that will force it to use count?

The documentation says:

  The field used to calculate the trade areas. This is based on either
  the number of customers (count) or the calculated weighted value
  assigned to each customer.

A fallback, I guess, would be to add a weight field and set it to 1, but I don't want to have to rely on modifying the input data.


texts: Sentinel-2 image size in Google Earth Engine

texts:
I am interested in the spatial coverage of Sentinel-2 and noticed that for an image of size 10980x10980 px (á 10*10 m²), the computed area covered is not 12056.04 km² as expected, but 12052.267583 km². I also tried the establish a band with the resolution in m² and got 12043.810046 km².

This is the code:

// Load single image
var img = ee.Image("COPERNICUS/S2/20180613T000139_20180613T000140_T09XWL")

// Define native pixelArea band (100 m**2 for S-2 B2)
var nativePixelArea = img.select('B2')
                         .multiply(0)
                         .add(100)
                         .rename('nativePixelArea');
// Compute pixel area as new band
var computedPixelArea = img.select('B2')
                           .multiply(0)
                           .add(1)
                           .multiply(ee.Image.pixelArea())
                           .rename('computedPixelArea');

// Add bands to image.
img = ee.Image(img).addBands([nativePixelArea, computedPixelArea]);

print(img.select('computedPixelArea').reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: img.geometry(),
  maxPixels: 1e12,
  scale: 10,
  }));

print(img.select('nativePixelArea').reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: img.geometry(),
  maxPixels: 1e12,
  scale: 10,
  }));


texts: Alternative backup File or Personal Geodatabase

texts:
I'm wondering if there is a way around or an alternative way to backup a File or Personal Geodatabase.

For example:

If I have a very large dataset within the File or Personal Geodatabase and I want to backup it up to several media storage devices such as External Hard drive, flash drive, DVD etc. Can this process be automated instead of performing it manually.

I work under the Government, so I may be limited to what I can do here in the office.


texts: Creating complementary graph from GIS data (e.g., shapefiles)

texts:
I am new to GIS.
I have some data about shapes, e.g., voting-districts, in several formats such as WMS, JSON, SHP. I want to create a graph (e.g., graphML) such that each district corresponds to a node, and edges indicate that two districts share a border. The nodes should be located in the "mass-point" of a district and edges should be associated with the distance between the nodes. Something like the red graph below.

Nodes should further store some properties, e.g., longitude, latitude, custom-data.
What are efficient ways to achieve this?
Let me know the pros and cons of your suggestions.


texts: Set zoom level to display full extent of specific layer when starting a project in QGIS

texts:
I have project in QGIS 3.2 with different shapefiles and WMS layer. When I start the projectfile (.qgs) I want the extent of the map to be the same as one specific WMS layer (same as the option "zoom to layer" in QGIS).

I created the following macro, where I get the right scale of the layer but the extent is not adjusted to the specific layer.

def openProject():
    from qgis.utils import iface
    iface.mapCanvas().zoomScale(4096000)

Is there an macro or even an option in QGIS to accomplish this?


texts: How to improve the precision of KML files on Google Maps?

texts: I am adding kml file to the google maps (api v3). Now the problem is that on google earth and google maps the accuracy of my kml decreases to kilometers. I am designing constituency boundaries and i need that to be accurate. Files are properly rectified, projected to WGS1984. In ArcGIS they give proper locations but when i plot their kml, it somehow is not that accurate. In addition, I am unable to go with point data in kml files. I mean they do not show the points neither go to its location.


texts: gdal/python WriteArray working for 8-bit but not 16-bit

texts:
I'm trying to use GDAL/python API to open a raster, do some math, then write the two resulting arrays as rastes. One is to be an 8-bit raster and the other 16-bit signed. However, only the 8-bit raster is coming out correctly. The 16-bit raster comes back as a raster of 0's. Here is the relevant code:

import numpy as np
import gdal
gdal.UseExceptions() # enable exceptions to report errors
drvtif = gdal.GetDriverByName("GTiff")

# open raster, read as array, do some math, return two arrays: arr1 and arr2

# arr1 and arr2 are both float32, though not necessary

print np.unique(arr1)
# prints [-15000, -10000, ..., 10000]
print np.unique(arr2)
# prints [1, 2, ..., 11]

# now write arr1 and arr2 as .tif:

outNDVI = '/path/to/outputs/output1.tif':
ndvidrv = drvtif.Create(outNDVI, ncols, nrows, 1, 3) # 3 = 16 bit signed
ndvidrv.SetGeoTransform(hc_gt)
ndvidrv.SetProjection(proj2)
ndvidrv.GetRasterBand(1).SetNoDataValue(-15000)
ndvidrv.GetRasterBand(1).WriteArray(arr1)
del ndvidrv

outQA = '/path/to/outputs/output2.tif'
qadrv = drvtif.Create(outQA, ncols, nrows, 1, 1) # 1= 8 bit
qadrv.SetGeoTransform(hc_gt)
qadrv.SetProjection(proj2)
## qadrv.GetRasterBand(1).SetNoDataValue(-15000)
qadrv.GetRasterBand(1).WriteArray(arr2)
del qadrv

outQA turns out as expected, with the correct range of values and correct projection, but outNDVI is just a raster of 0's (in the correct projection). The values should range from -10000 to 10000 (with NoData= -15000).

Any thoughts?

Deleting ndvidrv seemed to do the trick and I've updated above to reflect that.


texts: QGIS python script adds populates a column with the data, but the data doesn't appear when the shapefile is opened in ArcGIS

texts:
I know basic python, but am struggling to learn how to use it with QGIS. I have a polygon shapefile that maps the presence/absence of various species. Each species has a column in the attribute table a- if the species is present in the polygon then the species name is in the cell, otherwise "null".
I have a dictionary of fullname:abbreviation which I use to create a string for each row of abbreviations so that all the data for each row is represented in a single field.

The script seems to work, but when I tried to use the shapefile in ArcGIS the new column was empty. I saved the edits afterward. Can anyone tell me what I'm missing?

from PyQt4.QtCore import QVariant

subs = {'Argent & Sable':'AS',
'Barbastelle':'Bbl',
'Barred Tooth-striped':'BTs',
}
#etc etc

layer = iface.activeLayer()
layer.startEditing()
iter = layer.getFeatures()

for feature in iter:
    newStr = ''
    for field in feature.attributes():
        if field == 'NULL':
            pass
        else:
            if field in subs.keys():
                newStr += subs[field] + ', '

    feature['Species'] = newStr
    layer.updateFeature(feature)


texts: Is it possible to create a dynamic map in QGIS?

texts:
My table looks like this:

I’ve already done a Choropleth map with the first attribute of my table.

I would like to do the same Choropleth map with the rest of my attributes. I can’t replicate it manually because I have lots of attributes in my table.

Is there any possible way I could make a dynamic map where I select the category and the map gets coloured by that category? Qgis2web, Atlas?


texts: Performing Esri NAD27 to WGS84 datum transformations with ST_Transform in POSTGIS

texts:
I am trying to align with some data points that have been projected in Esri from NAD27 State Plane feet (SRID=32026) to Lon/Lat using the "NAD_1927_To_WGS_1984_6" datum transformation using ST_Transform in POSTGIS.

When using something like

ST_Transform( ST_GEOMFROMEWKT( 'SRID=32026;wkt_geom_string' ), 4326 ) )

I believe the default behavior is a NADCON grid-shift transformation, resulting in an approximate 5m shift to resulting features compared to their position from Esri resulting from above datum transformation being applied.

Therefore I thought I might need to use a custom string, and possibly even incorporate a pipeline, but so far, I have been unable to create a string that is syntactically correct.

The NAD_1927_To_WGS_1984_6 transformation is defined by Esri as Geocentric Translation with parameters dx=-8, dy=159, and dz=175.  It seemed that might equate to a Helmert 3 position transformation using Proj, but that presumption may be incorrect on my part.

Here's my benchmark point from Esri displaying it's resultant long/lat, using the mentioned datum transformation, that I wish to align with:

Source Esri PCS Point= 'SRID=32026; POINT(1443665.55622 684595.615785)' in units ft
Esri-projected GCS Point= 'SRID=4326; POINT(-122.672081324 45.523768294)' in dd using esri 'NAD_1927_To_WGS_1984_6' datum transformation

One approach I tried with POSTGIS was as follows:

WITH data AS ( SELECT '+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=-8,159,175 +nodefs'::text AS dtrans )
SELECT ST_AsEWKT( ST_Transform( ST_GEOMFROMEWKT( 'SRID=32026;POINT ( 1443665.55622 684595.615785 )' ), dtrans ) )
FROM data

Resulting in POINT(-122.670918714227 45.5218796715454).

Plugging the above two resultant long/lats into:

SELECT ST_Distance(gg1, gg2) As spheroid_dist
FROM (SELECT
    'SRID=4326;POINT(-122.672081324 45.523768294)'::geography as gg1,
    'SRID=4326;POINT(-122.670918714227 45.5218796715454)'::geography as gg2
    ) As foo;

Results in 228.71471422, which is way off.

I also tried command line cs2cs, attempting to build a pipeline that incorporated a helmert transformation, but was never able to get the syntax correct.

When I try:

echo 1443665.55622 684595.615785 | cs2cs -f "%.8f" +proj=pipeline +step +proj=lcc +lat_1=44.33333333333334 +lat_2=46 +lat_0=43.66666666666666 +lon_0=-120.5 +x_0=609601.2192024384 +y_0=0 +datum=NAD27 +units=us-ft +no_defs +step +proj=helmert +convention=coordinate_frame +x=0.67678 +y=0.65495 +z=-0.52827 +rx=-0.022742 +ry=0.012667 +rz=0.022704 +s=-0.01070 +to +proj=longlat +datum=WGS84 +no_defs

where

"+proj=lcc +lat_1=44.33333333333334 +lat_2=46 +lat_0=43.66666666666666 +lon_0=-120.5 +x_0=609601.2192024384 +y_0=0 +datum=NAD27 +units=us-ft"

is the proj4text for SRID=32026.

The resultant error is malformed pipeline, cannot instantiate the source coordinate system. I have tried several variations of this syntax with the same error. I have also scoured the Proj documentation and this board for similar issues with documented resolutions, but have come up empty.

Based on this background, my main questions would be:
Am I attempting to apply the correct transformation (helmert 3 pos) in Proj corresponding to a geocentric translation defined by Esri?
If so, does someone have an example of a Proj pipeline string so that I can see the correct syntax, or could help me correct the syntax in the examples I have provided?


texts: I am looking for some help drafting code to define the coordinate system for multiple las files. I have seen plenty on how to batch define projections with rasters and shapefiles, but nothing on las files. Is there a lis command for las data? Could that be used in a for loop?


texts: compute nearest distance point to point inside a polygon pyqgis

texts:
I need to compute the distance from an input point file to a series of target points in different shapefiles. Everything is in Africa. I want to do this in a pyQgis script. The script I have shown below does the job.

However, I would prefer if the aerial distances computed did not cross the sea. Is there a way to impose a restriction so that the distances are only computed to points inside a polygon (the African continent)?

import os
import glob
from qgis.core import * import qgis.utils

  processing.runalg("qgis:distancematrix","path-to-my-input-file","inputfileid",atts,"path-to-my-target-file",0,1,"path-to-output-file")


texts: Extract boundary from detailed layers in ArcGIS Desktop or QGIS?

texts:
Is it possible to return the boundary (outside) perimeter of a polygon layer?

Let's say you have the map of the USA, for simplicity, broken down to the county level, is it possible to return the boundary map of the whole country from that map?


texts: Running QGIS script in external Python code gives TypeError?

texts:
I have one piece of Python code that I wrote in QGIS 3.6 (I'll call this code1)
and I have a second piece of Python code that is a bigger project, which I wrote just in a normal text programme and which I run through the terminal of my xubuntu (I'll call this code2)

They both work seperately in QGIS / the terminal.
Now I want to combine them - code2 should at some point call code1 and run it as well. I don't want to open QGIS and run the code, everything should be done as one single process.

So I copied code1 into code2. Then I wrote this at the beginning of code2

from qgis.core import *
from qgis.gui import *
import qgis.utils
QgsApplication.setPrefixPath("/usr/share/qgis/",True)
qgs = QgsApplication([],False)
qgs.initQgis()

When the code reaches the first QGis function

params_proj = {
    "INPUT": layerinput_proj,
    "TARGET_CRS":"EPSG:4326",
    "OUTPUT":("/media/vali/project/geodaten/gem_proj.shp")
    }
feedback=QgsProcessingFeedback()
res_proj=QgsProcessingAlgRunnerTask(u"native:reprojectlayer", params_proj, feedback=feedback)
res_proj["OUTPUT"]

I get this error:

TypeError: QgsProcessingAlgRunnerTask(): argument 1 has unexpected type 'str'

Why does the code work in QGIS but not when I run it through the terminal?

I searched the Internet and could not find any solutions that work for me or that I can understand, as I am quite a beginner.

Since someone asked, here is a larger bit of my code:

#! /usr/bin/python3
import openpyxl, os, csv
from openpyxl.utils import get_column_letter, column_index_from_string
import xml.etree.ElementTree as ET, pandas as pd, geopandas as gpd, csv
from qgis.core import *# supply path to qgis install location
from qgis.gui import *
import qgis.utils

#------------------------------------------------------------------------
QgsApplication.setPrefixPath("/usr/share/qgis/",True)
# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([],False)
# load providers
qgs.initQgis()

#sys.path.append('/docs/dev/qgis/build/output/python/plugins')

path="/media/valentina/0bc36d5e-1c60-4266-a1ef-a284ac0b5315/project/geodaten/"

#-------------- change projection of gem shape -------------
layerinput_proj=QgsVectorLayer((path+"STATISTIK_AUSTRIA_GEM_20190101Polygon.shp"), "proj", "ogr")

params_proj = {
    "INPUT": layerinput_proj,
    "TARGET_CRS":"EPSG:4326",
    "OUTPUT":(path+"gem_proj.shp")
    }
feedback=QgsProcessingFeedback()
res_proj=QgsProcessingAlgRunnerTask(u"native:reprojectlayer", params_proj, feedback=feedback)
res_proj["OUTPUT"]

#-------------- calculate centroids -------------
layerinput=QgsVectorLayer((path+"gem_proj.shp"), "Gemshape", "ogr")

params = {
    "INPUT": layerinput,
    "ALL_PARTS":False,
    "OUTPUT":(path+"gem_centroid.shp")
    }

res=processing.run("native:centroids", params, feedback=feedback)
res["OUTPUT"]

#-------------- reproject centroids -------------
layerinput_proj=QgsVectorLayer((path+"gem_centroid.shp"), "proj", "ogr")
params_proj2 = {
    "INPUT": layerinput_proj,
    "TARGET_CRS":"EPSG:4326",
    "OUTPUT":(path+"gem_centroid_proj.shp")
    }

res_proj2=processing.run("native:reprojectlayer", params_proj2, feedback=feedback)
res_proj2["OUTPUT"]

#-------------- calculate geometry attributes -------------
layerinput_geom=QgsVectorLayer((path+"gem_centroid_proj.shp"), "gem_centroid", "ogr")

params_geom = {
    "INPUT": layerinput_geom,
    "CALC_METHOD":0,
    "OUTPUT":(path+"gem_centroid_geom.shp")
    }

res_geom=processing.run("qgis:exportaddgeometrycolumns", params_geom, feedback=feedback)
res_geom["OUTPUT"]

qgs.exitQgis()


texts: How to add legend in Qgis2Web plugin while exporting maps?

texts:
currently i am working with the publishing web maps with the Qgis2web plugin. but i didn't find any kind of legend window or event to add legend. Is there any other source to add Legend?below is snapshot of Qgis2Web parameters window. there i didn't find any event to add legend.

I need to display my all legend. e.g

to the final web map. suggestions are welcome


texts: Trajectory Similarity Measurement with Shapely or other Python libraries

texts:
Are any of the famous trajectory similarity measurement methods such as Dynamic Time Warping(DTW), Edit Distance based methods or Longest Common Subsequence based Measures implemented in shapely or any other python based libraries?

My end goal is to find similar sections between two trajectories. I do not care about the time stamps and when the trajectories are traversed. The lengths may be different. A short part of one trajectory may overlap with another trajectory.


texts: Nearest neighbor distance for Arctic area that spans 2-3 UTM zones

texts:
Using QGIS3.10, I have tried Distance Matrix and NNJoin (distance results 0 - 10.32 for both, likely in degrees). I have 93 points of 2 types in one shapefile. I have the same distance issue if the points are in 2 files. UTM projections are problematic in this region because the study area spans at least 2 and often 3 UTM zones (e.g., x: -63 to -78). The CRS is WGS84 (EPSG: 4326). The points (caribou locations) fall into 3 clusters 150-300 km apart, but one point is in a central position between the 3 clusters and it is about 154 km from its nearest neighbor (i.e., QGIS's Measure Line tool readily provides km for WGS84).
How could I resolve this to get nearest neighbor distances in km?


texts: Trying to define a "rotated regular grid" coordinate system in QGIS

texts:
I have some raster data whose metadata defines its coordinate system as:

rotated regular grid, virtual North Pole at 39.25 N, 162.00 W (rotated coordinates)

Its spatial resolution is 0.044 degree.

I am trying to get this data into regular longitude, latitude coordinates and to do so I have defined a custom projection with these parameters:

+proj=longlat +datum=WGS84 +lat_0=39.25 +lon_0=-162 +no_defs.

It doesn't really seem to help though, this is the raster layer and OSM basemap layer, both visualized in WGS84:

and this is the same layers, both visualized with this defined coordinate system:

I feel like I must be missing something obvious here.  Is this not how you would define a coordinate system for a rotated regular grid?

The source data, as well as the metadata PDF is here


texts:
I'm quite new to postgresql and postgis so this error will probably due to some rookie mistake.

I want to run the spatial query ST_within. When I describe the query like this

SELECT nw.entity_id AS NavElem_id, nw.featuretype, nw.processingstatus AS NavElem_ProcessingStatus, nw.frc, nw.geometry AS NavElem_geom, aaa.processingstatus AS AA_processingstatus, aaa.networktype, aaa.geometry AS AA_geom, aaa.entity_id AS AA_id
FROM sml.network nw, sml.aa_geo_att aaa
WHERE ST_Within (nw.geometry, aaa.geometry);

It works fine, however I do not want to provide the schema name everytime (in this case sml.), I want to make it generic so that I can use it in an sql file.

When I try to set the schema path in a first query and the ST_within in second query

SET search_path TO sml;
SELECT nw.entity_id AS NavElem_id, nw.featuretype, nw.processingstatus AS NavElem_ProcessingStatus, nw.frc, nw.geometry AS NavElem_geom, aaa.processingstatus AS AA_processingstatus, aaa.networktype, aaa.geometry AS AA_geom, aaa.entity_id AS AA_id
FROM network nw, aa_geo_att aaa
WHERE ST_Within (nw.geometry, aaa.geometry);

I get the following error:

ERROR:  function st_within(public.geometry, public.geometry) does not exist
LINE 5: WHERE ST_Within (nw.geometry, aaa.geometry);
          ^
HINT:  No function matches the given name and argument types. You might need to add   explicit type casts.
********** Error **********

ERROR: function st_within(public.geometry, public.geometry) does not exist
SQL state: 42883
Hint: No function matches the given name and argument types. You might need to add explicit type casts.
Character: 311

Even though I set the schema, ST within is looking in the public schema, which I don't want.


texts:
I have 1,000 point locations that represent hospitals, and I want them to be symbolized as different colored points. Since a map reader isn't going to tell the difference between 30 shades of red or blue, I can't have the same colors for neighboring points. Neighboring points, which represent different hospitals would imply the same organization or hospital institution (which isn't the case). What I've tried is using the Near tool in ArcGIS Pro to find the closest point for each point and what the distance is to it. Where I'm stalled is applying a staggered symbology, so that points within, say, 20-30 miles aren't given the same color.

I've also explored trying this with Python where I iterate through each feature point, use a buffer to define the area that I don't want to repeat the symbology, then apply a symbology, and move on to the next point. I'm not sure this is the correct path.

Does anyone have any other ideas for applying symbology based on proximity of features to each other?


texts: Add an event listener on a Marker in Leaflet

texts:
I am using Leaflet in order to render a map. I created a map with markers and I don't know how to implement the Event Listener 'onClick' on each Marker.

My code

var stops = JSON.parse(json);
var map = new L.Map('map', {
  zoom: 12,
  minZoom: 12,
  center: L.latLng(41.11714, 16.87187)
});
map.addLayer(L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token={accessToken}', {
  attribution: 'Map data',
  maxZoom: 18,
  id: 'mapbox.streets',
  accessToken: '-----'
}));
var markersLayer = new L.LayerGroup();
map.addLayer(markersLayer);

//populate map from stops
for (var i in stops) {
  L.marker(L.latLng(stops[i].Position.Lat, stops[i].Position.Lon), {
    title: stops[i].Description
  }).addTo(markersLayer).bindPopup("<b>" + stops[i].Description + "</b>").openPopup();
}

Example

map.on('click', function(e) {
    alert(e.latlng);
});

Leaflet deals with event listeners by reference, so if you want to add a listener and then remove it, define it as a function:

function onClick(e) { ... }

map.on('click', onClick);
map.off('click', onClick);


texts: How to load *.shp.xml file into shapefile for ArcGIS for Desktop?

texts:
I am using DRASTIC model and ArcGIS mapping for my thesis of groundwater assessment. I need to generate maps which includes Net Recharge map.

For that I need to export a landuse map and from the data I gathered, I cannot load its as seen below.

the only files in the folder are the following:


texts: Timeline and clockview model did not appear in cesiumjs

texts:
I have performed coding to the particular html file and it contain the code to access the wms layers from the Geoserver. However, when I run the html file, I am not able to see a clock view model and the timeline. Here is my code below.

<!DOCTYPE html>

  <meta charset="utf-8">

  <meta http-equiv="X-UA-Compatible" content="IE=edge">

  <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, minimum-scale=1, user-scalable=no">
  <title>Blue Marble</title>

Can you help me?


texts: Is it possible to enter start/stop points non-interactively in Road Graph plugin?

texts: I need to find the shortest paths in a network of roads for many pairs of start/stop points?  The Road Graph plugin in QGIS allows me to interactively select a starting and ending point on a map, but I can't do this interactively with hundreds of such points.  Is there a way to run this plugin in a non-interactive way?


texts: Changing Mouse Wheel zoom direction in QGIS?

texts:
I'm just starting to learn QGIS. I'm starting with the basic simple things. Coming from an ArcGIS background, I am used to the mouse wheel action zooming in and out the default ArcGIS way. QGIS default is the opposite. The only option I seem to find about changing mouse action options is to enable or disable the mouse wheel action, but nothing about the altering the action.

Can I do this?


texts:
I have a layer showing elevation curves and I'm trying to create an elevation profile using the Terrain Profile plugin. However, I only get a flat line. My guess is that I need to add Z values to the curves. The layer attribute table does have a column with height values(hoyde).

Can this be done?


texts: Field calculator fails when computing First X coordinate in ArcGIS10

texts:
I have this line in python to compute First X coordinates into the Field Calculator from each polyline:

!Shape!.firstPoint.X

Then, when I inspected the attributes I see this rare thing:

You may see that I have two fields (StartX and InicioX), both with the same Python code, one is a string field and the other a Float one. But the computed VALUE is wrong  :(

You can clearly see that "spatially" the features are situated in different positions, but on its computed X Value are the same.

Could someone explain me that?

Does anybody know any other way to compute this X,Y start and end coordinates from lines?

Thanks for your time.


texts: How to return all minimum paths for one to one with pgr_dijkstra?

texts:
Is it posssible?
I need to calculate all shortests paths for the same one to one, but in my topology exists severals solutions, not only one (a simply example is when all costs are 0).
How can I get all of them?
pgr_disjkstra always returns the same one solution and I need all of them.


texts: How to measure areas of features in OpenLayers plugin background maps correctly?

texts:
I can't figure out how to make areal calculations (eg. in km2). I have tried the 'field calculator', I have tried the 'measure area' function on the "attributes" toolbar, and i have tried fTools in the "Vector"-menu.

How I do:

Open Qgis
Add Google Hybrid Layer (OpenLayers Plugin)
Find tenniscourt (should be 200-260m2)
Measure

attribute toolbar->measure area - result: 65.060 m² hmm?
layer->new shapefile layer->polygon, save, "add feature" and drawing polygon. Then field calculator $area - result: 118.48828125
repeat with new shapefile - result: 4.05503328693158e-08
vector->geometric tools->export/add geometric columns = another random result

I can't figure out what I am doing wrong! For me it seems like there is some error in the program? I have tried saving the shapefile and opening it again in a "blank" project - same result - and new result every time?


texts: Dynamically adding in-memory QgsVectorLayer triggers a popup in the UI - how to eliminate

texts:
I use this function to create layers that display results to my queries. However a dialog always pops up asking me to confirm the crs. Is it possible to suppress this ?

    def addResultLayer(name, feature, symbol, groupname):
        res_abs_path = self.getStubVectorFile()
        result_layer = QgsVectorLayer("Polygon", name, "memory")
        QgsMapLayerRegistry.instance().addMapLayer(result_layer, False)
        result_layer.dataProvider().addFeatures([feature])
        result_layer.rendererV2().setSymbol(symbol)
        result_node = QgsLayerTreeLayer(result_layer)
        self.findGroup(groupname).insertChildNode(0, result_node)


texts: How would one format an Excel sheet so numeric fields are set to numeric field types when importing to ArcGIS 10.2?

texts:
I have a large Excel spreadsheet that contains lots of numeric data that I would like to import to ArcGIS 10.2 as a table using the Excel To Table (Conversion) tool.

The frustrating thing is that the numeric fields in Excel are repeatedly not recognised and are assumed as being text. I have set the cell format to numeric, cleared any possible blank cells, tried importing as CSV, removed invalid characters in the field headers. yet any possible iota of success evades me.

Does anyone know why this may be happening and/or how it can be resolved?


texts: Non-Square Pixel on georeferencing of images using GDAL

texts:
I have a image in .jpg format. I have some GCP and also know the pixels values of these GCp in that image. Now i trying to convert this image into geotiff. for that I have coded the following lines. But the output seems to be weird. Although the geo information are tagged in the images but the pixels in the output seems to have non-square pixel. Anybody help please???

My code :

     GDALDatasetH hSrcDSarray;
     GDALAllRegister();
     hSrcDSarray = GDALOpen("test.jpg", GA_ReadOnly );
     hDriverarr[i] =  GDALGetDriverByName( "GTiff" );
     GDALDatasetH hDstDSarray;
     char **papszOptions = NULL;
     GDALProgressFunc    pfnProgress = GDALTermProgress;

    hDstDSarray= GDALCreateCopy( hDriverarr, "geo.tif", hSrcDSarray, FALSE, papszOptions, pfnProgress, NULL );

    double w=GDALGetRasterXSize( hSrcDSarray ); //get width
    double h=GDALGetRasterYSize( hSrcDSarray ); //get height

    // The GCP and pixels values come from external file

    GDAL_GCP *gcps = NULL;
    gcps = (GDAL_GCP *) malloc(((10000-2)/4)*sizeof(GDAL_GCP) );

    for(int i=0;i<image_output.size();i++){
        gcps[i].dfGCPPixel = image_output[i].x;
        gcps[i].dfGCPLine  = image_output[i].y;

        gcps[i].dfGCPX = GCP_points[i].x;
        gcps[i].dfGCPY = GCP_points[i].y;
        gcps[i].pszId = (char *) malloc( 1000*sizeof(char) );

    }

    double adfGeoTransform[6];
    GDALGCPsToGeoTransform(image_output.size(),gcps,adfGeoTransform,0);

    free(gcps);

    OGRSpatialReferenceH hSRS;
    char *pszSRS_WKT = NULL;

    GDALSetGeoTransform( hDstDSarray, adfGeoTransform );

    hSRS = OSRNewSpatialReference( NULL );
    OSRSetUTM( hSRS, UTM_zone, TRUE );
    OSRSetWellKnownGeogCS( hSRS, CS );
    OSRExportToWkt( hSRS, &pszSRS_WKT );
    OSRDestroySpatialReference( hSRS );

    GDALSetProjection( hDstDSarray, pszSRS_WKT );
    CPLFree( pszSRS_WKT );

    GDALClose( (GDALDatasetH) hDstDSarray );


texts:
QGIS crashes without error or log massage as soon I try to delete a column in a GeoPackage using the attribute table:

How could I find out, what's going on in background of QGIS?

Installed QGIS3 on Mint 18.3 (Xenial) with ubuntu-gis-unstable repos


texts: How to get exact featureinfo in wms layer Geoserver

texts:
I got neighbor polygon feature info instated of clicked polygon feature info.
In My Code:

getfeaturenfo :function(e) {

            coordinates=e;
            var params = {
                        REQUEST: "GetFeatureInfo",
                        projection: "EPSG:4326",
                        EXCEPTIONS: "application/vnd.ogc.se_xml",
                        BBOX: map.getExtent().toBBOX(),
                        SERVICE: "WMS",
                        INFO_FORMAT: 'text/html',
                        QUERY_LAYERS: selectlayer(map.zoom),
                        FEATURE_COUNT: 50,
                        Layers: selectlayer(map.zoom),
                        WIDTH: map.size.w,
                        HEIGHT: map.size.h,
                        X:parseInt(e.xy.x),
                        Y:parseInt(e.xy.y),
                        CQL_FILTER:selectfilter(),
                        srs: map.layers[0].params.SRS};

                    // handle the wms 1.3 vs wms 1.1 madness
                    if(map.layers[0].params.VERSION == "1.3.0") {
                        params.version = "1.3.0";
                        params.i = e.xy.x;
                        params.j = e.xy.y;
                    } else {
                        params.version = "1.1.1";
                       params.y = parseInt(e.xy.y);
                        params.x = parseInt(e.xy.x);
                    }
                    OpenLayers.loadURL("http://myserver/tomcat/geoserver/india/wms", params, this, setHTML, setHTML);

In Layer Preview I got same issue,So is that problem with layer or with code?

N.B: If Polygon is big then it shows perfect featureinfo but when it's small then this issue occurs.
How Can I solve this issue?


texts:
I am running QGIS 3.4.14 on Ubuntu Linux 18.02. I cannot see the SAGA tool in the QGIS Processing window.  I can confirm that I have SAGA 7.3 installed and working as a stand-alone GUI application.  I have gone into the processing options in QGIS and enabled SAGA as a provider.
This QGIS website states  if you are running Linux, SAGA binaries are not included with Processing, so you have to download and install the software yourself. Please check the SAGA website for more information.
In this case, there is no need to configure the path to the SAGA executable, and you will not see those folder entries. Instead, you must make sure that SAGA is properly installed and its folder is added to the PATH environment variable. Just open a console and type saga_cmd to check that the system can find where the SAGA binaries are located.
I can confirm saga_cmd will open the command line version of SAGA.  I have added PATH="$HOME/bin:$PATH:/opt/saga-7.4.0" to ~/.profile.
Still, I do not see SAGA in the processing window.  How do I get the SAGA tool available in QGIS 3.4?


texts: Extract geometry from feature not working?

texts:
This works:
(with geometry being a Polygon and sent3 an ImageCollection):
var image = sent3.filterBounds(roi).first();
var mask = image.select("Oa07_radiance").mask().rename('mask');
var area = ee.Image.pixelArea().multiply(mask).rename('area');

var sumDictionarypolygon1 = mask.addBands(area).reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 300,
  maxPixels: 1e9
});
print('sum for pol 1 mask, scale=300', sumDictionarypolygon1);

While this does not (roi being a Feature containing a geometry)
var rois = ee.FeatureCollection("users/guillaumh/embalsescba"),
    sent3 = ee.ImageCollection("COPERNICUS/S3/OLCI");
var roi = rois.toList(99).get(26);
print(roi);

var image = sent3.filterBounds(roi).first();
var mask = image.select("Oa07_radiance").mask().rename('mask');
var area = ee.Image.pixelArea().multiply(mask).rename('area');

var sumDictionarypolygon1 = mask.addBands(area).reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: ee.Geometry(roi),
  scale: 300,
  maxPixels: 1e9
});
print('sum for pol 1 mask, scale=300', sumDictionarypolygon1);

The error is:
Dictionary (Error)
Feature, argument 'geometry': Invalid type. Expected type: Geometry. Actual type: Feature.

Why cannot I extract the geometry from the feature? It seems ee.Geometry(  ) does not do anything.