texts: How read a raster [Tiff] file to get the correct values for every band?

texts:
I am trying to read a raster tiff [RGBA] file with GDAL in [C++], I can load the file, get  the number of channels, the size of the bands and the metadata, but when I try to extract the Band [R, G or B] I get the images in gray scale but all of them has the same values; so when I try to merge them with opencv I get a Gray image and not the RGB, The code that I use to do this process is below... Could you give me an idea of what I am doing wrong?

I read the following answers and the GDAL doc but dont get a good result...
Link 1: Correct way to extract pixel data using GDAL. Current method gives repeated but changing value.

Link 2: https://stackoverflow.com/questions/33807705/create-an-rbg-color-qimage-from-a-gdaldataset?fbclid=IwAR1wIsteSrEyfZGGQW9s9gK4GyIYi6AIvFR_47FsDKbgfSEv5hrhsznT1ss  [I use this code and also try the GDALRASTERIO method]

Code:

    // Register all Drivers
    GDALAllRegister();
    // Load the raster
    GDALDataset *dataset = (GDALDataset *) GDALOpen(path, GA_ReadOnly);
    // Get some features of the raster
    int raster_cols = dataset->GetRasterXSize(),
            raster_rows = dataset->GetRasterYSize(),
            raster_channels = dataset->GetRasterCount();
    // keep data - Channels and Type of the image
   std::vector<std::string> raster_bandName(raster_channels),
                            raster_BandType(raster_channels);
   std::vector<cv::Mat> Mybands(raster_channels, cv::Mat( raster_rows, raster_cols, CV_8UC1 )); // Byte[Raster] -- 8U [OpenCV]
   cv::Mat img;
    // Walk over the bands ---
    for(int i = 1; i < raster_channels; i++)
    {
        // Get the data of the available bands
        GDALRasterBand *poBand = dataset->GetRasterBand(i);
        raster_BandType[i-1] = GDALGetDataTypeName(poBand->GetRasterDataType());
        // Get some descriptor of the band that was loaded
        raster_bandName[i-1] = GDALGetColorInterpretationName(poBand->GetColorInterpretation());

        std::string ban = raster_bandName[i-1];

        poBand->RasterIO(   GF_Read, 0, 0, raster_cols, raster_rows,
                            Mybands[i-1].data, raster_cols, raster_rows,
                            GDT_Byte, 0, 0 ); // Write the Image data in the OpenCV var

        std::cout<<raster_BandType[i-1]<<" "<<raster_bandName[i-1]<<std::endl;

    }

    cv::merge(Mybands,  img );

Update: The following code works for me, I don't see a real difference but It returns the desired image.....

GDALAllRegister();
GDALDataset *dataset = (GDALDataset *) GDALOpen(path, GA_ReadOnly);
int raster_cols = dataset->GetRasterXSize(),
                  raster_rows = dataset->GetRasterYSize(),
                  raster_channels = dataset->GetRasterCount();

    for(int i = 1; i < raster_channels; i++)
    {
        // Save the channel in var not in the vector of Mat
        cv::Mat bckVar(raster_rows, raster_cols, CV_8UC1  );

        GDALRasterBand *poBand = dataset->GetRasterBand(i);

        poBand->RasterIO(   GF_Read, 0, 0, raster_cols, raster_rows,
                            bckBar.data, raster_cols, raster_rows,
                            poBand->GetRasterDataType(), 0, 0 );
        // Image is a std::vector<cv::Mat> image;
        image.push_back( bckBar.clone() );
        // Clean mm
        bckBar.release();
    }

    cv::merge( image, usrImage2w );// RGB ok


texts: How can I include HTML in a Folium Marker popup?

texts:
I'd like to include an anchor tag in a Folium Marker popup so that I can click through to another page.
Here's my code:
    def poptext(row):
        return "" + row['Name'] + ""

    map = folium.Map(location=[52.620,-0.264], zoom_start=14)

    df.apply(lambda row: folium.Marker([row['Lat'], row['Lon']],popup=poptext(row)).add_to(map), axis =1)

    map

Unfortunately, the HTML doesn't seem to be evaluated. Any solution?


texts: Sorting polygons into a n x n spatial array

texts:
I have a square group of contiguous polygons:

By using the centroid coordinates of each polygon, I want to assign each polygon's spatial sequence in a row by row fashion:

I have tried this PostgreSQL statement:

select id,st_centroid(geom) as center
from work.qtrs_in_section
order by st_y(st_centroid(geom)),st_x(st_centroid(geom));

It does give me groups of sorted y-values, but the x-values are not sorted in each row group:

Ideas?


texts: Find polygon in grid line network PostGIS

texts:
I have a grid line table to use as the blade for splitting my polygon table in PostGIS with a ST_Split method. It works fine but now I want to find which parts of polygons are located within my grid.I thought about using ST_Centroid to find the centroid of split parts and then applying ST_Within. My question is: can I apply these two steps in one SQL query or should I compute centroid first and then use ST_whitin? Or is there anything else to do this?
Maybe something like this:

SELECT
    gridpoly.*
FROM
    (SELECT
        st_centroid(gridpoly.geom)
    FROM gridpoly) as center
JOIN
    ( SELECT geom FROM gridline As cells WHERE col=1 AND row=2 ) as line
ON ST_Within(center.geom,line.geom)

The gridpoly table is the result of the split function which contains polygons. But in this code I don't have any column as 'center'.


texts: Resampling hyperspectral data to Landsat 5 TM bands

texts: The convolution of hyperspectral data to Landsat 5 TM bands is done by using the spectral response function (SRF, or Relative Spectral Response [RSR]) of Landsat 5 TM sensor. The ranges of each valid band are determined by the full width at half maximum bandwidth (FWHM) of each SRF. The central wavelength is also determined from this range. However, when doing the spectral convolution, should I consider just the SRF whithin the valid band range (red line) of Landsat 5 TM sensor, or the whole SRF which also lies out of the valid band range (the response of red and black lines)? I have convolved using both alternatives, but I'm not sure which one is the right.


texts: Creating Legends for QGIS 2.0 Feature Blend Mode Outputs

texts: I am working on visualizing ~4100 buffered point features. The shading of the buffers is keyed to a field in the data. Because the buffers overlap, I turned on multiply in the new feature blend mode options. I like the visual result, but what seems to be missing is an ability to generate a meaningful legend for the resulting image. Are there any workarounds for this? Is the only alternative to venture down the path of polygon intersects and symmetrical difference operations?


texts: In QgsVectorDataProvider.deleteAttributes() set() problem

texts:
When I script the code:

def removeFields(provider, keepFieldNames, fileName, isBatch):
    """Remove fields other than keepFieldNames from shapefile fileName with provider."""
    toDelete = set()
    fields = provider.fields()
    for idx in range(fields.count()):
        name = fields.field(idx).name()
        if not name in keepFieldNames:
            toDelete.add(idx)
    if len(toDelete) > 0:
        OK = provider.deleteAttributes(fields)

I get this error:

File "C:/.../Topology.py", line 1028, in removeFields
    OK = provider.deleteAttributes(fields)
TypeError: QgsVectorDataProvider.deleteAttributes(): argument 1 has unexpected type 'set'


texts: How to create a single feature/segment from connected line segments

texts:
I want these 3 connected line segments to be a single feature/segment.

I used dissolve and unsplit tools but I am not getting desired output.

I dont want it done manually with merge tool, because I have lot of features like these.

I am using ArcGIS Desktop 10.1 (ArcInfo license).


texts:
I am using geoserver 16.0 stored in Tomcat and arcgis Runtime SDK 100.6. I am using a WMTS protocol on the communication and I am connection directly with the gwc:
"http://localhost:8070/geoserver/World/gwc/service/wmts?REQUEST=getcapabilities"
My problem is that when trying to get tiles from the gwc, it's always returning a
geowebcache-cache-result = MISS, even though the folder of gwc is growing up with the cached tiles, but it does not give me any further info on what's the problem.
When using a 12.something version of geoserver directly with the windows installer, not storing it in Tomcat, it is getting a HIT on the gwc result, but I need to use last version because of an improvement it has.


texts: Show all data of the WMS GeoServer layer on single click at web page using OpenLayers 6 WMS GetFeatureInfo (Tile Layer)

texts:
I have a vector layer in my GeoServer, which I want to publish on the web page.

When I open this layer manually from GeoServer, it shows the attribute table of the place/polygon, while clicking on the vector layer,
But I want to show it on web page, for doing so, I am using OpenLayers 6 (v6.4.3)
I have followed the instructions from following link
Example - WMS GetFeatureInfo (Tile Layer)
but the developer console on Chrome is showing just 1st row of the table

Where I am doing it wrong?
This is my HTML script for the web page,

    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <meta http-equiv="X-UA-Compatible" content="ie=edge">
    <title> Test AGRO DSS page - WMS GetFeatureInfo (Tile Layer) </title>

    <link rel="stylesheet" href="http://localhost:8080/agrodss_test/styles.css">
    <link rel="stylesheet" href="http://localhost:8080/agrodss_test/libs/v6.4.3-dist/ol.css">

 GFS District Forecast

&nbsp;

This is my JavaScript for OpenLayers,
    window.onload = init;

function init(){

    //Geoserver layer
    // wmsSource = indialayersource

    var indialayersource = new ol.source.TileWMS({
        url: "http://localhost:8080/geoserver/agrodss/wms",
        params: {"LAYERS":"agrodss:GFS_Forecast_Agromet_DSS_2020-12-02", "tiled": true},
        serverType: "geoserver"
    })

    // wmsLayer = indialayer

    var indialayer = new ol.layer.Tile({
        source: indialayersource,
        visible:true,
        title: "ForAgrodss"
    })

    // view = myview

    var myview = new ol.View({
        center: ol.proj.fromLonLat([80, 22]),
            zoom: 4.5,
            maxZoom: 8,
            minZoom: 4,
      });

      var map = new ol.Map({
        layers: [indialayer],
        target: 'js-map',
        view: myview,
      });

      map.on('singleclick', function (evt) {
        document.getElementById('info').innerHTML = '';
        var myviewResolution = /** @type {number} */ (myview.getResolution());
        var url = indialayersource.getFeatureInfoUrl(
          evt.coordinate,
          myviewResolution,
          'EPSG:3857',
          {'INFO_FORMAT': 'text/html'}
        );
        if (url) {
          fetch(url)
            .then(function (response) { return response.text(); })
            .then(function (html) {
              document.getElementById('info').innerHTML = html;
            });
        }
      });

      map.on('pointermove', function (evt) {
        if (evt.dragging) {
          return;
        }
        var pixel = map.getEventPixel(evt.originalEvent);
        var hit = map.forEachLayerAtPixel(pixel, function () {
          return true;
        });
        map.getTargetElement().style.cursor = hit ? 'pointer' : '';
      });

    // layer group

     baseLayerGroup = new ol.layer.Group({
        layers: [
            indialayer
        ]
    })

    map.addLayer(baseLayerGroup);

    map.on('click', function(e){
        console.log(e.coordinate);

    })
}


texts: GeoServer LDAP Authentication 'Test Connection' succeeds, but login fails

texts:
I am attempting to use Active Directory for GeoServer authentication.  When I click the 'Test Connection' button on the Admin > Authentication > Providers > LDAP settings page, it says "Connection Successful."  However, when I log out of Geoserver, and try to log in with the same credentials, I get an "invalid username/password combination" error.

My LDAP settings are as follows, where ADserver=my Active Directory server, and domain=my company domain:

Server URL: ldap://ADserver/dc=domain,dc=com
TLS: [Not checked]
User lookup pattern: [empty]
Filter used to lookup user: (userPrincipalName={0})
Format used for user login name: {0}@domain.com
Use LDAP groups for authorization: [Not checked]
User Group Service: default

And I have selected this Authentication Provider in my Provider Chain.

Why does the Test Connection succeed, but actually logging in fails?


texts: How to join different line segments between points into one entity in QGIS

texts: I have a shapefile of pipe lines. I want to join all the pipes between two valves to be considered as one segment in QGIS.


texts: Performance degradation in arcpy.mapping when setting text elements

texts:
I'm generating a map book of about 300 pages using arcpy, instead of using the toolbox, because I need to do a bunch of scripted actions - on each page I'm setting the extent and rotation of 3 DataFrames, setting the value of 8 TextElements, and then saving the layout using arcpy.mapping.ExportToPDF.

What I'm seeing is that the loop gets exponentially slower the longer I leave it running, and I've narrowed it down to setting the value of TextElements. I've got a simplified example where all I do is

Open the MXD
Grab a reference to the TextElement
Within an endless loop:

Set the value of the TextElement
Print out the time

Processing the times results in a graph like this:

...and it just keeps getting slower.

Has anyone seen this, and have a solution? Or do I just break up generating my map book into batches of 50-100 pages and restart the script? (which is my current workaround)

Simple example script:

import arcpy
import datetime
import os

MXD_MAPBOOK = r'E:\workspace\data\N&I\N&IMapBook_landscape2.mxd'

def findTextElement(doc, elementName):
    lstElements = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT")
    for elm in lstElements:
        if elm.name == elementName:
            return elm

mxd = arcpy.mapping.MapDocument(MXD_MAPBOOK)
dTime = datetime.datetime.now()
elem = findTextElement(mxd, 'baseRailLine')

while True:
    elem.text = 'Test number'
    print(str(datetime.datetime.now()))


texts: Lambert conformal conic projection - Europe

texts:
I am trying to calculate map of Europe in Lambert projection. I am using formula from Wikipedia (Lambert wiki), but calculation is very slow in comparison with Mercator. I know, there are sin, cos, cotg and pow to calculate (Mercator has only one tan and ln), but the difference is almost 30 minutes (Mercator is done in about 3 minutes).

Is there some kind of speed up, or similar projection, that can be used ? I also need to put result into image with width and hight. For Mercator, I am using this: https://stackoverflow.com/questions/5983099/converting-longitude-latitude-to-x-y-coordinate. Can it be altered to Lambert (or Lambert similar) projection ? Because if I use same approach only witch changed formulas for projection, Lambert is moved to upper part of image and I can see only part of Europe)

EDIT: Added code

double longitude = lon * (MPI / 180.0);
double latitude = lat * (MPI / 180.0);

double lambertPhi = this->lambertF * pow(cot(0.25 * MPI + 0.5 * latitude), this->lambertN);

lonOut = lambertPhi * sinf(this->lambertN * longitude);
latOut = this->lambertPhi0 - lambertPhi * cosf(this->lambertN * longitude);

All this-> variables are precomputed only once


texts: Handling Popups in Leaflet with Overlapping Features?

texts:
I am new to Leaflet as well as Javascript.

Currently, I am trying to create a spatial index of library holdings for old topographic maps that people can access and download online. The source of the features is a GeoJSON file.

Sample Map

The problem is that the school library has multiple types of the same map over the years, but when I click on a map, only one popup shows up. I want the option to cycle through multiple maps, but I am unsure on how to approach the problem.

Is there a special term in the local jargon for cycling through overlapping polygons, or is there a stronger approach to this problem?


texts: Software to easily go CSV - map - CSV for editing

texts:
I am looking for a piece of software that will easily let me display lat long points from a CSV file on a map, then be able to select a group of points (via a box or drawing), delete the points and resave as a CSV.

ArcGIS requires way too many steps to do this over and over again
Power map in Excel shows the points but can't be edited
Tableaux shows the points on a map, allows me to filter the bad points out but have not found a way to then resave the filtered data to a new CSV.

This should be simple.


texts: Can a Bluetooth barcode reader update a Geodatabase via ArcGIS Server on iPad?

texts:
I need the ability for a field operator (using an iPad) to select a feature in a geodatabase served via ArcGIS Server, and scan its barcode, which should be entered into the geodatabase as an attribute.

1) Will Esri's ArcGIS app handle this capability - can it talk to an external barcode reader, and use the output of the barcode reader to edit an attribute?

2) Assuming that it can't and we need to write an our own app, are there any mechanisms within Esri's iOS API which allow talking to hardware devices?

The ideal scenario would be a wizard-based approach which prompted the user to scan the barcode when a suitable feature was clicked, then entered the resulting barcode text as an attribute. A more manual, step-by-step workaround would be acceptable if necessary.

I've assumed that an external barcode reader, connected via Bluetooth, will be necessary as it's impractical to use the iPad's own camera due to accessibility issues (eg pits).

Thanks


texts: Python conditional on current and previous row?

texts:
I need to calculate a new field with a python conditional where some values only occur if a sequence in a separate attribute field (X) occurs.
For instance my data set looks like this. Y is the attribute I want to fill out:

ID   X  Y
1    1  -
2    1  -
3    0  -
4    0  -
5    0  -
6    1  -
7    1  -
8    1  -
9    0  -
10   0  -
11   1  -

Values for Y:
If X = 1, Y = 1
If X = 0 and previous X = 1, Y = 2
If X = 0 and previous x = 0, Y = 0

So my data would look like:

ID   X  Y
1    1  1
2    1  1
3    0  2
4    0  0
5    0  0
6    1  1
7    1  1
8    1  1
9    0  2
10   0  0
11   1  1

I have seen previous lines stored using:

prev_X = ""

But don't know how to incorporate it into my script which looks like:

fc = r'.....'

rows = arcpy.UpdateCursor(fc)

Y = 1
prev_X = ""  #<=========== DOESN'T WORK
for row in rows:
    if row.X == 1:
        Y = 1
        row.Y = Y
        rows.updateRow(row)

    while row.X == 0 and prev_X == 1:
        Y = 2
        row.Y = Y
        rows.updateRow(row)

    while row.X == 0 ande prev_X == 0:
        Y = 0
        row.Y = Y
        rows.updateRow(row)
        rows.next()

    row.X = X
    rows.updateRow(row)

del row, rows


texts: Is it possible to use the L8 OLI/TIRS COLLECTIONS LAND SURFACE REFLECTANCE Landsat8 product for calculation of the surface temperature?

texts:
Is it possible to use the L8 OLI/TIRS COLLECTIONS LAND SURFACE REFLECTANCE Landsat8 product for calculation of the temperature of the Earth surface?

I examined the following sources about this product:

Landsat 8 OLI/TIRS Level-2 Data Products - Surface Reflectance
https://lta.cr.usgs.gov/L8Level2SR

Landsat Surface Reflectance Level-2 Science Products
https://landsat.usgs.gov/landsat-surface-reflectance-data-products

PRODUCT GUIDE LANDSAT SURFACE REFLECTANCE-DERIVED SPECTRAL INDICES
https://landsat.usgs.gov/sites/default/files/documents/si_product_guide.pdf

But I didn’t find information about the Surface Temperature Data.
Also I download an example of this product and can’t find the surface temperature there.

However, one of the main purposes of TIRS bands of Landsat-8 – is getting the surface temperature. Also, it’s rather difficult to calculate the surface temperature on its own, since many procedures have to be done by users (corrections, calculations).

Thus, I suppose that L8 OLI/TIRS COLLECTIONS LAND SURFACE REFLECTANCE product should have the surface temperature and I just can’t find it.

May somebody help me to understand this subject.


texts: Spatial join points and polygons in BigQuery?

texts:
I have two tables; one has points (e.g., specific addresses with one geography POINT per row) another has areas (e.g., neighbourhoods with one geography POLYGON per row).
I'd like to join them, to have the neighbourhood associated with the address. Following this question I've tried the following:
SELECT
  pp.full_address,
  gn.name,
  ST_CONTAINS(gn.geography, pp.geopoint) as is_contained,
  pp.geopoint,
  gn.geography
FROM
  addresses pp,
  neighbourhoods gn
WHERE
  ST_CONTAINS(gn.geography, pp.geopoint)

But it doesn't seem to work. I get points that clearly aren't contained in the polygon. Note the is_contained column in the result is always True, but if I copy & paste a polygon/point pair into a standalone query like SELECT ST_CONTAINS(polygon, point) it returns False.
What am I missing?


texts:
Due to problems with creating AGS Services I have been told I should do a "repair" of ArcGIS 10.3.1 for Server.

Does anyone know if there is a utility or tool to do this?


texts: Ok to measure and union a projected-on-the-fly shapefile in ArcMap?

texts:
I have a dataframe in arcMap 10.0.4 that is of Africa Albers Equal Area Conic PCS. I bring in a census shapefile of WGS 1984 CS and it is projected on the fly.
I do 2 things:

I use the calculate geometry option to overwrite the "Shape_Area" field in the census shapefile so that I may obtain the planar measurements of the shape area (i.e. square Km as opposed to Decimal Degrees).
I 'union' an existing polygon grid projected as Africa Albers Equal Area Conic with the WGS'84 Census shapefile. The result of this union is a shapefile with the Africa Albers Equal Area Conic PCS.

I ask both these questions because they seem to me to stem from the same issue: error propagation. Am I creating measurement and geometric errors because I should be instead projecting the census shapefile to Africa Albers Equal Area Conic PCS? What is the capabilites of ArcMap to handle this issue?


texts: Edit multiple tables with the same structure at the same time?

texts:
does anyone know how to edit multiple tables at the same time? I have multiple .tab files and I want to change the order of the fields. Is there any way to do this for all of them at once. They all have identical table structures. I have both MapInfo and QGIS available to use.

Thanks, Alan


texts: Leaflet Transform draw Marker from data Ionic 2

texts:
I want to draw markers from data on Leaflet map

I'm currently doing with that way

                 L.marker(
            [48.840346,2.319415]).addTo(map).bindPopup("X"
 ",{minWitdh:256} ");

But i have a Data.json

  {
      "title": "A",
      "lat": 48.841491,
      "long": 2.355989,

   },
    {
      "title": "B",
      "lat": 48.878849,
      "long": 2.312855,
      },

i want to do something like this

 L.marker(
        [lat,long]).addTo(map).bindPopup("title"
 ",{minWitdh:256} ");

To get all The Data and render it on my Map ?


texts: Creating point layer from latitude/longitude table in ArcMap?

texts:
I want to create points from lat/long data. I used "File" -> "Add Data" -> "Add XY Data..." to create the points. I am 100 % positive that a I used the right columns for X and Y and I selected the correct CRS. I checked multiple times.
However, when creating the points they are just centimeters apart from each other, using the Measure tool to measure distance. That cannot be correct.
I also used "Add Layer" -> "Add Delimited Text Layer" in QGIS and got the correct locations of the points. It would be better for me to do the action in ArcMap though, since I need to load over 300 of the files and want to build a model later on to do that.

The coordinates are:

The CRS is ETRS_1989_UTM_ZONE_33N
I use ArcGIS 10.3 on Windows 7

What am I doing wrong/what could be the reason for the misplacement of the points in ArcMap?

I know that the points should be around 570 m apart from each other.


texts: Python and Overpass: Cannot get a complete address list even though I can see it on the map

texts:
I am new to OSM and running a Python script like this https://pypi.org/project/osmxtract/ to get building addresses and polygons of my street via Overpass.

I get a list. The list changes depending on my buffer but in weird ways: a small buffer yields no addresses, a large buffer generates a long list but does not include my address or most other local streets.

Is the buffer also changing the zoom so the server has to dump less detail for a bigger area?
Should I be specifying values to keep the building list short and limit my features?

lat, lon = location.geocode('Stretford, Manchester')
bounds = location.from_buffer(lat, lon, buffer_size=1000)

query = overpass.ql_query(bounds, tag='building')
response = overpass.request(query)

feature_collection = overpass.as_geojson(response, 'polygon')
with open('building.geojson', 'w') as f:
    json.dump(feature_collection, f)

geodataframe = gpd.GeoDataFrame.from_features(feature_collection)


texts: Strange behaviour of piechart function of QGIS

texts:
I am using QGIS to visual health status of elderly in an elderly home, whether they have medical devices or not is one attribute that I am trying to visualize.

I have create the map using mapwindows, and use R to create a tab-delimited txt file to carry the medical device information (7 columns, 6 columns of 0 or 1, and one uid that will be used to join the bed layer in QGIS.

I used the diagram tab in the properties menu of the bed layer, and pie-chart is being used for visualization: if the elderly have 3 devices, the pie will be presented with 3 parts with 3 different colours. And strange things are visualized.

A majority of elderly do not have any devices (i.e. that should be straight zeros for the 6 column in the txt file). I support no pie chart should be visualized, but QGIS still show a solid pie using the last attribute I have inserted in the diagram tab of the properties menu.

I have tried to insert a csvt file to explicitly state that those medical device columns are integer, real, but it still fails. QGIS visual an attribute with value 0 as 1, if all of the attributes used for plotting are straight zeros.


texts: Symbolizing multiple polygons as one by dropping internal boundaries?

texts:
Using ArcMap (ArcGIS 10.2 for Desktop), is there a way to symbolize two features within the same feature class as one feature without editing the feature and merging them?

For example, if I have a map of the United States (feature class) and all the states are individual features within this feature class and I want to symbolize the United States as a whole without the state divisions, could I do this?

I realize I could remove the outline, but that would also remove the outline from the United States as a whole, and I don't want to do that.

I feel like there is a simple answer for this, that I just have not discovered.


texts: How to obtain nugget parameter from fit.variogram in gstat?

texts:
I have the following code where I fit a variogram using gstat.

vA = variogram(dep_var~elevation=x+y, data, cressie=TRUE)
vmodel = vgm("Sph")
vB = fit.variogram(vA,vmodel)

How can I get the nugget parameter from the 'vB' variable?

I am using vB$psill[1] as the nugget but it is large.


texts: How do I change extent of shapefile into lat lon of NetCDF file?

texts:
I have loaded a shapefile and NetCDF into R, and I want to extract data based on the extent of the shapefile. I have done the following, but when I want to extract it I get NA's as the extent differs between the file:
library(rgdal)
library(sf)
library(raster)
library(ncdf4)

shp = readOGR("shapefile.shp")
net = brick("netCDF.nc")
extent(net)
extent(shp)

crs(shp)
crs(net)

net.mask = mask(net.brick, shp)
net.df = as.data.frame(net.mask[[1]], xy = TRUE)
head(temp.df)

> shp
class       : SpatialPolygonsDataFrame
features    : 1
extent      : -2939454, -1610437, 3261153, 4333895  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs
variables   : 8
names       : NA_L2CODE,              NA_L2NAME, NA_L1CODE, NA_L1NAME,                    NA_L2KEY, NA_L1KEY,    Shape_Leng,   Shape_Area
value       :       3.1, ALASKA BOREAL INTERIOR,         3,     TAIGA, 3.1  ALASKA BOREAL INTERIOR, 3  TAIGA, 11989512.2816, 423471100811

> crs(shp)
CRS arguments: +proj=longlat +datum=WGS84 +no_defs
> crs(net)
CRS arguments: +proj=longlat +datum=WGS84 +no_defs

> extent(net)
class      : Extent
xmin       : -180.125
xmax       : -51.875
ymin       : 40.875
ymax       : 90.125

> extent(shp)
class      : Extent
xmin       : -2939454
xmax       : -1610437
ymin       : 3261153
ymax       : 4333895

How can I change the extent of the shapefile so it gets latitude and longitude rather than X and Y coordinates?
EDIT:


texts: ArcGIS Tool to Delete only Geodatabase content (feature classes, tables, etc)?

texts:
Situation:

Need to add a Tool to my ModelBuider to Delete only the content from an Enterprise (SDE) Geodatabase (containing features classes, tables, relationship classes).

Tool "Delete" (from Toolbox "Data Management") deletes the Geodatabase and only want to delete its content.

Any help?


texts: ArcGIS for Desktop is rotating attachment pictures from feature service of ArcGIS Online?

texts:
I have a web map on my ArcGIS Online account.

In that map are features (line polygon points) with a picture attached to each feature.

After saving the feature service as a local copy on my desktop the attached images are rotated 90 degrees counter clockwise.  This is when displayed on the info tool or on data driven pages.

I'm using ArcGIS 10.2.1.


texts: Raster copy multiple rasters keeping name but dropping extension

texts:
I have ~100 rasters I need to clip to an area. All rasters are saved as .bil files, and for use in later processing I was trying to make these have 'no extension for esri grid', but keep the same name.
i.e. depth.bil becomes depth, hazard.bil becomes hazard etc. (There are multiple named files in each folder and in this model build example I am calling via %name% on iteration)

I tried to change this in 'calculate value' to take off the extension of the file:

expression = changeName("%Name%")
code block =  def changeName(name): return str(name[:-4])

This worked once, but when I changed folder it put it out of sync so it was naming files incorrectly to the raster it was clipping.
If I don't remove this it says as a .bil file given it needs no extension in the clip tool.
I can do some python but from what I know I would still be in the same predicament as this in ModelBuilder.


texts: Splitting a layer containing multiple polygons into multiple layers each containing one polygon

texts:
So firstly, I used the buffer tool in Qgis3 to create (twenty) polygons with a certain width from (twenty) points.
This lead to the situation below, where o1 and o2 are two original points and 1 and 2 are the resulting polygons.
Now, I want points (like A in the picture below) from another point layer to have an attribute stating in which polygon they are. However, in overlapping situations I would like to see both polygon id's to be added, which does not happen when using the 'add polygon attribute to point' (which only adds one polygon id).

To resolve this, I need to have a different layer for each polygon, so I can use dummy variables. How do I create a different layer for each polygon? I have heard of the Split Polygon plug-in in Qgis2, but I am working in Qgis3. To clarify, I would like for both polygons to be kept intact, as I've seen others wanting the overlapping section belonging to only one of the polygons.


texts: How could I make this PostGIS query run faster?

texts:
SELECT st_band(rast, 1), st_width(rast), st_height(rast),
st_bandpixeltype(rast, 1), st_bandnodatavalue(rast, 1),
st_scalex(rast), st_scaley(rast), st_upperleftx(rast),
st_upperlefty(rast) FROM public.elevation WHERE st_intersects(rast,
st_polygonfromtext('POLYGON((-94.00010416666670210
67.50010416666670210, -93.00010416666670210 67.50010416666670210,
-93.00010416666670210 67.00010416666670210, -94.00010416666670210
67.00010416666670210, -94.00010416666670210 67.50010416666670210))',
4269))

There's a GIST index created over rast column:

CREATE INDEX elevation_rast_gist
  ON elevation
  USING gist
  (st_convexhull(rast) );

This is the 'explain analyze' output

Seq Scan on elevation  (cost=0.00..6.30 rows=1 width=32) (actual
time=2039.622..2614.787 rows=2 loops=1)
   Filter: (((rast)::geometry &&
'0103000020AD1000000100000005000000511BE8B4018057C0511BE8B401E05040511BE8B4014057C0511BE8B401E05040511BE8B4014057C0511BE8B401C05040511BE8B4018057C0511BE8B401C05040511BE8B4018057C0511BE8B401E05040'::geometry)
AND _st_intersects('0103000020AD1000000100000005000000511BE8B4018057C0511BE8B401E05040511BE8B4014057C0511BE8B401E05040511BE8B4014057C0511BE8B401C05040511BE8B4018057C0511BE8B401C05040511BE8B4018057C0511BE8B401E05040'::geometry,
rast, NULL::integer))
 Total runtime: 2615.145 ms

The coordinates used in intersection are not important. The point is I
want to fetch all raster rows that intersect a given window. So:

Is there any other additional index I should create?
Is there any other way to write the query to make it faster?

Many thanks in advance, and best regards,


texts:
I am building an addin that will create a selection of multiple features that fit within the scope of multiple 'where' statements. All fields will be within the same layer, the layer will always be named "ParcelsWslivers" and will always have a field named "acres" and "VPA". Currently, when I run my script it returns the following error.

My code currently is

import arcpy
import os

class CreateLayerOfPotentialSlivers(object):
"""Implementation for Test2_addin.button (Button)"""
def __init__(self):
    self.enabled = True
    self.checked = False
def onClick(self):
    # Get the current map document and the first data frame.
    mxd = arcpy.mapping.MapDocument("CURRENT")
    lyr = arcpy.mapping.Layer("ParcelsWslivers")
    field = 'Sliver'

    # select parcels with acres less than 0.01 acres
    arcpy.SelectLayerByAttribute_management(lyr, "NEW_SELECTION", '[acres] <= 0.01')
    # select parcels with greater than 100,000,000 VPA
    arcpy.SelectLayerByAttribute_management(lyr, "ADD_TO_SELECTION", '"VPA" >= 100000000')

I feel that my issue may be that I'm not being specific with the code about my workspace. However I want it to be conditional so that it does'nt matter which computer this is opened on it will still work within the current document.


texts: Detecting geometry change between two line layers in QGIS

texts:
I have two line layers of the exact same features but I suspect there is a difference in some of geometry between the two.
Is there a way to simply detect change between the two layers?
I am not interested in the attribute change; I simply want to see the areas where there is a difference between the two layers.


texts:
I need to burn a DEM surface to generate an open channel. (all dimensions and slopes are known)
I thought I should use raster calculator to subtract a raster represent channel depths from the DEM raster. but there is only one matter, the land surface is varied and hence I don’t know where the channel side slopes will intersect with the DEM.
Is there a function to find the intersection line between two rasters (channel side slope and DEM)


texts: Calculating Sentinel derived NDVI within a Polygon and exporting to a CSV in Google Earth Engine

texts:
This question uses code based off @Kel Markert https://code.earthengine.google.com/349615d7802d59f677181bef0badad9f
I am attempting to get a maximum monthly NDVI value from 6 small polygons over a number of years from Sentinel 2 in Google Earth Engine and export to CSV. The main difference between my code and his is the use of Sentinel 2 data instead of landsat. When I run my exact code with landsat image collection instead of Sentinel it works, but I cannot figure out why.
When I run it with the Sentinel 2 data I get the error

"FeatureCollection (Error) Error in map(ID=00000000000000000001):
Dictionary.get: Dictionary does not contain key: NDVI."

This leads me to believe that there must be a problem with my band names, but when I inspect them after selecting out only the NDVI band into var NDVI_only, there is only 1 band named "NDVI" just like I would expect...
Link to my code
Link to my feature class
var geometry = ee.FeatureCollection("users/marshallthewolf/valley_bottoms");
print(geometry);

Map.centerObject(geometry);

// Filter by Geo and Growing days
var S2_SR = ee.ImageCollection('COPERNICUS/S2_SR')
            .filterDate('2018-05-01', '2020-10-01') //first BDA year
            .filterBounds(geometry) // filter to ROI
            .map(function(image){return image.clip(geometry)}) //iterate and clips roi over whole collection
            .filterMetadata('CLOUDY_PIXEL_PERCENTAGE',"less_than", 25);
// Add NDVI band
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};
// Apply across whole collection
var S2_NDVI = S2_SR.map(addNDVI);

// Select out only NDVI band
var NDVI_only = ee.ImageCollection(S2_NDVI.select(["NDVI"], ["NDVI"]));
print(NDVI_only)

// Update table and export -----------------------------------------------//
var startDate = ee.Date('2018-05-01'); // set analysis start time
var endDate = ee.Date('2020-10-01'); // set analysis end time

var bandName = ee.Image(NDVI_only.first()).bandNames().get(0);
print(bandName)

// calculate the number of months to process
var nMonths = ee.Number(endDate.difference(startDate,'month')).round();

// get a list of time strings to pass into a dictionary later on
var monList = ee.List(ee.List.sequence(0,nMonths).map(function (n){
  return startDate.advance(n,'month').format('YYYMMdd');
}))
print(monList)

var result = geometry.map(function(feature){
  // map over each month
  var timeSeries = ee.List.sequence(0,nMonths).map(function (n){
    // calculate the offset from startDate
    var ini = startDate.advance(n,'month');
    // advance just one month
    var end = ini.advance(1,'month');
    // filter and reduce
    var data = NDVI_only.filterDate(ini,end).mean().reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: feature.geometry(),
      scale: 1000
    });
    // get the value and check that it has data
    var val = ee.Number(data.get(bandName));
    val = ee.Number(ee.Algorithms.If(val,val,-999));
    // return zonal mean
    return val;
  });
  // create new dictionary with date strings and values
  var timeDict = ee.Dictionary.fromLists(monList,ee.List(timeSeries));
  // return feature with a timeseries property and results
  return feature.set(timeDict);
});

// print to see if it is doing what we expect...
print(result);


texts: Getting data on area where shapefile overlaps raster file using QGIS?

texts: I have attached an image to help explain what I am asking. I want to get the average value of where the heat map crosses over each municipality in the shapefile and export those values into a spreadsheet. The heatmap is a raster image. I use QGIS 3.0.


texts: Aggregating attribute data with one-to-many relationship class in ArcMap?

texts:
I am using ArcMap 10.6.

I am joining a table of public housing developments to a parcel map layer, and there is a one-to-many relationship between the parcels and developments: multiple buildings may occupy any given parcel, and each has its own row of attribute data, such as number of housing units. I'm hoping to aggregate the data from the housing developments to the parcel layer, e.g. so that I can display the sum of all housing units contained in each parcel.

I could do this quite easily with an Excel pivot table or with the R aggregate() function, but if possible I'd like to avoid bringing in a whole new table (i.e. with the aggregated data) and keep the relationship within the geodatabase, which I'm less sure how to do.

In very simplified terms, my data looks something like this:

And I'd ultimately like the parcel layer's attribute table to look like this:

How I might accomplish this with a relationship class/ArcMap relate?

It's not the end of the world if I have to aggregate the data outside of ArcGIS and bring in a separate table, but I feel like that should be avoidable, and I'm curious as to how.


texts: How to (programmatically) determine if a WMS endpoint belongs to a GeoServer?

texts:
Given a WMS endpoint URL (and the ability to make calls to it), what is a reliable way to work out whether it is a GeoServer, ArcGIS Server, or something else? The reason we want to know is due to variable support for things like styles and legends, regardless of what the spec says.

For instance, does an ArcGIS Server WMS endpoint:

always have /MapServer/WMSServer in the URL?
always contain a <esri_wms:GetStyles> element in its GetCapabilities' Capabilities section?

Does a GeoServer WMS endpoint GetCapabilities:

always contain a <Keyword>GEOSERVER</Keyword> in its <KeywordList>
always contain <OnlineResource xlink:type="simple" xlink:href="http://geoserver.sourceforge.net/html/index.php"/>

EDIT

To be clearer about why we want to do this, we run a catalogue that connects to many different servers of different types. Esri and GeoServer treat the width and height parameters to GetLegendGraphic differently (the former treats it as the width of the whole graphic, the latter as the width of each box). Similarly, if you want to change the background color of a legend in GeoServer, you pass ?LEGEND_OPTIONS=bgColor:0055ff while for Esri you pass bgColor:0055ff without the LEGEND_OPTIONS. GeoServer respects LEGEND_OPTIONS=fontColor:..., Esri doesn't (as far as I can tell). And finally, Esri returns a pretty useless response to request=getStyles.

So, all up, it's better if we can make this assessment automatically, before making any detailed requests.


texts: Getting coordinates by clicking on QGIS Canvas with PyQGIS?

texts:
Are there any procedures when working with QgsMapCanvas (pyqgis, qgis) in order to obtain a pair of coordinates by clicking on the map canvas?

Are there any sample Python codes that demonstrate this?


texts:
Recently I upgraded my qgis from 2.6.1 to 2.14.9. In the old version I used the tool 'polygonize' (not the plugin!) to transform lines into polygons. The results were always very good. Sadly, this tool no longer exists in the 2.14 version. Instead there is the 'lines to polygons' tool. However, this new tool does not seem to work properly. I did some tests with line layers which I had succesfully transformed into polygons in v2.6 but in the new version I do not get the same result! See the difference in results below:

I did some research and I see that more people are struggling with the lines to polygons tool. Is there a way to get the old polygonize tool back or is this problem resolved in more recent versions of QGIS?


texts: How to view metadata of ArcGIS World Imagery Clarity?

texts:
I want to see the year and resolution for certain areas for the World Imagery Clarity.
Clarity: https://www.arcgis.com/home/item.html?id=ab399b847323487dba26809bf11ea91a
Now, it seems there is some kind of existing layer, but when I load it into ArcGIS Pro, it doesn't seem to work. Nothing comes up besides the imagery and the street overlay.
Existing (broken) metadata: https://www.arcgis.com/home/item.html?id=298f03f7f02c465a87bcf817c3b8859b
I was hoping for something more like this: https://support.esri.com/en/technical-article/000018129
The above link is for metadata for the World Imagery Classic, not for Clarity. I wish I could get a polygon like that for Clarity. That is my ideal format.
What I did find was ArcGIS Wayback, which is helpful, but it's not exactly what I'm looking for. I'm looking for some kind of polygon metadata later for World Imagery Clarity.
ArcGIS Wayback: https://livingatlas.arcgis.com/wayback


texts: Data update from File Geodatabase to PostGIS

texts:
We are using PostGIS to serve our geographic data by WMTS with GeoServer, and the origin data are produced with ArcMap (using the file geodatabase format).

We use the ArcMap database connection tool to export data to the database.  However all the data take up almost 20GB, it may cost 10 hours to export all of them to the database.

Now once there are some modification of the origin data, how to update the database effectively?

It is not a wise option to export them all again.


texts: Troubleshoot with "Split selected features with selected line from another layer"

texts:
I am trying to use the "Split selected features with selected line from another layer" tool from the plugin Digitizing Tools in QGIS.  Sometimes, this tool works with my shapefiles but sometimes it won't for reasons I can't figure out.  So far, my troubleshooting skills have led me to check:

a) the projection for the project, the polygon and the polyline is all the same

b) the polyline extends beyond the polygon I want to cut

c) if I create a new polyline in my polyline .shp file, I can use it to cut the polygon (but I need to use the pre-existing polylines for cutting) [I tried this after reading the bug report http://hub.qgis.org/issues/15312 ]

I had wanted to add this question as a comment to a previous Question on this tool but my StackExchange standing isn't high enough: Is there a way to split polygons by several lines in QGIS?


texts: projectRaster() brick from NetCDF changes all values to NA

texts:
I have loaded a NetCDF brick like so...

library(raster)
library(rgdal)
library(ncdf4)
nctemp <- brick("C:/nco/out.nc")

Which has the following dimensions...

class       : RasterBrick
dimensions  : 103, 200, 20600, 1680  (nrow, ncol, ncell, nlayers)
resolution  : 0.125, 0.125  (x, y)
extent      : 244, 269, 40, 52.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\nco\out.nc
names       : X1960.01.16, X1960.02.14, X1960.03.15, X1960.04.14,
X1960.05.15, X1960.06.14, X1960.07.15, X1960.08.15, X1960.09.14,
X1960.10.15, X1960.11.14, X1960.12.15, X1961.01.16, X1961.02.14,
X1961.03.16, ...
Date        : 1960-01-16, 2099-12-16 (min, max)
varname     : tmin

Further, print(nctemp) returns the following...

2 variables (excluding dimension variables):
double time_bnds[bnds,time]
double tmin[lon,lat,time]
    long_name: Downscaled Minimum Temperature in Degrees Celsius
    units: degreesC
    _FillValue: 1e+30
    missing_value: 1e+30

4 dimensions:
lon  Size:200
    standard_name: longitude
    long_name: longitude
    units: degrees_east
    axis: X
lat  Size:103
    standard_name: latitude
    long_name: latitude
    units: degrees_north
    axis: Y

Ultimately, I have a shapefile that I want to use to summarize the pixels. I want to project the raster brick to Albers Equal Area, like so...

nc <- projectRaster(nctemp,crs="+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

But when I attempt this or any other transformation, all of the values are changed to NA. I tried to follow along with the suggestions in this question:
https://stackoverflow.com/questions/38356326/projectraster-returns-all-na-values

...and moreover, the linked question on GIS Stack Exchange. However, when I project my shapefile to Albers and then change the extent of the raster to match my shapefile, the projection fails and gives me the error of missing values. Some other possible clues to something weird going on...

it can't rotate (raster::rotate gives me an error saying it doesn't "look like an appropriate object for this function")
a function like raster(inputfile, varname="lat")doesn't work; when I try that,
it says that "lat" is not a variable (though it's a dimension for the actual variable of interest, "tmin" in this case, which seems to work for other intents and purposes).

What am I missing?


texts: making image to geometry with mobis images

texts:
I am trying to do the classification of farming areas and I do not want to include the city and forest areas. so I want to use land without forest and city.
I want to make forest_and_city to geometry so I can use it in my study area.

///////////////////////////////////////////////////////////////
    //// whole world land data
    ///////////////////////////////////////////////////////////////
    var dataset = ee.ImageCollection('MODIS/006/MCD12Q1');
    var igbpLandCover = dataset.select('LC_Type1');

    ///////////////////////////////////////////////////////////////
    //// Cut  bordor and separated image yearly
    ///////////////////////////////////////////////////////////////
    var multiband = igbpLandCover.toBands()
    // print(multiband)

    ///////////////////////////////////////////////////////////////
    //// all propery(17)
    ///////////////////////////////////////////////////////////////
    var property1 = multiband.eq(1)
    var property2 = multiband.eq(2)
    var property3 = multiband.eq(3)
    var property4 = multiband.eq(4)
    var property5 = multiband.eq(5)
    var property6 = multiband.eq(6)
    var property7 = multiband.eq(7)
    var property8 = multiband.eq(8)
    var property9 = multiband.eq(9)
    var property10 = multiband.eq(10)
    var property11 = multiband.eq(11)
    var property12 = multiband.eq(12)
    var property13 = multiband.eq(13)
    var property14 = multiband.eq(14)
    var property15 = multiband.eq(15)
    var property16 = multiband.eq(16)
    var property17 = multiband.eq(17)

    var forest_and_city = property1.max(property2)
    forest_and_city = forest_and_city.max(property3)
    forest_and_city = forest_and_city.max(property4)
    forest_and_city = forest_and_city.max(property5)
    forest_and_city = forest_and_city.max(property13)

    forest_and_city=forest_and_city.select(17)

    forest_and_city= forest_and_city.neq(1).clip(geometry)
    forest_and_city = forest_and_city.updateMask(forest_and_city.neq(0));

    print(forest_and_city)
    Map.addLayer(forest_and_city,{},'forest_and_city')


texts:
I am getting "keyword can't be an expression" error while running the following code:

>>> import arcpy
>>> arcpy.MakeFeatureLayer_management("Point","P1","Layer" = 'WM-VALVE-ID')

I ran "Make Feature Layer" tool from the toolbox which worked. Below is the python snippet:

arcpy.MakeFeatureLayer_management("Point","Point_Layer","/\Layer/\ = 'WM-VALVE-ID'","#","OBJECTID OBJECTID VISIBLE NONE;Shape Shape VISIBLE NONE;Entity Entity VISIBLE NONE;Handle Handle VISIBLE NONE;Layer Layer VISIBLE NONE;LyrFrzn LyrFrzn VISIBLE NONE;LyrLock LyrLock VISIBLE NONE;LyrOn LyrOn VISIBLE NONE;LyrVPFrzn LyrVPFrzn VISIBLE NONE;LyrHandle LyrHandle VISIBLE NONE;Color Color VISIBLE NONE;EntColor EntColor VISIBLE NONE;LyrColor LyrColor VISIBLE NONE;BlkColor BlkColor VISIBLE NONE;Linetype Linetype VISIBLE NONE;EntLinetype EntLinetype VISIBLE NONE;LyrLnType LyrLnType VISIBLE NONE;BlkLinetype BlkLinetype VISIBLE NONE;Elevation Elevation VISIBLE NONE;Thickness Thickness VISIBLE NONE;LineWt LineWt VISIBLE NONE;EntLineWt EntLineWt VISIBLE NONE;LyrLineWt LyrLineWt VISIBLE NONE;BlkLineWt BlkLineWt VISIBLE NONE;RefName RefName VISIBLE NONE;LTScale LTScale VISIBLE NONE;Angle Angle VISIBLE NONE;ExtX ExtX VISIBLE NONE;ExtY ExtY VISIBLE NONE;ExtZ ExtZ VISIBLE NONE;DocName DocName VISIBLE NONE;DocPath DocPath VISIBLE NONE;DocType DocType VISIBLE NONE;DocVer DocVer VISIBLE NONE;ScaleX ScaleX VISIBLE NONE;ScaleY ScaleY VISIBLE NONE;ScaleZ ScaleZ VISIBLE NONE;ID ID VISIBLE NONE;ID_ ID_ VISIBLE NONE;B_P B_P VISIBLE NONE")

So, when I copied the Python snippet into my code and tried running it, I got another error "ERROR 000230: Failed selecting with /\Layer/\ = 'WM-VALVE-ID' Failed to execute (MakeFeatureLayer)."


texts: SQLite error - no such function: InitSpatialMetadata

texts:
I am porting a QGIS plugin to QGIS3.

I get a problem when I try to create a new SpatiaLite DB, in particular when I initialize spatial metadata (GEOMETRY_COLUMNS and SPATIAL_REF_SYS).

My code is:

import sqlite3
con = sqlite3.connect(dbname)
con.enable_load_extension(True)
cur = con.cursor()
cur.execute("PRAGMA foreign_keys = ON")
cur.execute("SELECT InitSpatialMetaData(1);")

I get the following error :

  "no such function: InitSpatialMetadata"

My SQLite version is 3.29.0.

Any idea about this?


texts: GDAL/OGR Determine if dataset is raster or vector

texts:
I am writing a script using the GDAL/OGR Python bindings and it is possible that the user may provide a raster dataset or a vector dataset as input.  If the user provides a vector dataset, the script will convert it to raster.

Is there a way to determine whether the dataset is raster or vector before opening them?

I could just try to open the dataset with OGR first and if it returns None I could try to open with GDAL, but I am looking for a more simple way like gdal.IsRaster(filepath) if it exists.


texts: QGIS 2.18.4 having errors because of my preexisting Python 2.7.12 installation

texts:
I am an experienced Python developer, but new to QGIS. I just installed QGIS 2.18.4 and am getting the following error message when I start it.

Couldn't load QGIS utils.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line 21, in
    standard_library.install_aliases()
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\future\standard_library\__init__.py", line 465, in install_aliases
    from future.backports.urllib import request
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\future\backports\urllib\request.py", line 93, in
    import hashlib
  File "C:\Python27\Lib\hashlib.py", line 138, in
    _hashlib.openssl_md_meth_names)
AttributeError: 'module' object has no attribute 'openssl_md_meth_names'

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.18.4 'Las Palmas', 59c1d21

Python path:
['C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python',
u'C:/Users/462974/.qgis2/python',
u'C:/Users/462974/.qgis2/python/plugins',
'C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\matplotlib-1.3.1-py2.7-win-amd64.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\nose-1.3.3-py2.7.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\tornado-4.0.1-py2.7-win-amd64.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\backports.ssl_match_hostname-3.4.0.2-py2.7.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\certifi-14.05.14-py2.7.egg',

'C:\\Python27\\Lib',
[Rest of my preexisting PYTHONPATH],

'C:\\PROGRA~1\\QGIS2~1.18\\bin\\python27.zip',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\DLLs',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\plat-win',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\lib-tk',
'C:\\PROGRA~1\\QGIS2~1.18\\bin',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\PIL',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\jinja2-2.7.2-py2.7.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\markupsafe-0.23-py2.7-win-amd64.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\pytz-2012j-py2.7.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\win32',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\win32\\lib',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\Pythonwin',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg',
'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg']

I already had Python 2.7.12 installed on my machine before installing QGIS; it seems to have come with a prepackaged version of Python, 2.7.5. It appears to be modifying my PYTHONPATH environment variable, adding new directories to it on both sides. The problem is that because it adds C:\Program Files\QGIS 2.18\apps\Python27\lib to the end of the path, after C:\Python27\Lib, it finds the Python 2.7.12 version of hashlib first, which apparently contains code the rest of QGIS' Python distribution can't handle (like the line causing this AttributeError). I scanned the Settings menu and don't see anything about how the PYTHONPATH is constructed. Is it normal for QGIS to not play nice with preexisting Python installations?


texts: Clipping with Voronoi polygon using QGIS

texts: I'm trying to clip a road network (polyline) with Voronoi polygons. The polyline don't clip but I don't know why. I try to separate my road network according to geographical points collected in the field. This seems to be the only method I know of to separate road sections based on the value of a point.


texts: Openlayers, add new ol.layer.Tile after AJAX 'done'

texts:
I'm trying to add an additional ol.layer.tile on my existing map view after a 'status' 'done' Ajax post call. I have tried to use "async" and "await" on initialize the ol.map. A sample code is below:

async function MapUI()
{

// Create Base Map and Layer using OpenLayers

  var view = new ol.View
  ({
  // change projection to geoserver's projection

    projection: "EPSG:4326",
    center: [-98.5284468774974,52.71391708013823],
          zoom: 7
  });

  // creating a openlayer base map

  var osm_map = await new ol.layer.Tile
  ({
    source: new ol.source.OSM()
  });

  // creating a map layer from geoserver

  var vector_source = await new ol.source.Vector
  ({
        url: MyGeoJsonURL,
        format: new ol.format.GeoJSON()
  });

  var vector_map = await new ol.layer.Vector
  ({
    source: vector_source

  });

  var raster_source = await new ol.source.TileWMS
  ({
    url:'http://localhost:8080/geoserver/Canola/ows?',
    params:{'LAYERS': 'Canola:Riskmap'},
    serverType: 'geoserver',
    crossOrigin:'anonymous'
  });

  var wms_map = await new ol.layer.Tile({

    source:raster_source
  });

  var map = await new ol.Map({

    layers: [osm_map, vector_map,wms_map],
    target: 'map',
    view: view
    });
}

Or the after "done" do something method inside of the call, like shown below

$.post("../lib/php/riskAnalysis.php",
{
     CropDensity:CropDensity,
     DiseaseHistory:DiseaseHistory,
     RegionRisk:RegionRisk,
})
.done(function(databack, status)
{
    $("#status").html("done");
    {
 var raster_source = new ol.source.TileWMS
  ({
    url:'http://localhost:8080/geoserver/Canola/ows?',
    params:{'LAYERS': 'Canola:Riskmap'},
    serverType: 'geoserver',
    crossOrigin:'anonymous'
  });

  var wms_map = new ol.layer.Tile({

    source:raster_source
  });

  map.addLayer(wms_map); // where map is the original map layer

}

})

None of the above code works... Not sure did I have any errors in my code or my approach was wrong.


texts: Are there sources for obtaining pre-made KML data from Google?

texts:
I have been asked to investigate methods for obtaining KML files.

All I know so far is that it is a format based on XML that Google uses to format their maps (or the metadata).

However, are there places or methods that I may avail myself of to download KML data without going through Google Earth?

For example, is it possible to get the KML for the vector data of a state in the US (roads, boundaries, etc.) ?

To restate it simply: I want Google's vector data so I can display it in ArcMap. Are there repositories for such things?


texts: rasterize a geography polygon using postgis 2.1

texts: I am still bit new to Postgis 2.1. I am trying to rasterize a polygon, meaning that I have a polygon of Syria and I want to create an empty raster that contains the shape and geography reference of syria. I have been looking at the ST_AsRaster() function, but I can't seem to figure out the settings for this. Any suggestions.


texts: Using Geotiffwrite or Worldfilewrite for global TIF in MATLAB?

texts:
In MATLAB I have written the code:

[image,geo]=geotiffread('africa.tif');
info=geotiffinfo('africa.tif');
geotiffwrite('Wafrica.tif',image,geo);

but there was error:

geotiffwrite>validateR (line 847) The input, R, is a
map.rasterref.MapCellsReference object indicating that you are
working in a projected coordinate system. If so, then specify a
projected coordinate system by setting the appropriate values for the
'CoordRefSysCode' or 'GeoKeyDirectoryTag' optional parameters.
geotiffwrite>validateInputs (line 348) R = validateR(R, size(A),
hasColorMap, type);
geotiffwrite (line 238) [filename, A, cmap, R, Params] =
validateInputs(filename, A, cmap, R, Params);


texts: Calculating latitude from X,Y-coordinates on Mercator projection

texts:
I have a task to create a map of a specified section of the world. In particular, I need to map everything between the 120E and 100W longitudes and everything between the 60S and 30N latitudes mapped (I am making a map of the Pacific). I will be using the Mercator projection (since I do not have to map areas of extreme latitude).
Suppose that my map will be A cm wide and B cm high.
Trivially, I can calculate the longitude from the x-axis coordinates by using the following formula:
Longlitude = ((x/map_width) * 220 ) - 120

(x=0 is at the far left of the map, corresponding to the 120E longitude, which is why I take away 120 at the end)
Now, I require a similar formula for calculating latitude from the y-axis coordinates, with y=0 corresponding to the 30N latitude
So far, I have used Xarinko's answer from https://stackoverflow.com/questions/2103924/mercator-longitude-and-latitude-calculations-to-x-and-y-on-a-cropped-map-of-the to create this function in Python:
map_width = 1000 #in mm
map_height = 700 #in mm

map_long_left = -120
map_long_right = 100
map_long_delta = map_long_right - map_long_left

map_lat_bottom = -60
map_lat_bottom_rad = map_lat_bottom * math.pi / 180

map_lat_top = 30
map_lat_top_rad = map_lat_top * math.pi / 180

def xy_to_longlat(x: float, y: float) -> tuple():
    world_map_radius = (map_width / map_long_delta) * (360/(2 * math.pi))
    map_offset_y = ( world_map_radius / 2 * math.log( (1 + math.sin(map_lat_bottom_rad) ) /
                                                    (1 - math.sin(map_lat_bottom_rad))))
    equator_y = map_height + map_offset_y
    a = (equator_y - y) / world_map_radius

    lat = 180 / math.pi * (2 * math.atan(math.exp(a)) - math.pi/2)
    long = map_long_left + (x / map_width) * map_long_delta
    return (long, lat)

This function assumes that the map is from the 60S to 60N latitudes. But I only want to map between the 60S and 30N latitudes. It also does not accurately calculate latitude. It produces a close value, but it can be significantly off in some instances.


texts: How should the bounds of decimal and degree coordinates be programmatically handled in respect to hemispheres?

texts:
Are the bounds inclusive of 0, 90/180?

Many systems and libraries simply clamp it to from -n to +n without considering zero. In many of these cases the handling of 0 is arbitrary (which hemisphere it lands on for example). In other cases -0 is supported (which offers more even distribution).

Some descriptions suggest that 0 doesn't really exist, values simply move closer to it but never reach it, however some implementations include 0.

If it's down to the person implementing the system which system is preferred?

0 < n < max
0 <= n < max
0 < n <= max
0 <= n <= max

Given the even options should 0 be negative or fall on a specific hemisphere?

I can't answer my own question as it's locked...

https://en.wikipedia.org/wiki/Signed_zero

From a programming perspective, standards wise signed zero is an anomaly, it is an artifact that is often ignored and not supported directly by programming languages. It's also mathematically broken. The different between -0 and 0 is still 0.

An algorithm wont care about human representation and segments.

Given a coordinate in the range -90 to 90, it will most likely have 90 added to it and then possible another 180 depending on hemisphere.

Thus simplifying things to two values of 0..360 (an internal representation simplifed as much as possible, removing the sign entirely would be preferable for machine use) starting at point 0. As far as the machine is concerned there are no such thing as hemispheres, it is a presentation concern.

This does still become a problem for humans, where which hemisphere 0 lands on is still arbitrary.

Problems can also arise where machines talk to one another using coordinate systems then contain a format intended for humans (they specify hemisphere).

However it's very unlikely for something to land exactly on 0,0 and the precision loss is usually so small that it very rarely, almost never, manifests as a real problem. Many systems will add some extra precision to either eliminate the problem by isolating it to within a tiny fragment of least significance or to slow it's accumulation.

Whatever method is used involving hemispheres, it should be well defined enough to convert both ways without precision loss (such as > being used for output but < on input).

Hemisphere replaces the sign. As it makes no sense for coordinates of exactly zero, the obvious option is to leave it out. As with the floating point representation, it's an artifact that can be ignored.


texts:
How can I add a legend to a raster to the Globcover land cover dataset in QGIS without plugins?
The data can be found here http://due.esrin.esa.int/globcover/
Unfortunately, as far as I know, QGIS doesn't support the use of .lyr files as legend for classified raster datasets. For an import dataset I created a QGIS compatile legend file.


texts: Is it possible get pixel value with only PyQGIS and not GDAL?

texts:
It seems that gdal can caused in some crashes in QGIS, specific in used for raster layers (such as links : QGIS crashes on close if data was read through GDAL in Python Console  or GDAL/QGIS broken Ubuntu 14.04 or Why close QGIS for extract pixel values by GDAL?.
So Is it possible to get pixel values only with pyqgis and without gdal?


texts: How to add z value vertices from a point layer to a line in QGIS?

texts:
I created a point layer from a line using Vector geometry/Extract vertices in the Processing toolbox.

The purpose is to extract the elevation from a DEM into the points corresponding the vertices and then into the line layer.

I managed to get the elevation values on the points but I can't really figure it out how to add the values to the vertices in the line layer.

My final objective is to export a LineString geojson file containing the x, y, z values.
There might be easier ways but I couldn't find any so far.

P.S. I created a z value field in the line layer, just need to copy the values from one layer to another.


texts: How to modify two polygons which have common point and line at the one time

texts:
I would like to modify two polygons at the one time because these two data have same points and lines. Is it possible?

In the present circumstance, I have to modify polygon one by one. First, I modify one polygon. After that, another polygon will be modified with snap function for overlapping to another polygon.


texts:
How do I define a custom CRS in QGIS from the Mapinfo grid definition below?

"MIM Grid",1008,33,7,139.48639389,-20.69008981, 1.000063,0,0,7,0.998796124,-0.049054077,12520.6,0.049054077,0.998796124,39952.907
For your information the Definition is:
Datum: Geocentric Datum of Australia (GDA94)
Projection: Universal Transverse Mercator
Central Meridian: 139.48639389 (decimal degrees) 139d 29’11.01802” E (degrees minutes seconds)
Origin Latitude: -20.69008981 (decimal degrees) 20d 41’24.32332” S (degrees minutes seconds)
Central Scale factor: 1.000063
Rotation Matrix: [ cos (theta)= 0.998796124 -sin (theta)=-0.049054077
sin (theta)= 0.049054077 cos (theta)= 0.998796124 ]
X offset: 12520.6
Y offset: 39952.907

I had a look at the QGIS documentation and PROJ4 but can't work out a way to input this information in a CRS definition in QGIS.


texts: Understanding why list object has no attribute name?

texts:
I have this script which gives me an attribute error and points out line 8.

The error says

  " 'list' object has no attribute 'name'"

Not sure what went wrong.

How do I fix this script to run?

import arcpy
from arcpy import env
env.workspace = "C:/EsriPress/Python/Data/MyFolder"
fc = "airports.shp"
rows = arcpy.SearchCursor(fc)
fields = arcpy.ListFields(fc)
for field in fields:
    if fields.name == "NAME":
        for row in rows:
            print "Name = ({0}".format(row.getValue(field.name))


texts: How to determine which tiles cover your Area of Interest?

texts:
I'm relatively new to using a GIS.

I've been wondering how could I determine tile(s) on a selected area/point at a specified zoom level?


texts: From a shapefile, copy specific rows (defined by a separate table) to a new shapefile

texts:
I have a polygon shapefile of origin zones.
I also have a polygon shapefile of destination zones.
I also have a table that lists origin-destination pairs.

I would like to create a multi-part shapefile where each row represents an origin-destination pair. So when you select a single row, two separate polygons would be highlighted (one polygon representing the origin zone, the other representing the destination zone).

I don't know how to go about creating the multi-part shapefile. Does anyone have any suggestions?


texts: How can I get the elevation of a point in openlayers?

texts:
I want to know if it is possible to obtain the elevation of a point with some type of control of others. I am painting a linestring with a series of coordinates but I need to know the elevation of every point in order to print a chart. I know that Google has an API for elevation, but I need to use only openlayers.

Note: My baselayer is a google map.

thanks!


texts: Modify sidebar information caused by marker click

texts:
First time poster here, hopefully someone can point me in the right direction! Pretty new to Leaflet, HTML, and CSS in general so this project has been a great learning experience for me but I've really hit a roadblock at this point.

I'm using Leaflet to create a WebMap of a walking tour for a local town. http://beattyre.github.io/WebMapTest/ I have pictures, points, and descriptions stored in a GeoJSON file and have created a webpage for testing purposes as I work on this project. Basically, what I hope to do, is to have the sidebar information to be populated with the GeoJSON attributes when I click on a point, rather than create a pop-up with the information. I've seen it done before on other webmaps such as this one http://benjmyers.com/ratmap/#/, but am having trouble figuring out how it's done. I think I've discovered that AngularJS is the library is what I need to accomplish this, but I'm not quite sure how. Can anyone point me in the right direction or give me some tips on how to create this feature?


texts: line from min to max of a digital elevation model

texts:
I have some geometries as multipolygons and a dem (digital elevation model).
I need a line from the max value to the lowest value of the dem for every multipolygon.
Is qgis a tool i could use or do i need PostGIS ?
How could i realize this problem ?


texts: Proper handling in PostgreSQL/PostGIS only after many queries?

texts:
I am trying to analyze GPS data of a parade on a fixed route with a certain variance due to the urban surrounding (canyoning). I'd like to project the GPS waypoints to the nearest points on the route underlying linestrings.
I add a new column in the same table in which the data is located and update this with the calculated data. But during the first few queries I get for all projected waypoints the projection of the first point.
Only if I repeat the query a few times the desired output is generated.
Because I'm a beginner in programming I can not figure out what seems to be the trouble.

SELECT AddGeometryColumn ('vehicle_2015_02','projected_point',4326,'POINT',2, false);
UPDATE vehicle_2015_02 v02
SET projected_point = ST_LineInterpolatePoint(
    romo2015_way.st_makeline,
    ST_LineLocatePoint(
        romo2015_way.st_makeline,
        vehicle_2015_02.geom
    ))
FROM romo2015_way, vehicle_2015_02;

vehicle_2015_02:      table with GPS waypoints+timestamp as points
projected_point:      new column in table "vehicle_2015_02" as points
romo2015_way:         table with route as points


texts: Filter features by geometry in OpenLayers 2

texts: I have three different type of features on my layers: points, lines, and polygons. I need to create a rule based style based on these three types of layers. So, for example, all points could be green, lines red, and polygons blue. I've been looking through the docs but I can't seem to really find anything that would allow OL2 to do this. Can anyone provide an example, or point me in the right direction?


texts: Importing and viewing polygon layer on map?

texts:
How is it possible to import and view a polygon layer (I have it in a shapefile) on a map inside MIT App Inventor? Maybe with LoadFromURL method, but in this case which file I should paste associated with the URL?

Usuful resources I found are:

https://developers.google.com/maps/documentation/javascript/importing_data
https://developers-dot-devsite-v2-prod.appspot.com/maps/documentation/javascript/examples/layer-kml
https://puravidaapps.com/geolocation.php

I don't know much about GeoJSON or JavaScript API, so: now that I have a map in the app, how should I proceed?


texts: We have ArcGIS Enterprise in SQL and I see a lot of iNumber tables. I can't find information about what is the purpose, how they got created and do I need to delete them?


texts: Unsuccessfully trying to iteratively run gdal_polygonize in a batch file

texts:
Thanks in advance for your help with this.

I've been reading, experimenting, and banging my head against the wall for a couple days on this, and could really use some assistance. This is my first batch file attempt.

I'm trying to loop through .tif files in a folder and run a gdal process on them. I have been calling the .bat file via command line in the OsGeo4W terminal.

I can successfully run the gdal_polygonize process on individual files, but have not succeeded in running it iteratively.

SETLOCAL EnableDelayedExpansion

SET MYPATH = C:\Users\mkcarte2\Desktop\Polygonize\

FOR /F %%i IN ('DIR /B %MYPATH%*.tif') DO (

    SET MYPATHFILE=%%~nI ::Trying to parse only the filename, so that I don't stack extensions in the output file
    SET infile=%%i
    SET outfile=%MYPATHFILE%.shp!

    gdal_polygonize %MYPATH%!infile! -F  "ESRI Shapefile" %MYPATH%!outfile!
)

My Error Messages:

For first file:

  Warning 1: Layer name 'ESRI Shapefile' adjusted to >'ESRIShapefile'
  for XML validity. Warning 1: Field name '%~nI.shp' adjusted to
  '_nI.shp' to be a >valid XML element name.
  0...10...20...30...40...50...60...70...80...90...100 - done.

For subsequent files:

  ERROR 4: `!infile!' does not exist in the file system, and is not
  recognised as a supported dataset name. Unable to open !infile!


texts:
I've become a big fan of a visualisation technique called grid maps. Here are some examples:

http://flowingdata.com/2015/05/12/the-great-grid-map-debate-of-2015/
http://i.telegraph.co.uk/multimedia/archive/03297/2010-vs-2015_3297055b.gif
https://flowingdata.com/2016/09/26/the-spread-of-obesity/

There is a tutorial on automating the creation of these using D3.js and force directed layouts but I'm not skilled in this area (plus I'd feel more comfortable handling a shapefile than I would SVG, even if it is grossly distorted in terms of geographical accuracy):

https://medium.freecodecamp.com/creating-grid-map-for-thailand-397b53a4ecf#.y7yf7m61f

Are there any QGIS / R / Python techniques which could be used to convert a shapefile to a grid map based on evenly sized squares or hexagons (whilst retaining each feature's data)?

I realise that this would completely destroy/distort the geographic accuracy of the shapefile, my interest is being able to generate these for use as high level visualisations using ggplot in R.

Edit 1: With a single tag edit, iant opened my eyes to googling for 'cartogram' (d'oh) which has yielded this site's consideration on the issue: https://richardbrath.wordpress.com/2015/10/15/equal-area-cartograms-and-multivariate-labels/

Also, I thought I had duplicated another post but their question is a bit different as they want an equal number of a measure to fall within each feature: How to create contiguous area cartograms in QGIS or R?.

Edit 2: In response to stev-k's answer, here are the results from running the suggested script. Red points are centroids of the polygons (in grey) with the hexagons overlayed to show how some are empty and some are crowded.


texts: GRASS r.viewshed Command Line Error: Raster map  not found

texts:
I'm struggling to make sense of the grass command line syntax for specifying a raster geotif as an input parameter for a watershed analysis using a DEM tif.  I may also have the user folder setup, database, mapsets completely wrong too.
The commands are to be run from the CMD in windows 10 using GRASS 7.6.
My command:
grass76 --text C:\Temp\grass\newLocation\PERMANENT --exec r.watershed -s elevation="C:\Temp\grass\database\dem.tif" threshold=10000

Message:
Starting GRASS GIS...
WARNING: Concurrent mapset locking is not supported on Windows
Executing <r.watershed -s elevation=C:\Temp\grass\database\dem.tif threshold=10000> ...
SECTION 1a (of 5): Initiating Memory.
ERROR: Raster map <C:\Temp\grass\database\dem.tif> not found
WARNING: Subprocess failed with exit code 1
Execution of <r.watershed -s elevation=C:\Temp\grass\database\dem.tif threshold=10000> finished.
Cleaning up temporary files...
Press any key to continue . . .

I've tried without " around the filepath, with forward slashes, within <>.
I've copied the dem.tif to all the folders: 'C:\Temp\grass\newLocation', 'C:\Temp\grass\newLocation\command', 'C:\Temp\grass\newLocation\PERMANENT', 'C:\Temp\grass\database'
Any obvious syntax problems?
It seems like the --exec is successfully recognising I want to run a command put claims there is a raster file path or file access issue...
I want to avoid any setting up of a GRASS GIS Database directory, GRASS Location or a Mapset...whatever that is.
I just want to be able to run GRASS commands from the command prompt - in reality probably automated via python and subprocess and to change to input parameters such as the data source and output files.  So I would prefer to do that via within the command itself.  I want to avoid user input of setting up an environment.
Is that possible?
Have I mucked up the setup which is causing issues? (I set EPSG of 27700 when creating the user folder).


texts:
I have a data set of points from which I'm trying to make a selection of the points whose location is between two line segments. What would the way to do it?

Here's a picture to illustrate what I'm trying to get:


texts:
i have a college project that searches content in map. I did  all instruction written in this site:
http://kartoza.com/how-to-build-and-debug-qgis-with-qtcreator

then built my sample project.
this is my .pro file:

QT+= core gui xml

greaterThan(QT_MAJOR_VERSION, 4): QT += widgets

TARGET = test
TEMPLATE =lib

SOURCES += main.cpp\
    mainwindow.cpp

HEADERS  += mainwindow.h

FORMS    += mainwindow.ui

INCLUDEPATH +=/home/amin/Documents/qgis/QGIS-master/src/core
INCLUDEPATH +=/home/amin/Documents/qgis/QGIS-master/src/core/symbology-ng
INCLUDEPATH +=/home/amin/Documents/qgis/QGIS-master/src/gui

DEFINES += GUI_EXPORT= CORE_EXPORT=

and this here is my source file:

QgsApplication app(argc, argv,true);

 QString myPluginsDir        = "/home/timlinux/apps/lib/qgis";
 QString myLayerPath         = "/home/timlinux/gisdata/brazil/BR_Cidades/";
 QString myLayerBaseName     = "Brasil_Cap";
 QString myProviderName      = "ogr";
 // Instantiate Provider Registry
 QgsProviderRegistry::instance(myPluginsDir);
 // create a maplayer instance
 QgsVectorLayer * mypLayer =new QgsVectorLayer(myLayerPath, myLayerBaseName,       myProviderName);
 QgsSymbolV2 * x= QgsSymbolV2::defaultSymbol(mypLayer->geometryType());
 QgsSingleSymbolRendererV2 *mypRenderer = new QgsSingleSymbolRendererV2(x);
 QList <QgsMapCanvasLayer> myLayerSet;
 mypLayer->setRendererV2(mypRenderer);
 if (mypLayer->isValid())
 {
    qDebug("Layer is valid");
 }
 else
 {
   qDebug("Layer is NOT valid");
 }

  // Add the Vector Layer to the Layer Registry
 QgsMapLayerRegistry::instance()->addMapLayer(mypLayer, TRUE);
 // Add the Layer to the Layer Set
 myLayerSet.append(QgsMapCanvasLayer(mypLayer, TRUE));

// Create the Map Canvas
QgsMapCanvas * mypMapCanvas = new QgsMapCanvas(0, 0);
mypMapCanvas->setExtent(mypLayer->extent());
mypMapCanvas->enableAntiAliasing(true);
mypMapCanvas->setCanvasColor(QColor(255, 255, 255));
mypMapCanvas->freeze(false);
// Set the Map Canvas Layer Set
mypMapCanvas->setLayerSet(myLayerSet);
mypMapCanvas->setVisible(true);
mypMapCanvas->refresh();

// Start the Application Event Loop
return app.exec();

it built successfully!
I solved all the built issue,but I dont know how to run it and set the run environment run path?


texts: Giving toggle functionality to custom control in OpenLayers

texts:
I created a custom OpenLayers 4 control named "quitarLayers" that removes two layers when I click on it.

How can I make it toggle so that it adds the two layers or removes them, depending if the layers are visible or not?

quitarLayers = function(opt_options) {

var button = document.createElement('button');
button.innerHTML = '<img class="quitalayers" src="img/layerwhite.png"
alt="image">';

var queHaceQuitarLayers = function() {

  map.getLayers().remove(layerPrincipal);
  map.getLayers().remove(layerSecundario);

};

button.addEventListener('click', queHaceQuitarLayers, false);

var element = document.createElement('div');
element.className = 'remove-layers ol-unselectable ol-control';
element.appendChild(button);

ol.control.Control.call(this, {
  element: element,
});

 };
ol.inherits(quitarLayers, ol.control.Control);


texts: Setting Instance of MULTIPLE LAYERS using "mapLayersByName" in QGIS Python Scripting?

texts:
I'm currently using "mapLayersByName" successfully to set one active layer at a time as follows...

vl = QgsMapLayerRegistry.instance().mapLayersByName('My_Layer')[0]
iface.setActiveLayer(vl)

Can an instance of multiple layers be set using the "mapLayersByName" method?

I have tried...

vl = QgsMapLayerRegistry.instance().mapLayersByName('My_First_Layer', 'My_Second_Layer')[0]
iface.setActiveLayer(vl)

As you can see my end goal in this case would be to set the layers within that instance of MULTIPLE layers (vl) as the active layers using the "iface.setActiveLayer()" method.  I notice this method is not plural (i.e. not "setActiveLayers()").  Significant?


texts: Map is missing when I "export as image" via print composer

texts:
I have point data on a vector layer that prints accurately via print composer (when I
"export as image"), however my map below either does not show up or shows up at an incorrect scale.   I used openlayers plugin yahoo street map for the map.

I have examples to illustrate what is happening, but can't post them here since this is my first time.  I'd be happy to email them to someone if you have ideas as to what might be going on.

I'm new to the world of QGIS, so any help on this would be greatly appreciated.  Thanks!


texts: How can I identify lines that intersect a certain raster layer area?

texts: How do I identify certain roads ( of a line shapefile) that intersect with a layer that is produced from the reclassify tool (for raster data)? I am trying to identify which roads specifically are in the same area as the one created by the reclassify.


texts: Segmentation - problem with mosaic tiles

texts:
Could you please tell me how to solve my problem with segmentation in ERDAS ?

I am using ERDAS 2013 and the tool Image Segmentation to segment my multiband picture (with high resolution).

Unfortunately, the results show the mosaic tiles, that come as a result of incorrect color matching in the mosiac dataset.

I do not have the original raster catalog to redo the color matching.
How should I perform my segmentation to get rid of the visible tiles on the segmentation reult?

See the image:


texts: LeafletJS How to set custom map scale for a flat image map (CRS.Simple)

texts:
I've been trying to get a flat image fantasy map working within LeafletJS. I have pretty much everything working well except for the map scale legend where two things are not working:

I can't seem to change the scale of the control to reflect the scale of the map. The scale is trying to set itself for earth LatLng scales yet the map is of a very localised place.
As the map is repositioned in the window the scale changes to reflect the curvature of the earth (so the scale increases the closer you get to the poles).

All I want is a flat scale across all areas of the map and wanting to have the scale accurately reflect the scale in the image. For example, if in the image 1km = 100 pixels how would I go about making that work within the map?

My basic code is as follows:

// ---- VARIABLES ---- //
// Map bounds (pixels)
var mapSW = [0, 8192]; // Deepest Z/X/Y.png +1 x 256.
var mapNE = [8192, 0]; // Deepest Z/X +1 x 256.

// ---- MAP SETUP ---- //
// Setup the map view
// Initial coordinates as well as starting zoom factor

var map = L.map('map').setView([0, 0], 2);

// Reference the map tiles output from Photoshop with properties
// for zoom factors and attribution.
//
// Use CRS.Simple so that the map can be referenced in terms of
// pixels rather than latitude and longitude

L.tileLayer('map/{z}/{x}/{y}.png', {
    minZoom: 2,
    maxZoom: 5,
    nowrap: true,
    crs: L.CRS.Simple,
}).addTo(map);

// Set the map bounds

map.setMaxBounds(new L.LatLngBounds(
      map.unproject(mapSW, map.getMaxZoom()),
      map.unproject(mapNE, map.getMaxZoom())
));

// Scale legend
// TODO. Not working properly yet

L.control.scale().addTo(map);

EDIT / UPDATE:

I don't consider this a valid answer because this is just a hack over existing code. I know, I know - I'm lazy. ;)

Rather than just use Leaflet's Scale control I hacked this into the nicer looking graphicscale add-on (https://github.com/nerik/leaflet-graphicscale).

Added in a couple of options (as different maps will require different scales). Sorry, not sure why the code isn't wanting to format.

L.Control.GraphicScale = L.Control.extend({
  options: {
    ...,
    unitsPer1000px: 1000, // x units per 1000 pixels at max zoom. Metres is default
    scaleUnit: null, // override with your own unit designation. null is default and will revert to m / km
},
...
Replaced the _update method to just use the zoom factor to tweak the scale set in the options.

_update: function () {
    var mapScale = Math.pow(2, this._map.getZoom()) / Math.pow(2, this._map.getMaxZoom()),
        dist     = this.options.unitsPer1000px / mapScale;

this._updateScale(dist, this.options);

},
Replaced the _getDisplayUnit method to allow custom units. If custom units not specified revert to m / km.

_getDisplayUnit: function(units) {
    if(this.options.scaleUnit) {
        // Defined scale unit
         return {
            unit: this.options.scaleUnit,
            amount: units
        };
    } else {
        // Not defined. Use metres/km
        var displayUnit = (units<1000) ? 'm' : 'km';
        return {
            unit: displayUnit,
            amount: (displayUnit === 'km') ? units / 1000 : units
        };
    }
}
Then just call it with the options in the javascript.

var graphicScale = L.control.graphicScale({
    doubleLine: true,
    fill: 'hollow',
    showSubunits: true,
    unitsPer1000px: 6.8,
    scaleUnit: 'parsecs',
}).addTo(map);


texts: Move feature(s) causes points to be duplicated in QGIS 2.8.9

texts:
I have a point shapefile in QGIS where I am manually adjusting some of the points due to overlap. I am using Toggle editing --> Move feature(s) --> Save for all layers. Instead of moving the points these are duplicated. The original point stays where it is, in addition an extra point is created wherever I move it to. This does not seem to lead to any any additional entries in the attribute table.

Any ideas about what I might be doing wrong?


texts: Calculating correlation between points and raster values with missing points (need missing as 0s) (QGIS)

texts:
I feel like this should be easy but I'm struggling. I have a raster data set of Topographic Wetness Index (TWI) and a point data set of crop yield. I want to see how TWI and yield are correlated by point sampling the raster, but I have holes of missing yield monitor data because the field flooded and the crops died in these spots. So, these missing areas are important and should be filled as zeroes, and not interpolated from surrounding point data. Here are some screenshots using QGIS to give a better view of the problem.
My yield monitor data with missing areas from flooding.

My topographic wetness index.

And the yield overlaid on TWI.

How should I fill these missing areas with zeroes so I can point sample and calculate correlation?


texts: Syntax to load wkt or (x,y) formatted geometry from .csv files on PyQGIS

texts:
I'm trying to load some points through a vectorlayer on my QMainWindow, part of my pyqgis application. The code is:

    uri = 'file:///C:/Users/.../data/Targets_com.csv?delimiter=,&crs=epsg:4326&xField=x&yField=y'
    # uri = 'file:///C:/Users/.../data/Targets_wkt.csv?delimiter=,&crs=epsg:4326'#&wktField=wkt'
    # Loading the layer and check it:
    self.target_layer = QgsVectorLayer(uri, 'Targets' , 'delimitedtext')
    if self.target_layer.isValid():
        QgsMapLayerRegistry.instance().addMapLayer(self.target_layer)
        print('Layer OKKK')
    else:
        print('Layer Fail')
    self.target_layer.updateExtents()

When I run my application with that, I get "Layer Fail", so for sure there is something wrong in my "uri" definition, isn't it?

If I change uri like:

uri = 'file:///C:/Users/.../data/Targets_com.csv?delimiter=,&crs=epsg:4326'

that is, cutting the geometry definition (on both cases: (x,y) and wkt), I get the "Layer OKKK", but I can't see my points on the map.

I think the .csv files are well written, cause I can load them on QGIS. Anyway, this are my files.csv (and I tried with one more using ";" like delimiter):

I have read all the posts I found about csv files and syntax with qgis, but I don't know how to solve this, anyone have any idea? i'm using the OSGeo4W package (QGIS 2.16, its own python interpreter, PyQt4,...) and using eclipse-PyDev for development.


texts: Red lines appearing in QGIS Georeferencing?

texts:
I am georeferencing some old aerial pictures in QGIS. The area encompass mountains, hills and a flat area.

For some pictures the error is quite low and I am satisfied with the overlapping (error is around 2 using the polynomial 2). For others, the error is very high (as much as 100!). What bothers me most, is that I cannot understand the red lines that appear in the screen. When I delate these points with long lines, sometimes the error decreases. Other times the nearby points increase even more the lengths of these lines.

I am using as SR the 3003 and as method for transformation polynomial 1 or 2 or 3 and nearest point. I usually put more than 30 points of control.


texts: How to include a quotation mark in a text field from field calculator

texts: I'm trying to attribute pipe size and need to include a quotation mark in a text field.  How do I include that in the field calculator.  "4"" doesn't work.


texts: How to suppress polygon boundaries shared by objects of same type in QGIS?

texts:
Given polygons - like landcover in OpenStreetMap or street areas in cadastre - how can one suppress those polygon boundaries which are shared by objects (features) of same type in QGIS (see figure below)?
AFAIK there exists no such symbolization rules in QGIS 2.6 (yet).
It doesn't seem useful to aggregate (union) the geometries.
But preprocessing could still be an approach, like splitting the common boundaries first, then "suppressing" these when rendering.
Any ideas?


texts: Labelling each single diagram sector with value in QGIS

texts:
I'm working with QGIS 3.10, using the diagrams function for vector layers. Is it possible to label each single diagram sector with an attribute / calculated percentage?
I read this article: Label QGIS diagram sector with percentage etc, but it´s 4 years old.
Is there another way to achieve this today (QGIS version 3.10 or 3.12)? Especially directly in the QGIS.


texts: MapServer's UNION layer - mapfile settings to use attributes from two GeoJSON point files at the same locations?

texts:
First note that we're using MapServer 7.2.2 to display our data from AWS S3 buckets.
We're successfully using MapServer to display point data from GeoJSON files stored in AWS S3.  What we would like to do next is to use the wind speed attribute from one file and the wind direction attribute from another to create a layer that displays wind barb information at those points.  We have used MapServer's uvraster capability for a such a display but these are at grid points versus at selected locations.  So, it appeared that using MapServer's UNION layer feature could work.  The UNION layer would need the speed value for determining which wind barb symbol to use and the direction value for the angle of the symbol.  As a test, just displaying the speed value and attempting to rotate the text with the direction value.
While I can set up the mapfile for the UNION layer to just display values from either the speed or direction layers as a test to see if the data is being accessed from S3 correctly, I can't seem to access both attributes within the UNION layer.  Relative parts of the mapfile:
LAYER

  NAME wspdvalues #wind speed layer
  TYPE POINT
  STATUS OFF

  CONNECTIONTYPE OGR

  CONNECTION "/vsis3/bucketname/windspd_2020-09-11T00:00_2020-09-11T01:00.json"

  DATA "0"

END

LAYER

  NAME wdirvalues # wind direction layer
  TYPE POINT
  STATUS OFF

  CONNECTIONTYPE OGR

 CONNECTION "/vsis3/bucketname/winddir_2020-09-11T00:00_2020-09-11T01:00.json"
  DATA "0"

END

LAYER

NAME pointtest # union layer
TYPE POINT

STATUS DEFAULT
CONNECTIONTYPE UNION
CONNECTION "wspdvalues,wdirvalues"
PROCESSING "ITEMS=spd,dir"

CLASS
  STYLE
   COLOR 0 255 0
   SYMBOL "circle"
   SIZE 3
   MAXSIZE 6
  END
  LABEL
   TEXT (tostring([spd],"%.0f"))
   ANGLE [dir]
   OUTLINECOLOR 0 0 0
   COLOR 255 255 255
   TYPE truetype
   FONT arial-bold
   SIZE 13
   POSITION uc
   BUFFER 4
    PRIORITY [priority]
   END
  END
  PROJECTION
   "init=epsg:920910"
 END
END

Again, this is just a first test rather than using symbols at this time.  If the GeoJSON file has spd and dir attributes in the same file and used in one layer, then it displays fine.  As we're displaying the data values individually, we'd rather not have to duplicate this data by combining into another file.  Hence why UNION looked like it might work.
The geometry sections of each GeoJSON file are identical and in the same order.  Hopefully, just missing something.


texts: Change shapefile encoding to ASCII in QGIS

texts:
I want to convert my shapefile encoding source to ASCII but unfortunately, it is not given in the encoding source type. Maybe it is given with a different name or it is not supported in QGIS?
Unfortunately there is no guide available.


texts: Simplify GeoJSON by retaining external boundaries (clustering)

texts:
I created a website (http://www.cartescolaire.paris) that basically loads a GeoJSON and displays it on a map using Leaflet.

This geoJSON is pretty large (over 2 Mb), the loading time can be very long (it doesn't even load on IE 11). More importantly the resulting map is not very responsive when zooming / navigating.

There are around 110 zones/multipolygons (clicking on a point in the map highlights the zone it belongs to), each of them made from dozens of polygons.

However the only important information that I want to visualize is the external boundaries of each zone. Such a compressed geometry would be much more efficient performance-wise.

The complexity arises from the constraint that the zones shouldn't overlap. The final result should be disjoint clusters.

Any idea how I could do that?

@Stephen : The boundaries are pretty accurate and the orange sections are correctly located (they are actual buildings). I am prepared to do some manual processing if need-be ...

@Nathan : I didn't know turf.js and the buffering technique. I am a newbie in GIS space.
I've already tried Mapshaper, but although it is very efficient in simplifying each single polygon boundary, it doesn't help much in building the external boundaries for a group of polygons. Did I miss something ?

As you properly spotted, several zones have clusters outside the main area and that's one of the main issues.

I will also have a look at leaflet-beta to check if the rendering is better than in v0.7.7

My initial idea was to :

Run DBSCAN algo on each zone to identify clusters of polygons, fine-tuning  epsilon parameter visually
Run a concave Hull on each cluster to identify the boundaries
Finally, simplify the geometry

But I am sure there are GIS tools that can help instead of coding everything myself.

I finally followed Nathan's advice and buffered each geometry (+25 meters then -25 meters) then simplified the geometries up to 20%.
The result is pretty acceptable and as expected the web page loads must faster and the navigation is fluid.

I used QGIS with a bit of Python scripting for buffering and mapshaper to simplify the geometries.

I haven't used turf.js for buffering since I noticed a bug in turf-buffer current version by which some features disappear after buffering.
The issue has already been reported here.