texts: Uploading shapefile to geoserver using cURL

texts:
I'm building the curl command in a python script, but when I run the script I get this error:

* Connected to localhost (127.0.0.1) port 8081
* Server auth using Basic with user 'admin'
> PUT /geoserver/rest/workspaces/csEvo/datastores/ABSM/file.shp HTTP/1.1
> Authorization: Basic YWRtaW46Y3NlbWVyZ2Vuc3lz
> User-Agent: curl/7.15.1 (i386-pc-win32) libcurl/7.15.1 OpenSSL/1.0.1c zlib/1.2.3
> Host: localhost:8081
> Accept: */*
> Content-type: application/zip
> Content-Length: 16302
> Expect: 100-continue
>
< HTTP/1.1 100 Continue
HTTP/1.1 500 Internal Server Error
< Content-Type: text/plain
< Transfer-Encoding: chunked
< Server: Jetty(6.1.8)
Error occured unzipping file:error in opening zip file* Connection #0 to host localhost left intact
* Closing connection #0

but if I run the same built curl command directly in the cmd line it works.

My py script is as follows (zips the shapefile and uploads it into geoserver) and I'm running it from inside the OSGEO4W shell in order to inherit environment settings:

#!/usr/bin/env python
import os
import zipfile
from subprocess import call
# Geoserver Variables
gshost = "http://localhost:8081"
gsuser = "admin"
gspassw = "geoserver"
gsworkspace = "csEvo"
gsstore = "teste"
AB_shpfile = "C:/data/theshpfile.shp"
#zip file for Geoserver upload
thezipfile = zipfile.ZipFile(AB_shpfile.replace('.','').replace('shp', '.zip'), "w")
for filenametozip in [ABSM_shpfile, ABSM_shpfile.replace('.shp', '.dbf'), AB_shpfile.replace('.shp', '.shx'), ABSM_shpfile.replace('.shp', '.prj')]:
       arcname=os.path.basename(filenametozip.replace('.',''))
       arcname=str(arcname).replace(str(arcname)[-3:], '.' + str(arcname)[-3:])
       thezipfile.write(filenametozip, arcname)
# Send to Geoserver
thefileforgs = AB_shpfile.replace('.','').replace('shp', '.zip')
call("curl -u "+gsuser+":"+gspassw+" -v -XPUT -H " + '"Content-type: application/zip"' + " --data-binary @"+'"'+str(thefileforgs)+'"' + ' ' + gshost +"/geoserver/rest/workspaces/"+gsworkspace+"/datastores/"+gsstore+"/file.shp", shell=True)

What am I missing?


texts: What does "TGM" vs "TKG" mean from Sentinel 2 Images (Google Earth Engine)?

texts:
I extracted sentinel 2 images from Google Earth Engine using the Python API with something like this:

# Sentinel 2 Image Collection
s2_image_collection = ee.ImageCollection("COPERNICUS/S2_SR")

filteredCollection = s2_image_collection.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))

When I loop over the images in filteredCollection, I get the following ids:

COPERNICUS/S2_SR/20191002T162111_20191002T162929_T16TGM
COPERNICUS/S2_SR/20191002T162111_20191002T162929_T17TKG
COPERNICUS/S2_SR/20191005T163131_20191005T163713_T16TGM
COPERNICUS/S2_SR/20191005T163131_20191005T163713_T17TKG
...

I was surprised to see two images for each date. What does "T16TGM" vs "T17TKG" mean?


texts: GDAL OGR API GeoJSON from other datasource as string NOT file

texts:
Using the GDAL/OGR API, I need to be able to convert various types of GDAL datasources into GeoJSON.  However, I do not want to write the GeoJSON to a file, but rather I need to just return the text string within the application.

For example, if I have a set of shapefiles, I want to be able to extract the data from one of these shapefiles as a GeoJSON text string.

I want to avoid the slow I/O penalty of writing it to a file and then reading the file and then deleting the file.  Such I/O would likely slow things down too much, as I would be running such tasks very frequently (if I can ever get it to work).

Is there some way to get GDAL data as GeoJSON text without writing it to a file first?

Note that I'm using the C API for GDAL/OGR within an application.  So a C or C++ solution would be preferred.  However, I'm also open to other suggestions which may give me other ideas.

PS .  This is on iOS.


texts:
I have a table planet_osm_roads. The table has one column ref that indicates the code of the road.

Let's say that I have a road coded 'AA1' which is represented by 5 segments (or rows in the database)

osm_id  |   ref | way(geometry)
--------------------
 1         AA1    ...
 2         AA1    ...
 3         AA1    ...
 4         AA1    ...
 5         AA1    ...

I want a new table where the 5 segments are collapsed into 1, based on their common value in ref.

Is this the correct way?

CREATE TEMP TABLE ref_roads AS
SELECT ref,
       ST_Multi(ST_Union(f.way)) as the_geom
     FROM planet_osm_roads As f
GROUP BY ref;


texts:
I noticed that the structure and size of the file was changed

pattern https://s3.amazonaws.com/elevation-tiles-prod/geotiff/{z}/{x}/{y}.tif
example https://s3.amazonaws.com/elevation-tiles-prod/geotiff/1/1/1.tif

Old files described here https://gis.stackexchange.com/a/252936/97778

But now (on 11-12 November 2017) I get 512x512 tiles with no tiling.

These files were previously located at

pattern https://s3.amazonaws.com/elevation-tiles-prod/v2/geotiff/{z}/{x}/{y}.tif
example https://s3.amazonaws.com/elevation-tiles-prod/v2/geotiff/1/1/1.tif

Who can tell me. Is there correct place to get terrain tiles?
Or some notification channel about new changes


texts: Not able to update polyline geometry via update cursor

texts:
I am trying to read polyline from the feature class and then add a point to the geometry and the update it in the same feature class....
I am able to read it and add the new point...but not able to save  the geometry?
with arcpy.da.UpdateCursor(near_features, ["objectid","SHAPE@"]) as cursor:
    for row in cursor:
        lnObjid= row[0]
        array = arcpy.Array()
        geom = row[1]
        for part in geom:
            for pnt in part:
                array.append(pnt)
        for ptid,lineid,ptX,ptY in zip(lst_pt_objIds,lst_Ln_objIds,lst_Pt_X,lst_Pt_Y):
            if lineid == lnObjid:
                new_point = arcpy.Point(ptX,ptY)
                new_array = arcpy.Array([array.getObject(x) for x in range(1,array.count)]+[new_point])
                array = new_array
        if array.count !=0:
            geom = arcpy.Polyline(array)
            cursor.updateRow([row[0],geom])
            array.removeAll()


texts: QGIS custom function error when trying to get color of feature on another layer

texts:
In QGIS 3.16.2 (same on 3.16.0) on Win 10, I created a custom function get_color to get a feature's color by layer name\id and feature ID as proposed in this answer by @Dror Bogin: https://gis.stackexchange.com/a/383753/88814
The expression was successfully created as it appears under custom expressions. However, I'm not able to get the function to run - using this data (project + layers). I try to define the color of the points as the color of the nearest polygon, setting the point's color to data driven override with this expression: get_color('Buildings',array_first(overlay_nearest('Buildings',$id))). This returns an error: Evaluation error: $id. See my screenshot:
Using array_first(overlay_nearest('Buildings',$id)) however works fine, it returns a feature no. So $id seems to be a problem together with the custom function get_color in my case.
What is wrong - is there some setting I have to change? I restarted QGIS as well as Windows.
That the solution in principle works is demonstrated by @Dror Bogin who used the very same data as me (linked above) on his Win/QGIS 3.16.1 successfully, see his screenshot, where I get the error mentioned:


texts: Oracle search against intersect of a region that covers the whole world doesn't find any records

texts:
I am new to Oracle Spatial, I am trying to implement a service that will take coordinates from a map, and search our database to find any geometries that intersect with that map view. I works fine for the most part, but when I have the map zoomed out all the way where we can see the entire world, Oracle search returns no records at all.

Can someone see why the following SQL would return 0 records when I have over 20 records in the geometry table?

select g.id, g.name, SDO_UTIL.to_kmlgeometry(g.geom) as kml from geography g where SDO_RELATE(g.geom, SDO_GEOMETRY(2003,8307,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1), SDO_ORDINATE_ARRAY( -180,-90, 180,-90, 180,90, -180,90 )), 'mask=ANYINTERACT') = 'TRUE';

Several sample records inserted into the database:

INSERT INTO geography VALUES( GEOGRAPHY_ID_SEQ.nextVal, 'd2', SDO_GEOMETRY( 2003, -- 2d poly 8307, NULL, SDO_ELEM_INFO_ARRAY(1,1003,1), -- polygon SDO_ORDINATE_ARRAY(-77.01204711074818,38.88772105792418, -77.01200419540395,38.89196328335449, -77.00831347579947,38.89249647861828, -77.00827056045522,38.88768765356904, -77.01204711074818,38.88772105792418) ) );

INSERT INTO geography VALUES( GEOGRAPHY_ID_SEQ.nextVal, 'd3', SDO_GEOMETRY( 2003, -- 2d poly 8307, NULL, SDO_ELEM_INFO_ARRAY(1,1003,1), -- polygon SDO_ORDINATE_ARRAY(-77.05552889816111,38.86892085934577,-77.05329022237119,38.87065233921285,-77.05467677541158,38.87288138804671,-77.05786820241124,38.87248696633578,-77.05840989692072,38.87006762135151,-77.05552889816111,38.86892085934577) ) );


texts: Google Earth Engine unequal histogram values

texts:
I am creating histograms in Earth Engine of pixel values of cloud masked EVI from different years using Landsat 8. For each of the layers I create, I use an identical geometry to clip the output, but the pixel counts in my histograms vary from each other considerably. Below is a picture of several of these histograms exported to excel and displayed as line graphs to make them easier to see.
An extract of my code is below, the full thing consists of the same code repeated several times but with different years instead of '2018' in each section.
Why is this happening? I understand that because of cloud masks the total pixels in each output will vary slightly, but that would produce nothing on this magnitude as it is only very small areas which are masked. I also made certain that my pyramiding policy was the same throughout. I really have no idea why this is happening!

var forest = BROADLEAF.merge(MIXED).filterBounds(geometry4);
////////////////CLOUD MASK////////////////
function maskL8sr(image) {
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;
  var qa = image.select('pixel_qa');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));

  return image.updateMask(mask).divide(10000)
      .select("B[0-9]*")
      .copyProperties(image, ["system:time_start"]);
}
//////////////2018/////////////////////////
var collection18 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
    .filterDate('2018-06-01', '2018-08-31')
    .filterBounds(geometry4)
    .sort('CLOUD_COVER', false)
    .map(maskL8sr);

var composite18 = collection18.median();

// Display the results.
Map.addLayer(composite18, {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3});
//////////////////////////////////////////////////////////
var evi2018 = composite18.expression(
  '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
    'NIR': composite18.select('B5'),
    'RED': composite18.select('B4'),
    'BLUE': composite18.select('B2')
});

var evi2018real = evi2018.clip(forest);

var histogram5 = ui.Chart.image.histogram({
  image: evi2018real,
  region: geometry4,
  scale: 30,
  maxPixels: 1e13})
    .setSeriesNames(['18']);

print(histogram5);
Map.addLayer(evi2018real, {bands: ['constant'], min: 0, max: 0.75}, "2018");


texts:
So I followed the instruction here: Can QGIS read an ODBC connection? for adding a .mdb into QGIS. There are great instructions about 1/3 of the way down. But I still can't get it to open. I on Windows 7 with QGIS 2.4
I restarted QGIS and my computer. I am still getting the error: "Invalid Data Source: ODBC:co2014ProdSumm (btw this is the .mdb) is not a valid or recognized data source. "

These are not vector layers nor is this a ESRI personal geodatabase. This is an Access database with 2 tables I need out of it.

I also found a few other posts
How to access a .mdb file (Microsoft Access) in Qgis (with windows 7)?
How to import info from mdb file into a project?

and this one might be the answer but I don't know what it's talking about:
Microsoft Access mdb "is not a valid or recognized data source" in QGIS

Is there any way to add a .mdb into QGIS? Also I do not have Access on my computer.

Edit 1:I hope this clarifies my issue.
Each time I get Invalid Data Source: C:...Co 2014 Annual Production Summary-xp.mdb is not a valid or recognized data source.
Here is where I added the access connection

Trying to add as a vector file

Trying to add a database, I also tried local:host 3000 The error in the pic is the same no matter the source, name etc... I've tried changing all fields around.


texts: How to get sublayer infowindow to close in Carto.js

texts:
I have a map with 5 different sublayers and am having an issue with the infowindows between sublayers.  If I click on a sublayer point to open an infowindow and then click on a different sublayer to open another infowindow that first infowindow does not close (I'm assuming it's because it's not part of the same sublayer?).  Clicking on individual points of the same sublayer does have the desired behavior of having the first infowindow close when the second one pops up.  Here is the code I'm using:

var layerSource = {
        user_name: 'myuser',
        type: 'cartodb',
        sublayers: [{
            sql: "SELECT * FROM cartodataset WHERE site_type = 'LUST CLEANUP SITE'",
            cartocss: $("#site-styles").text()
          },
          {
              sql: "SELECT * FROM cartodataset WHERE site_type = 'CLEANUP PROGRAM SITE'",
              cartocss: $("#site-styles").text()
            },
            {
              sql: "SELECT * FROM cartodataset WHERE site_type = 'MILITARY UST SITE'",
              cartocss: $("#site-styles").text()
            },
            {
              sql: "SELECT * FROM cartodataset where site_type = 'MILITARY CLEANUP SITE'",
              cartocss: $("#site-styles").text()
            },
            {
              sql: "SELECT * FROM cartodataset where site_type = 'MILITARY PRIVATIZED SITE'",
              cartocss: $("#site-styles").text()
            }]
        }

 var sublayers = [];

 cartodb.createLayer(map, layerSource)
       .addTo(map)
       .done(function(layer) {
           for (var i = 0; i < layer.getSubLayerCount(); i++) {
              sublayers[i] = layer.getSubLayer(i);
              cartodb.vis.Vis.addInfowindow(map, layer.getSubLayer(i), ['site_name','address','site_type','status']);
           }
       })`


texts: I need to display ImageLayer in OpenLayers. Which image formats are supported in OpenLayers?


texts: Can ArcGIS Server render Mapbox Vector Tiles?

texts:
I am planning to build an application using Mapbox GL JS. All the data sources I need to import are stored on an ArcGIS server 10.5.

There are several feature layers on this server which I'll be able to render in Mapbox GL by requesting GeoJSON from the server using code like from this example.

map.on('load', function () {
    // Add a layer showing the city parks
    map.addLayer({
        'id': 'parks-layer',
        'type': 'fill',
        'source': {
            'type': 'geojson',
            'data': 'https://services.arcgis.com/V6ZHFr6zdgNZuVG0/ArcGIS/rest/services/City_of_Redlands_Parks/FeatureServer/0/query?where=1%3D1&outSR=4326&f=pgeojson'
        },
        'paint': {
            'fill-color': 'rgba(200, 100, 240, 0.3)',
            'fill-outline-color': 'rgba(200, 100, 240, 1)'
        }
    });
});

This will work fine for datasets that are smaller and won't slow down the application. There are also larger datasets on the server such as land-ownership which are quite large. Is there any way to request mapbox vector tiles from an ArcGIS server for rendering in Mapbox GL JS?


texts: QSpatialite broken after update to QGIS 2.18

texts:
Lately I updated from QGIS 2.14 (LTS) to 2.18 on Mac OS X. Since then, the QSpatiaLite plugin has ceased to work, yielding following error message:

I have tried the following to repair the problem, without success:

Removing the SpatiaLite DB from the DB Manager Plugin and re-adding it. I am then normally prompted by SpatiaLite to create a new DB, however SpatiaLite can't connect to it, I get the following two errors:

Reinstalling QGIS 2.14 with all package dependencies from http://www.kyngchaos.com/software/qgis
Installing QGIS from brew as

brew install osgeo/osgeo4mac/qgis-lts
Removing all QGIS files with terminal, like:

sudo rm -R /Library/Frameworks/SQLite3.framework/
sudo rm -R /Library/Frameworks/GEOS.framework/
sudo rm -R /Library/Frameworks/GDAL.framework/
sudo rm -R /Users/[myusername]/.qgis2
sudo rm -R /Applications/QGIS.app/


texts: Importing QGIS3 processing tools in standalone script

texts:
Inspired by this question, but for QGIS 3 (and therefore Python 3), I would like to use the QGIS processing algorithms in a standalone Python script without a GUI, with a standard installation of Ubuntu Linux.
A simple import processing is insufficient:
import processing
processing.run(
   "qgis:rasterlayerstatistics",
   {
       'INPUT': .'/somefile.tif',
       'BAND': 1,
       'OUTPUT_HTML_FILE': './result.html'
   }
)

Results in:
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'processing'

How do you get access to the QGIS processing module to use the QGIS processing algorithms in a standalone Python script (not the QGIS interpreter)?


texts: How to Dynamically Calculate Fields in ArcGIS Report?

texts:
I came across How to calculate custom fields in ArcGIS Reports? and it's very similar to what I'm trying to accomplish.  I have two fields that are sums and would like to divide them in order to get a percentage.

Based on the response, I'm assuming a C# expression is needed, but I haven't a clue how to write a C# expression.  I was hoping someone could point me in the right direction, or provide a general C# expression example using sum and/or divide...?

Based on the screenshot below, my expression would be something like: sum(Field1) / sum(Field2)
 --- green / red


texts: Export the arc-node topology table in QGIS

texts:
How to export the arc-node topology table in qgis?

I will illustrate how that would be the idea.

The figure below shows a polygon composed of arcs and nodes. Each color represents an arc. Each arc has a number, just like the nodes. To the right of the figure I present a table that has the arcs and also all the nodes that make up the arcs. Note that the numbers of the nodes are unique. That is, no node has the same number as the other. This is very important to me.

The tool Convert polygon / line vertices to points located in SAGA, in fact, stores the arcs and nodes. However the numbers of the nodes are not unique. And I need the numbers of nodes to be unique. Any suggestions / ideas?


texts: Select by expression doesn't work - QGIS

texts:
I want to make a "selection by expression", as shown on the screenshot.

I choose the FIELD and the VALUE from the side panel, but it's selecting 0 rows (there are actually more than 100...).

I just don't see why.


texts: Using Web AppBuilder Developer Edition with HTTP/HTTPS?

texts:
I am using;

ArcGIS for Server Standard Enterprise 10.4.1
Portal for ArcGIS 10.4.1

Also i am using Web AppBuilder Developer Edition 2.3 and deploying webapps in our server.

Normally when i deploy a webapp to server root directory and reach site a http web adress my app waits at startup loading screen then redirect to https format and starts.

When i create an app and deploy to my server with WAB 2.7 version, with http web adress my app do not start and stuck at startup screen with this error;

 Failed to load resource: the server responded with a status of 404 (Not Found)
main.js?wab_dv=2.7:199 jimu.js init...
init.js:379 Uncaught Error: Invalid template:

    at Function.e.getCachedTemplate (init.js:379)
    at Object.buildRendering (init.js:378)
    at Object.create (init.js:355)
    at Object.postscript (init.js:354)
    at new <anonymous> (init.js:200)
    at Object.loadConfig (ConfigManager.js?wab_dv=2.7:130)
    at Object.initApp (main.js?wab_dv=2.7:229)
    at init.js:211
    at ja (init.js:28)
    at init.js:29
e.getCachedTemplate @ init.js:379
buildRendering @ init.js:378
create @ init.js:355
postscript @ init.js:354
(anonymous) @ init.js:200
loadConfig @ ConfigManager.js?wab_dv=2.7:130
initApp @ main.js?wab_dv=2.7:229
(anonymous) @ init.js:211
ja @ init.js:28
(anonymous) @ init.js:29
ka @ init.js:28
ha @ init.js:28
g @ init.js:31
(anonymous) @ init.js:35
bar7cbf8c34e03e2da1a4f946e2ce708e5c.js Failed to load resource: net::ERR_CONNECTION_TIMED_OUT

But when i add https in front of adress web aplication starts with no problem.

what changed in 2.7 version that i get this error?


texts: how to make interactive map on sphere in unity3d

texts:
I am working on a Unity3d Project and want to show world map using google on a sphere Gameobject or Plane but actually the problem is Google provide static API which don't provide image greater then 640 x 640.
How can I get High defination interactive image or else solution. I heard about tiles but I didn't find it on Google Map.

How you will make map in Unity3d using any Map system like google or OSM. Or QGIS can provide any help in this regard>


texts: Classification Wizard: Object-Based Supervised Classification - should the segmented image create a segmented classified product?

texts:
My goal is to use the values of my 12-band input mosaic, the segmented image I've already created from the 12-band, and create a classified object based image. I'm using the Classification Wizard in ArcGIS Pro 2.6.0 using Object Based and Supervised classification. The 12-band mosaic was the image that initiated the Wizard. I have my segmented image, my schema, my training samples and testing samples ready for the Wizard.
I thought was going to get a product that was a classified image but in the segmented image (hence the objects). But It appears from the pixelated salt-n-peppery look that it classified my 12-band raster (yes that is the image that I started the Wizard from).
How do I get the segmented image to be the classified image, but using pixel info from the 12-band mosaic since it has all my magical indices?
I've done this in eCognition and it was the segmented image that was classified, alas we don't have eCognition any more.


texts: Is ArcGIS archiving based on ObjectID or GlobalID?

texts:
My question is about the "archive" functionality that is available for the enterprise geodatabase tables in ArcMap.

Does it use the ObjectID to determine that an object is the same as an archived one?

Or does it use the GlobalID?

Or is it possible to use any of them?

Do you have any documentation on the subject? I was hoping to find something in ESRI's Archiving documentation, but I had no luck.

I believe this question is pretty much version and platform independent, but if you think differently, let me know and I can provide more details.


texts: Calculating average values of condition or receive values from a table for calculation in another table

texts:
I have a table with data (in this example the field "TSS"), a field time corresponds to the data for each point. I want to create a field with the mean values for each point all time values, without creating new table. Alternatively I would like to be informed of new table the mean values for each point if I fill the table with the measurements and present. Ηere create a numeric field (Double) named "avg" with the condition to enter a function (like AVERAGEIF of excel).

I am relatively new to GIS.


texts: OpenLayers 3 Map center is not getting set at New Jersey

texts:
Map is not getting set at center location.

var vectorLayer = new ol.layer.Vector({
    source: new ol.source.GeoJSON({
        projection: 'EPSG:3857',
        url: 'njcitypolygonv2.geojson'
    }),
    style: function(feature, resolution) {
        style.getText().setText( feature.get('city') );
        return styles;
    }
});

var map = new ol.Map({
    layers: [
        new ol.layer.Tile({
        //source: new ol.source.MapQuest({layer: 'sat'})
        source: new ol.source.OSM()
        }),
    vectorLayer
    ],
    target: 'map',
    view: new ol.View({
        center: ol.proj.transform([40.0000, 74.5000], 'EPSG:4326', 'EPSG:3857'),
        zoom: 3
    })
});

Can anyone please help?


texts: Building a WebGIS with search functionality for a charity

texts:
I am playing with the idea of building a Basic Web based GIS system for a charity.
The problem:

I work for a charity that deals with people with Cleft Palates.  They currently have a website, but need individuals to contact them to identify which branches of health teams / services / charities they have access to.

I am wanting to develop an open sourced system that allows a user to enter a Zip / Postcode, and for a list of the organizations that cater for there area to be produced.

Because this is a small charity finances are tight, so it has to be free / open source.

Hence i am doing the leg work.

I can digitise all the boundaries in question into .shp files.
Also i can probably develop the front end.

Thanks in Advance


texts: How to assign a categorical value to a polygon based on another vector layer

texts: I want to assign a biome (from TNC's ecoregion layer http://maps.tnc.org/gis_data.html) to every country in the world.  For countries that contain multiple biomes, I want to assign the biome that represents the greatest proportion of that country. What is the easiest way to do this on QGIS?


texts: Can we remove notes layer on Openstreetmap?

texts:
Have a look at the following link :

A piece of OSM map

You can see many violet points which is added by people, I want to know if I can load only the road maps without these  violet points ( these points and related icons (such as shoes ,bicycle, climbing, ...) representing restaurants,banks,shops, ... and I do not need them. They make my map busy).

And maybe good to say that I dislike the style of humanitarian/transport/cycle map. Standard style is the best for my work.


texts:
I am using this answer to calculate some basic statistics of some points that fall within the bounds of a polygon (a vector grid), such that:
gridfile = 'grid.shp'
pointfile = 'points.shp'

point = gpd.GeoDataFrame.from_file(pointfile)

poly  = gpd.GeoDataFrame.from_file(gridfile)

pointInPolys = sjoin(point, poly, how='left')

grouped = pointInPolys.groupby('index_right')['X','Y','Z'].agg(['mean'])

grouped.columns = ["_".join(x) for x in grouped.columns.ravel()]

The input point data has X, Y and Z columns. However, it is only returning statistics (mean) for X and Y and the stats for the Z column are not being returned:
            X_mean     Y_mean
index_right
1221        -64.781242  32.439396
1902        -64.781206  32.439096
2412        -64.781169  32.438777

The data is definitely available in the prior step by checking:
pointInPolys.keys()

Index(['X', 'Y', 'Z', 'geometry', 'index_right', 'DN'], dtype='object')

Is there a reason why the Z column stats are not being calculated?


texts:
Does anyone know how to get road coordinates data from OpenStreetMaps?

I know its possible I just don't know how.


texts: Currently, I have two attribute tables my vector table and an additional table for vector interaction. I was wondering if I could merge the attribute tables within the software. As I would like to use the 'delete duplicate geometries' as some of the data would be reoccurring and would like to delete them without manually doing that. Is this possible?


texts: PyQt QTableWidget for QGIS python plugin: infinite loop issue

texts:
I'm using PyQt4 and Python 2.7 to develop a plugin for QGIS 2.

I have QTableWidgets used to display/compare existing information in the database with newly entered information pertaining to the same feature. There are two QTableWidgets on a dialog box--when selecting a row in one table I would like to connect to one method, selecting a row in the other table would connect to another method. Pressing okay on the button box should accomplish tasks in each method, pressing cancel should roll back new edits.

I have two problems:

The first is an infinite loop, connecting repeatedly to the functions, and failure to close the dialog box and complete task whether okay, cancel, or close dialog is pushed. What's puzzling is that something different seems to happen each time--my guesses as to what may be going wrong: something to do with mixed up signals and slots, perhaps my function should return something so I can exit it.

The second, is that cancel doesn't destroyEditCommand(). There are times when destroyEditCommand doesn't work when okay is pushed for checkData() method.

In case this is important information, the buttonBox emits an accepted() signal to my plugin module in the accept() slot and rejected() signal to plugin module in the reject() slot. This was set up in PyQt4 designer.

Here are some pertinent parts of my code:

def updateTableWidget(self, fid):
    ...
    self.dlg.ui.maintTable.itemSelectionChanged.connect(self.updateData)
    self.dlg.ui.existTable.itemSelectionChanged.connect(self.checkData)
    …

def updateData(self):
    QMessageBox.information(None, 'Connection Result', "In updateData method")
    result = self.dlg.exec_()

    if result == 1:
        self.myLayer.endEditCommand()
        QMessageBox.information(None, 'Whoopie PIE', "You have arrived at the terminal")
    else:
        self.myLayer.destroyEditCommand()

def checkData(self):
    QMessageBox.information(None, 'Connection Result', "In checkData method")
    result = self.dlg.exec_()

    if result == 1:
        QMessageBox.information(None, 'Info', "Turn around and try again later")
        self.myLayer.destroyEditCommand()
    else:
        self.myLayer.destroyEditCommand()


texts: How can we display an animated output for watershed management in openlayers

texts:
Here i'm using Tomcat(7.0)+GeoServer(2.1.3)+Openlayers(2.11).

How can we display an output like Flow of water for the forest watershed management application.

I have to show the water flow Based on the slope range range of the area, rain fall(qty in mm), duration of the rain, and so on.
For this right now I have only 2 methods to generate this output:
1)1 is to fill the color of layer by percentage wise.
2)2nd is to render the layer by percentage wise.

can Anyone guide me how it is possible in Openlayers?

All help is greatly appreciated...

Best Regards,

RamRamesh..:)


texts: Updating field from table used as source in network dataset

texts:
I have a Feature Dataset that contains a Network Dataset and a Feature Class (as a source of Network Dataset) in enterprise geodatabase (PostgreSQL). The Feature Dataset is versioned because we want to edit/update and build the network dataset.

In that feature class, I have a field called condition that will be updated if a road segment damaged or in good condition. I set a restriction in the network dataset if the value is damage, the route should not passing that road segment. I planned that value will updated by a system that querying right to the Default version through a versioned view.

Until now, I successfully updated the default version. I checked in arcmap and see the changes of the data. But, the network dataset that I build is not changes. The route still passing through the damaged road.

How should I update a versioned data for a network dataset through SQL?


texts: Open a table for editing in a Console outside of ArcMap

texts:
I'm trying to write a C# console application to convert a CSV of Lat,Lon data into a Table using insert/update cursors. However, I am having the darndest time figuring out how to get a reference to an ITable without first opening an instance of ArcMap. How would it be done?

The purpose is to open a table inside a geodatabase that exists. After I have a reference to it, I can create an insert cursor and turn the CSV into a proper table inside the geodatabase directly without bothering with Access (the devil).


texts: Configuring QGIS new OAuth2 support for WMS or WFS

texts:
Has anyone successfully configured QGIS new OAuth2 support for a WMS or WFS service?

I have an ArcGIS 10.6.1 service published as a WFS service but cannot get past the OAuth2 authentication.


texts: aerial photo to georeferenced mosaic with X,Y, Z, omega, phi, kappa

texts:
It's worth mentioning that I'm not very experienced.

I have a bunch of aerial photos acquired with digital camera. Photos were taken by night to detect streetlights. Each has 60% forward overlap. I also have file with photo_id, GPS time, X, Y, Z, omega, phi and kappa.

I would like to georeference each photo (creating any world file out of text file) and than create mosaic for further investigation.

Any idea how to do that in QGIS?


texts: I have a feature class with mostly unpopulated fields (ArcMap Standard 10.1, in an SDE database, with attachments).  I exported the attribute table as a .csv file so that a co-worker without access to ArcMap could do the research.  He filled out the .csv file and now I need to add the data back into the feature class.  I am not having any success with "Load Objects" because it just adds all of the data to the end of the attribute table instead of joining the fields. I am not having any success with "Join" because it duplicates the fields.   Merge and append require the data to be of the same type (table and table or feature class and feature class).  I feel like I need a "link" tool... Does anyone know how to do this?  Thanks!


texts:
I'm trying to offset a symbol X feet due East from it's original location, where X is controlled by a real data column ("DBH").

I've gotten as far as making a new marker layer in my style, and changing the units to "map unit".  I can offset it a fixed distance, but am stuck on how to do it with a field.

make_point($x+"DBH", $y)

seemed logical to me, but it isn't working.  I've tried type casting DBH with to_real, I've tried geom_to_wkt.


texts: Subtracting features from one layer using another layer in ArcGIS Desktop?

texts:
I have 2 layers of land parcels for a particular area (shown on top of each other like so):

The darker parcel lines and road polygons are from an earlier parcel map. The land has been developed and many more small land parcels have been added. I want to create a layer that shows all of the new land parcels (i.e. the small parcels that are lightly highlighted) only.

Is there a certain tool in ArcGIS Desktop that I can use to create a layer that has only the new parcels in it (the lightly shaded boxes)?

I have tried performing different types of clipping and intersection operations but I have not been able to generate what I'm looking for.


texts: Merge (a lot of) features using QGIS or other open source

texts:
I have one shape-file with a lot of features (polygons) and if I try to merge these into one big polygon in QGIS it takes days to complete the merge. So I wonder if there is another way? I saw this question:

How to "Merge Selected Features" with python?

The problem is that I don't know python and I have no idea of what they are talking about... So, could someone either give me some advice to perhaps do this operation from command line (I'm using Mac) or maybe if someone has a ready python script that does the trick?


texts: Group polygon features to match a set of specifications

texts:
I have two different sets of polygon features (398 census tracts and 80 ZIP codes) that each roll up to a larger feature (a US county). Though census tracts are smaller than ZIP codes, they do not roll up (i.e. nest within) ZIP codes.

My question -- is there a method/tool using ArcGIS or QGIS (or any software) to separately group the 398 census tracts and the 80 ZIP codes to form 10 polygon features while minimizing the difference between two resulting sets of 10 polygon features?

To clarify, I want to group the 398 tracts -> 10 features, and then separately group the 80 ZIP codes -> 10 features, so that I have two disparate sets of 10 features each. I want to optimize this grouping so that the overlay between these two sets is maximized (i.e. minimize mismatch).


texts: Coordinates of raster image from .TAB file

texts:
I have got a raster image (tif) from a mapping site and an associated .TAB file.  I am trying to correctly map it into CAD/3d model - which I know is set up to the correct coordinates.

The .TAB file is as follows:

!table
!version 300
!charset WindowsLatin1

Definition Table
  File "Getmapping_Index_Map.tif"
  Type "RASTER"
  (336451.73,971812.515) (0,591) Label "Pt 1",
  (336451.73,974812.515) (0,0) Label "Pt 2",
  (339451.73,974812.515) (591,0) Label "Pt 3",
  (339451.73,971812.515) (591,591) Label "Pt 4"
  CoordSys Earth Projection 8, 79, "m", -2, 49, 0.9996012717, 400000, -100000
  Units "m"
  RasterStyle 4 1
  RasterStyle 7 16777215

I assumed that the coordinates of the corners of the image were (336451.73,974812.515) and (339451.73,974812.515) but while that gets it centred on the right point it is not the right scale.  I can't work out what the (0,0)(591,591) coordinates are for.  It seems to be coming out about 3 times the correct size.

(I can't see a similar question, although that may be due to not knowing the correct terminology.  If so my apologies.)


texts: Finding coinciding points in point shapefile using ArcMap?

texts:
I have a point SHP with some points overlapping, they have the exact same X & Y coordinates, I want to extract these points to a new table.

I tried FindSimilar python script but no results.


texts:
The Google Maps in my OpenLayers plugin in QGIS 2.12 does not show anymore. The other maps do show, that is:

OSM
Bing
MapQuest
Stamen
Apple Maps

I tried in QGIS 2.8 and Google Maps does not show there either. This is the case in both my home and work computers, which are on different networks. Did the Google Maps go away?


texts: what is the corresponding projected coordinate system to WGS84?

texts:
Possible Duplicate:
  Global DEM to Slope calculation

I want to calculate the slope from my global DEM?the DEM has WGS84 projection(degrees).In order to correctly calculate the slope from DEM,it has to be projected coordinate system not geographic coordinate system.it is really difficult to say which one of the  projections under the list :
projected coordinate systems is the corresponding to WGS84(world) under the list geographic coordinate systems.Any ideas?


texts: How can I pin all my labels in QGIS 2.2/2.4 (Data Defined Labels) without clicking every single label?


texts: Assigning multiple styles for same layer from GeoServer in OpenLayers 3?

texts:
This is my code, the function set style is not working osm background layer. I have many sld files for one layer and I want apply them on the same

layer:


texts: Creating Digital Elevation Model from Contour Data using QGIS?

texts:
In ArcGIS I used to create DEMs from contours with ease. However, we recently opted to stop using ArcGIS (for various reasons), and switching to QGIS.

The reason I use contour data is it is all that is available to me from our surveyor general (20m contour lines).

How can I convert them into DEMs?

There is no step by step GUI tut for GRASS on the matter, GDAL only uses point data, and IDW Interpolation does not produce a good quality DEM either.


texts:
Currently I am working on a GIS project where I have to update a field in a layer with some values to get my work done. Because the manual efforts are time consuming hence I added a custom addin to my project, but the addin is not updating the fields. Below is my code on the button class.

class C(object):
"""Implementation for Update_addin.btn5 (Button)"""
def __init__(self):
    C.enabled = False
    self.checked = False
def onClick(self):
    sel = ""
    mxd = arcpy.mapping.MapDocument("CURRENT")
    layers = arcpy.mapping.ListLayers(mxd)
    for layer in layers:
        desc = arcpy.Describe(layer)
        if desc.fidSet:
            sel = layer.name
            break

    if sel != "":
        if arcpy.Describe (sel).FIDSet:
            if selection == "Water Projects":         #for water projects.
                with arcpy.da.UpdateCursor (sel, "ACTION") as curs:
                    for row in curs:
                        row = ("C",)
                        curs.updateRow (row)

            else:          #for other projects.
                with arcpy.da.UpdateCursor (sel, "Remark") as curs:
                    for row in curs:
                        row = ("C",)
                        curs.updateRow (row)
        else:
            print("No selection")
        arcpy.RefreshActiveView()
    else:
        print("No selection")

After updating the "ACTION" field the selected feature should moves to another layer based on the query statement.
I am getting below error in the console window.

Traceback (most recent call last):
File "C:\Users\ks38460\AppData\Local\ESRI\Desktop10.3\AssemblyCache\{4EE0D60C-BD1E-4348-8B27-CAF3B9D13FF6}\Update_addin.py", line 41, in onClick
desc = arcpy.Describe(layer)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 1246, in Describe
return gp.describe(value)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing\_base.py", line 374, in describe
self._gp.Describe(*gp_fixargs(args, True)))
RuntimeError: Object: Describe input value is not valid type

Could someone help me understand and rectify the problem in the above code.

I am using ArcMap 10.3 and Python 2.7.


texts: Data Extract Widget in flex viewer 3.6 does not contain output data

texts:
I am currently trying to configure a map service with flexviewer 3.6 on ArcServer 10.2. I have already published my .mxd file as a map service and it displays correctly in my flexviewer application.
Now I want to use the extract data widget to download specified feature layers. The widget works fine and it downloads the .zip file, but it does not contain any data. Did I configured my map service correctly or what's the underlying problem?

I have already read the related questions and the help section of Esri, but I have not found an answer yet.

This is the relevant part of my config.xml file:

    <operationallayers>
        <layer type="dynamic" label="Thermochronology Kyrgyz Tien Shan" imageformat="png8" url="http://MYSERVERIP/arcgis/rest/services/Kyrgyzstan_Geochronology_TEST/Thermochronology_Kyrgyz_Tien_Shan/MapServer" visible="true"/>
    </operationallayers>
</map>
<widgetcontainer layout="float">
    <widget url="widgets/DataExtract/DataExtractWidget.swf" config="widgets/DataExtract/DataExtractWidget_1.xml" label="Datenextraktion" icon="assets/images/i_scissors.png"/>
</widgetcontainer>


texts: Merging two "text string" columns in QGIS?

texts:
I have edited a cvs/xls (2003 format) file to make sure certain columns are split in two, so I could make two drop down list filtered to each other. At least that was the plan, I have not done this yet. Now it is just some text string, so all the points of my shapefile what the table could have been joined to are filled in. But considering I am gonna create some database trigger, I want it back to the way it was before.

This is what it was originally:

And now it is like this

Column1-Colunm2

Begeleidingsnoei-Extra

Now I want to have it back like the original. Written out the same like that.

How do I do this?

Edit: I tried to use the case expression but no matter what this happens:

What I am doing wrong here?

Edit 2: Even using Case Null has the same effect.


texts: QGIS Dissolve/merge features that overlap

texts:
I have a vector layer that I created through classifying an image, isolating a specific class (raster calculator), polygonize, and then created a buffer. I am trying to merge/dissolve/union the polygons that are either touching or overlapping. I've tried dissolving and technically it works, but I am only left with one polygon and I need each one for a total count.
Any Suggestions?
Input

Input Attribute Table

Desired Output

Output Attribute Table. Need each polygon shown and not merged as one. So it should look like image 2, but will less features.


texts: What is the ArcGIS equivalent to WEED tolerance?

texts:
Example:
User is digitising new (line, poly) features and accidentally adds additional vertices within very small distance of previous vertex, resulting in redundant vertices. I want to prevent this from happening within the current editing environment (ArcGIS 10.x). Back in workstation this could be accomplished using WEED tolerance parameter in the ARCEDIT editing environment. It seems there is no equivalent in ArcGIS.

I've been offered a couple of answers to this question, both of which have been inadequate thus far:
Answer - Change the XY resolution or XY tolerance of the feature class.
I absolutely do not want to do either of these as it creates geometry problems in the data. I don't want to change anything about the storage of the data itself, this would be like changing the coverage precision.

Answer - Use the 'classic' snapping environment and snap to sketch vertices and set the tolerance to a minimum distance value.
While this might eliminate some of these occurences, it does not actually prevent the user from creating redundant vertices. It also requires additional validation work afterwards to identify and fix possible errors (check geometry, topology etc.). I am trying to minimise errors at the editing stage while also reducing the amount of work to validate the data.

So the question boils down to "I used to be able to do this? Why can't I now?". The answer could well be "There is no equivalent, deal with it". I'm hoping this isn't the case though.


texts: I discovered that Open Street Map in QGIS would not print properly because the street names were illegible. I looked through some questions that others had posted and saw that I needed to install the Tile Map Scale plug-in, but afterwards I still wasn't able to print and see the names of the streets at any scale. Can anyone tell me how to properly use this plug-in so that I can read the names of streets in OSM?


texts: raster2pgsql command line issue (Windows 10) - ERROR: syntax error at or near "raster2pgsql"

texts:
I'm using the OsGeo4W64 shell that comes with QGIS (I'm using Windows 10) and have managed to log into my postgres database. However I'm looking to import a raster into a new table, and all roads seem to point to raster2pgsql.
I've installed postGIS through the CREATE EXTENSION postgis; and CREATE EXTENSION postgis_raster; commands in my database and I have installed postGIS 3.1.1 through the 'Stack Builder' tool without issue.
Now when I connect to a database and type in raster2pgsql;, or if I type a full command like raster2pgsql -s 4326 -t 1x1 -f DN -I -Y D:/one.tif rasterSchema.RasterName > raster.sql ; , I get a syntax ERROR:  syntax error at or near "raster2pgsql" . I assume it's some kind of pathing issue?
I actually get the same result with shp2pgsql too...
I'm quite new to postgres and don't really know what to do about this


texts:
I work with QGIS and I tried to convert polygon into point. I need points all over the perimeter of the polygon. I saw several answers like:

Problems converting polygon to point

but I didn't see an answer.


texts: Using Python nested if statement for labels in ArcMap?

texts:
I am trying to use python to write a label expression that will display a specific set of labels based on nested if statments. The expression should result in a 4 line label if the value is greater than zero. It should be a 3 line label if the value is equal to zero, basically removing the zero from the label.

def FindLabel( [OPERATOR], [WELLNAME], [WELLNUM], [EURMBOE_LM], [PERMIT_DAT], [SPUD_DATE_], [COMP_DATE_], [WELL_STATU]):
    lbl = [OPERATOR] +'\n' + [WELLNAME] +'_'+ [WELLNUM]
    if float([EURMBOE_LM] >0):
        if [WELL_STATU] == "Permit":
            return lbl + '\n' + [PERMIT_DAT] + '\n' + [EURMBOE_LM]
        elif [WELL_STATU] =="DUC":
            return lbl + '\n' + [SPUD_DATE_] + '\n' + [EURMBOE_LM]
        elif [WELL_STATU] =="PDP":
            return lbl + '\n' + [COMP_DATE_] + '\n' + [EURMBOE_LM]
    else:
        if [WELL_STATU] =="Permit":
            return lbl + '\n' + [PERMIT_DAT]
        elif [WELL_STATU] =="DUC":
            return lbl + '\n' + [SPUD_DATE_]
        elif [WELL_STATU] =="PDP":
            return lbl + '\n' + [COMP_DATE_]

Is the expression wrong?

See pic of labels still displaying the 0


texts: Draw two separate lines from A to B and B to A

texts:
I have an Origin-Destination table like this

Origin Destination Value
A B 10
B A 100
A C 120
...

For all the locations A, B, C... I have another table indicating their location (Lat and Lon)

Is there any way to draw two separate lines between A and B so that I can visualize the directionality and distinguish the flow volume between these two places?

Thanks!


texts:
I am using ArcGIS 10.2 to add raster layers to a geodatabase. When I create a new 32-bit floating-point raster layer, I specify that it should use "Packbits".

After the layer is created, I check Layer -> Properties and find Compression = NONE. I expect it should say Packbits.

The same is true for LZW (Compression = None)

Meanwhile, RLE shows Compression = LZ77

Why am I getting these unexpected behaviours?

No errors/warnings appeared during layer creation.


texts: Producing GIS data from large amount of ASCII DEM files?

texts:
I have just received a large amount of ASCII DEM files. In total 564 Files each containing 40,000 point data (22.560,000 Points). I would like to convert this into a point cloud with a view to interpolating to raster for analysis and cartographic  output. I have both ArcGIS Advanced and FME 2014 to work with. I did this before a few years ago and had problems importing into Excel due to restrictions on the amount of rows excel can read/write. Does anyone do this kind of work regularly and could advise me on suitable workflows to produce raster data from these files including the prep-work to create point cloud data, and also the interpolation process to create realistic and accurate raster data.

The data format is one point per line...unfortunately each file has the file name written on line 1, and there are other data written in the columns such as a code for the data source as well as slope (gon) and aspect (gon).


texts:
Is it possible to create a definition query that will return a subset of values based on a specific field length?

I only want string values where length is less than or equal to 2 characters.


texts: Configuring GRASS 7 to work alongside QGIS 2.18?

texts:
All day I was struggling with QGIS 18 and GRASS 7.
My mission was simple: import a file in QGIS 2.18, run v.clean and then go home. Unfortunately it turned out that my mission was not that easy.

First of all, although I read everywhere that GRASS is now a component of QGIS 2.18, when I tried running the plugins, received an error telling me that I needed to configure the working directory for GRASS.

I tried different things (see the list below) doing that and I failed many times.

Installing GRASS 7 separately and moving the contents of the folder in the apps subfolder of QGIS, in order to be able to point it in the program
Found out that there already was a grass installation in my QGIS folder. Linked the application to that folder
Read what was written in the QGIS documentation about configuring grass, followed the steps, failed.
Installed MSYS, linked the folder to the corresponding one in QGIS (heard that GRASS 7 did not require MSYS)
Started the GUI of the software, got an error while trying to import the shapefile, a module (forgot the name) was missing and the command failed.

Can anyone guide me step-by-step on how I can configure the software in order to work properly?

Thinking as a last resort downgrading to QGIS 2.14 and hoping for the best.


texts: Generating Points Along Line at irregular intervals using ArcGIS Desktop?

texts:
What tool do I use to generate points along a line from an input list of lengths?

For example:

Length_km
0.12
0.2
0.38

etc.

ArcPy solutions are okay.

I tried:

testOut = outpath + '\\' + 'testOut'
arcpy.LocateFeaturesAlongRoutes_lr('pts.shp', 'polyline.shp', route_id_field = 'OBJECT_ID', radius_or_tolerance = "", out_table = testOut, out_event_properties = "", route_locations = "ALL", distance_field = "DISTANCE", zero_length_events = "ZERO", in_fields = "FIELDS", m_direction_offsetting = "M_DIRECTON")

& received the following error:

ERROR 000709: The type of out_event_properties does not match the type of in_features.
Failed to execute (LocateFeaturesAlongRoutes).

EDIT

Based on @jbalk 's suggestions, I tried the following code (note that my .py script and .shp, .dbf files are in the same directory, so just calling .shp without a path works here):

makeRoute = outpath + '\\' + 'make2route.shp'
arcpy.CreateRoutes_lr('polyline.shp', "routeID", makeRoute, "LENGTH")
testOut = outpath + '\\' + 'testOut.shp'
arcpy.MakeRouteEventLayer_lr(makeRoute, "routeID", 'pts.dbf', "routeID", testOut)

And it gives the following error:

File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\lr.py", line 528, in MakeRouteEventLayer
    raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000735: Event Table Properties: Value is required
Failed to execute (MakeRouteEventLayer).

Below are screenshots of the attribute tables for:

input 'polyline.shp'
output 'make2route.shp'
input 'pts.shp'

Just to reiterate my objective:

I have "polyline.shp"
I have "pts.shp" that are not located on the
line. However, there is an attribute field that gives the distance
along the line ("Dist_km"). I don't care to recalculate/collapse to the line; I must use the "Dist_km" field for this analysis.
I need coordinates for these distances ("Dist_km")
along this line ("polyline.shp")

#


texts: Temperature and Moisture indices in Geomorphometry and Gradient metrices toolbox

texts:
I have generated various temperature and moisture indices using toolbox in arcgis. But HLI output is showing values more than 1. values should be between 0 and 1. similarly IMI values should be between 0 and 100, but i am getting much higher values. What may be the reason?
Is there any need to do focal mean of these outputs using 3*3 windows?


texts:
I have several arrays obtained from K means, each array contains a variable number of classes, e.g. from -1 to 3, from 0 to 4, from -1 to n. I would like to be able to save each class of an array as a mask in a new array with the same shape as the original one. Bearing in mind that the number of classes for each array is variable and saving only classes other than -1.
with rasterio.open(cluster[0])as ds:
    arr=ds.read()
arr
array([[[-1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1.,
         -1., -1., -1., -1., -1.,  1.,  1., -1., -1.],
        [-1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1., -1.,
          1.,  1.,  1.,  1.,  1.,  2.,  2.,  1., -1.],
        [-1., -1., -1., -1., -1., -1., -1., -1.,  1.,  1.,  1.,  1.,
          2.,  2.,  2.,  2.,  0.,  2.,  2.,  1., -1.],
        [-1., -1., -1., -1., -1., -1.,  1.,  1.,  2.,  2.,  0.,  0.,
          0.,  0.,  0.,  0.,  2.,  2.,  2., -1., -1.]]], dtype=float32)

class = np.unique(arr)
array([-1.,  0.,  1.,  2.], dtype=float32)

#Saving classes on an array works, but I don't know how to do that for many arrays
cls1 = (arr==0)
cls2 = (arr==1)
cls3 = (arr==2)

#I tried something like that but it doesn't work.

for s in class:
    for k in arr:
        cls = s == arr
        print(cls)


texts:
I'm new to GIS, and I have some tasks that have to work with lat and lng values. The problem I have is this. I received some Excel files from a data source that supposed to contains some buildings' lat and lng in area A, and what I got was just have x and y (I don't know what does it mean). I know there are several different values like map xy, so I think those are one of the Geographic coordinate system that I don't know. I have to figure it out what is it and have to convert it into lat and lng that you can use it in Google Maps. I wonder if it is possible?

So, the numbers very strange, the area A supposed to have like lat 35.XXXXX and lng 127.XXXXX and what I got was something like x: 12834.6278 and y:3551.5214. At first I thought that those are map x and map y, because one time data source gave me that kinds of numbers and I had to find tool to convert it into lat and lng. However when I use that tool and put numbers I got into Google Maps, it pointed area where not even close to area A.

those are sample data

     x     |      y
12834.6278 |  3551.5214
12834.6689 |  3551.5184
12834.8509 |  3551.4267
12834.8004 |  3551.4605

and this is lat and lng inside of area A, the format of expected numbers.

     x     |      y
35.856770  |  128.576365
35.860748  |  128.570098


texts:
Wow I've just discovered the world of TRANSIMS, a transport model with ambitious aims: "to model all aspects of human behaviour related to
transport in one consistent simulation framework" (Nagel, 1999). I was sceptical at this claim at first, but further exploration has revealed that it is increasingly being adopted by Municipal Planning Organisations in the US, is completely open source (complete with a NASA license), and works on Linux as well as Windows.

All this is very exciting for someone researching modelling travel flows.

My question to the GIS community is has anyone got this to integrate with other OSGEO tools such as QGIS, PostGIS or PGrouter (although TRANSIMS has its own route planner)?

There is also a python-based visualiser for this software:
http://sourceforge.net/projects/transimsstudio/?source=recommended

Hope this set of links is of use to others, and hope that people who've encountered TRANSIM can provide insight into its integration with GIS.

Nagel, K., Beckman, R. J., & Barrett, C. L. (1999). TRANSIMS for urban planning. 6th International Conference on Computers in Urban Planning and Urban Management, Venice, Italy. (Available on arXiv).


texts: Joining .csv with shapefile in ArcGIS Desktop?

texts:
I'm trying to join a .csv file with a shapefile in ArcGIS. In both the files theres a column named BlockID that are going to be joined the joined attribute between the two files. My problem is that when I do the join validaton it gets stuck at "counting the number of matchingrecords for the join" and after a while it says that the join is cancelled. I have done joins before, big ones like the one I'm doing this time, and never had this problem.

Does anyone know what might be the reason why the validation just cancels?


texts: Is there a QGIS plugin that will allow me to embed an interactive web map on a website that requires limited coding? If so, what plugin(s) would be useful for this? I do not need to use open source layers. I already have all of the data.


texts: Using Boolean And to overlay layers in ArcGIS Pro?

texts:
I am working for a project related to a site selection in which I have 5 layers to overlay to have as result:

Zone1: area suitable
Zone2: area not suitable.

Boolean And seems to be the perfect function to use, but it only allows me to combine 2 layers.

Is there any other better function or should I proceed is steps (2 layers by 2)?

I am using ArcGIS Pro with an Advanced level license..


texts:
I am currently trying to learn how to use Leaflet / Mapbox (and JS, to some extent) by creating custom test maps.

My problem is that I want to add a logo on the top of my map, with absolute coordinates (in pixels).

I have already looked for out-of-the-box methods on the API reference and GIS stack, but couldn't find anything. I also tried some "homemade" solutions (with my limited skills), like adding the logo as a marker or custom layer with non-fixed geographic coordinates (using getBounds().getEast() for example). Unfortunately, nothing worked...

It can be easily done with the Google Maps API, using :

var controlUI = document.createElement('div')

and the logo passed as a background for the 'div' section, but I would like to avoid this API.

Here is an example of what I want it to look like :


texts: Google Maps Distance Matrix API - Queries per second?

texts:
I'm pretty familiar with Google's Distance Matrix API, but I'm unsure about one aspect of the Usage Limitations. The documentation states that you have access to 100 elements per query (element = origin * destination). So one could calculate up the distance of 10 origins and 10 destinations in one request, I understand this.

Where I'm stuck is where Google mentions that developers can request 100 elements per 10 seconds. What exactly does this mean? Is it saying that if I were to make a request of 10 origins & 10 destinations (100 elements), it would take 10 seconds to compute this?


texts: Choosing projection for running polygon intersections at global scale (i.e, geodesic intersections)?

texts:
How do I choose a projection for intersecting two global polygon datasets in WGS-84 (4326), and get accurate estimates of their overlap area?
I searched the internet for answers to this question, and all turned out to be wrong I think, except for the lovely projection table posted by Jochen Albrecht.
So, when you are intersecting two polygons in 4326, you will need currently in R (sf package, as of 12/2020) to reproject them into a common projection, run the intersection, and then back-transform them into 4326 to get a correct estimate of the overlap area (using sf::st_area).
Most sources (e.g. Chapter 6 Reprojecting geographic data of Geocomputation with R) will recommend an azimuthal equidistant or equal-area projection centered on the focal polygon, which if you have lots of polygons (I have >200K), means a for-loop.  Not that slow actually, but many of the suggested projections will do geodesic intersections between polygons with some bias.


texts: I tried this example Processing vector to raster faster with R on Windows, where the option "FORK" for makeCluster is not available. Defaults to "SOCK" or can be set "PSOCK". When running into the decisive line: system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2'))) I get the error: Error in splitIndices(length(x), ncl) : argument "x" is missing, with no default Timing stopped at: 0 0 0 Is the option "FORK" absolutely necessary to run this code? Or should it work with options "SOCK" and "PSOCK" as well? Any hint welcome!


texts: Geoserver "Layer Groups" with Raster Layer

texts:
I need to group multiple geotiff raster layers together. Each layer works well on its own when being previewed through the "Layer preview" utility.

However, when I group the layers together through the "Layer Groups" utility, the created grouped layer is previewed as an empty/white layer.

I am a bit suspicious that the problem lies in the coverage band details (palette_index). I am a web developer who is trying to integrate the layers in a website and therefore I don't have an experience in what does the coverage band details exactly mean.

I am attaching the configuration screen of two of the tiff files I am referring to.

Update 1:
I have now pasted the geoserver log, does this show a memory related issue? One line says: 2015-10-25 06:02:42,864 ERROR [geotools.map] - Call MapContent dispose() to prevent memory leaks. You can view the log of 1000 lines here: http://jpst.it/CFPL


texts:
I've come across this problem a few times and have not been able to find a satisfactory answer. I am new to Python though so perhaps I have not been looking in the right places!

Anyway, In ArcGIS I have been using the following code to select entries by attribute in the attribute table:

import arcpy

arcpy.env.workspace = "C:\Users\T0015685\Documents\WCGIS\Geog485\Lesson1\us_hydro.shp"

qry = '"COUNTRY" = \'CANADA\''
flayer = arcpy.MakeFeatureLayer_management("us_hydro", "us_hydro_Hudson6")
arcpy.SelectLayerByAttribute_management(flayer, "NEW_SELECTION", qry)
cnt = arcpy.GetCount_management(flayer)
print "The number of selected records is: " + str(cnt)

This script works fine. HOWEVER, when I alter the query entries from verbal to numeric attributes I get an error. So, for example, if I change the query to - qry = '"HYDRO_= \'1744\'' - the script fails. I think this is something to do with the SQL expression but I'm not sure.

I want to be able to use code such as this to search feature ID's but can't because every time I use numbers the code doesn't work. Could someone please explain why this is the case or perhaps offer a solution.


texts: How to find and delete identical points?

texts:
I recently performed Integrate on a merged point feature class in ArcGIS 10.1.  I specified "1 meter" XY tolerance for the Integrate command.  The resulting point feature appears to be what I want except points that were integrated are now stacked on top of each other, which will be problematic in future processing steps.  The attributes for my dataset consist of an OID, Shape, x, y, z.  Where x and y = UTM coordinates and z = buffer distance.  The command Delete Identical (Data Management) does not appear to work for my dataset (i.e based on Fields: OID, Shape, x, y, z and XY tolerance = 1m) as the output still has stacked points.  I am also working with big data--the integrate output has 37 million points.

There is a related thread, titled Remove spatially duplicate features using ModelBuilder, although I would like to open the question up to python solutions too.3

What is the best way to find and delete (spatially) duplicate points?

Is there a better solution to integrating such that the output does not contain overlapping points?


texts: How to keep feature ID within GeoJSON ogr2ogr import to PostGIS

texts:
I'm attempting to import a GeoJSON file into PostGIS using ogr2ogr
The basic command is
ogr2ogr -f "PostgreSQL" PG:"dbname=my_db user=postgres" my_file.geojson -preserve_fid -nln geo_table -append
Using this sample how do I rewrite it so the id becomes a id Colum in the database table?
{"type":"Polygon","id":123,"coordinates":[[[100,0],[101,0],[101,1],[100,1],[100,0]]]}
Produces
ERROR 1: ERROR:  setval: value 0 is out of bounds for sequence
The "real" file contains both Polygon & MultiPolygon


texts: Multi label with mulit placement in OpenLayers

texts:
I have a vector layer in OpenLayers. It has some attributes like left, right, top, bottom, etc.
I want every feature to be multi-label (to have more than one label).
The javascript code is as follows:

style = new OpenLayers.Style();
rule = new OpenLayers.Rule({
    filter: new OpenLayers.Filter.Comparison({...}),
    symbolizer: {
        fillColor: 'blue',
        label: ?????,
    }
});

style.addRules([rule]);

styleMap = new OpenLayers.StyleMap({
    'default': style
});

myLayer.styleMap = styleMap

I want my style to be like this:

How do I do?

  note: I know that I will be able to do this by "SLD", but I have to use style map in OpnenLayers.


texts: Installing external software for QGIS Plugin

texts:
I'm wondering if anybody has experience in installing external software for QGIS plugins.  I am developing one and I feel like I have hit a wall.  I need to use osm2pgrouting, but it is written in C and needs to be compiled on the client's computer in order to be used properly.

I am experimenting with subprocesses like so:

import subprocess
import os

os.chdir("/media/isaac/Data/1sandbox/osm2pgrouting")

params1 = ["cmake", "-H.", "-Bbuild"]
os.chdir("/media/USERNAME/Data/1sandbox/osm2pgrouting/build")
params3 = ["make"]
params4 = ["make", "install"]

subprocess1 = subprocess.Popen(params1, stdout=subprocess.PIPE)
subprocess3 = subprocess.Popen(params3, stdout=subprocess.PIPE)
subprocess4 = subprocess.Popen(params4, stdout=subprocess.PIPE)

But it is failing and I believe that it is because root privileges are required in order to install this software.  I would like to include a "first time setup" function to compile and install it so that the user doesn't have to bother, but I'm wondering if this is even possible.

Is there a smarter way of going about this or am I just stuck instructing the user to install osm2pgrouting themselves manually before using the plugin?


texts:
I have two sets of linestrings which have to perfectly overlap. The set A have to contains the set B (the A is a combination of multiple Bs which are touching each other).

I would like to build a query to easily spot the vertices where it doesn't match.

For now I developped a query based on ST_DumpPoint on each dataset, and I was checking for the point from A that don't match any point from B. It worked well on small datasets, but the one I am currently working on results in 15K points for A and 9K for B, so something like 150M match to compare. It takes forever (and I stop it after 2min of running).

I also tried to use ST_SymDifference, but it just give me everything when I run it directly on the geometries, and it takes forever if I try to compare St_Union(A) and ST_Union(B), so not a solution either.

I am guessing I am not the first one to have such a need, and hoping there is an easy way to do that which I didn't found yet, do you know such a way?


texts: Postgresql (soln would be basic SQL I think): ORDER by distance to a point and return those distances in array pairs w/ the full records

texts:
so I know there have been some questions asked about how to calculate distance using postgis and how to order results by their distance from a respective point.

I'm not very talented at SQL (at all) and so was wondering if I could get some help returning:
    a.) a list of records sorted by their distance to a point
    b.) with each record in a small array containing the record itself and the distance from the record to the point passed as the second geometric var

I've already learned how to order my results from this question: https://stackoverflow.com/questions/5111627/order-by-nearest-postgis-georuby-spatial-adapter

So I'm guessing maybe my solution might look something like:

matches =ReferencePoint.where(custom_sql).order("ST_Distance(reference_points.geom, ST_GeomFromText('POINT #{other_loc.x #{other_loc.y}', 3361))")

Which I assume works...haven't actually tested it....

...that would be perhaps amended like:

#aborted code

OK, I don't want to run a second select command, I want to add the sql to the code in my order command. Which means I probably have to use #select....don't know how to do that simultaneously after an ORDER command. So yeah, I need some help. Thanks

EDIT:
Example data:
- lets say I had 100 boats scattered about around a buoy. I want sql to return those boats in order of their distance to the buoy, closest first. I'm going to use pseudo-sql code to represent the last part of this (text inside of select method):

buoyshape = Buoy.find(2341).projected_geometry_3361 #I document srid codes in my columns
Boat.order{st_distance(:proj_shape_3361, buoyshape)}.select("* + st_distance(proj_shape_3361, ?)", buoy_shape)

OK so I don't know how to write sql so the last #select component doesn't do what its supposed to do. But basically I want the distance between the buoy and boat to be returned as an additional column for each record.


texts: How to compare a 'predicted' layer with an 'actual' layer for the same study extent?

texts:
I have developed a vector layer which predicts which polygons are more or less likely to experience coral bleaching.
I now have another vector layer with the same extent with actual coral bleaching data.
Each of the layers are split into the same distribution of polygons and each polygon has a value of 1 (high), 2 (medium) or 3 (low) for both datasets.

I want to validate my predicted layer with the actual results in a quantitative manner but I'm unsure of the best way to go about this.


texts:
What are the conventions for the talker ID in e.g. GSV messages for SBAS satellites; is it the ID of the positioning system to which the satellites in the message belong?
Is it allowed to mix SVIDs from different positioning systems, say GPS and GLONASS, that belong to an SBAS in a single GSV message and what would the talker ID be then; GNGSV seems to be deprecated but not forbidden yet.


texts:
Because OGC SFS is relatively strict when it comes to mixing Geometries with their MULTI- Counterparts, we came up with a solution to store everything as MULTI... and convert incoming geometries via BEFORE INSERT OR UPDATE Triggers using ST_Multi(NEW.geom).

So our tables are defined as

CREATE TABLE IF NOT EXISTS basins (
 id serial primary key,
 geom geometry(MULTIPOLYGON,31287)
 ...

);

And a BEFORE INSERT OR UPDATE Trigger ensures that POLYGON Features are also ST_MUlti´d

Now I recognized, that this won't work, because PostGIS won't accept any INSERTs with the wrong typmod, i.e. geometry(POLYGON,31287) altough these would be implicitely convertible to each other.

When I cannot solve this on the database side, I have to solve it on the application side. So the only way to solve it would be defining Views which accept all geometries and convert the geometry in the INSTEAD OF Trigger. Because QGIS for example relies heavily on the correct specific typmods to work in a performant manner, I end up using the Table (with correct metadata) for viewing and the views for writing, which I think is very awkward to use.

Are there any better patterns to streamline this a bit?


texts: Is it possible to efficiently get all the polygons inside a polygon from a shapefile without using PostGIS?

texts:
I have a huge shapefile (over 400MB) of an entire city that contains many polygons of the residential area of this city. I can check it on QGIS as the following:

This shapefile contains a huge number of polygons, and I'd like cut a small shapefile from it in a circular or squared area around a point. So, let's say I give a point to a system, I want to return all the polygons around this area inside an offset (just so I can send only this segmentation of the data to the frontend of my application).
I know I could convert this shapefile to GeoJSON and find what are the polygons that exist inside the area that I have created (using turf.js for example), but that would be a brute force solution and it wouldn't be efficient.
On the other hand, I know I could install PostgreSQL and PostGIS and implement solutions like this one.
My doubt is that if I can achieve the goal of finding all polygons inside a polygon efficiently using only the command line, Python or JavaScript without implementing a spatial index for that, working directly with the shapefile. Is it possible?


texts: How can I use Android to add some additional data to geotiff images?

texts:
I would like to create an application which will use a geotiff image and add some additional data on it. All these with the Android Platform.
I've done some research and found GDAL but I don't know if there is something else which exists?

Thanks

Martin


texts: ArcGIS Online trouble with 3D extrusion in ArcScene

texts:
I'm having difficulty viewing my Feature Class in Scene Viewer where my saved symbology is not appearing as 3D Extruded, but rather as 2D.  I have a Feature Class uploaded to my ArcGIS Online account and opened it in Map Viewer.  Here, I was able to apply my desired symbology and save.  When I view this Feature Class in Scene Viewer, by symbology is correct, but I cannot view it as 3D Extruded.  If I change the Symbols option from "Original Symbols" to "Change Symbols", I'm able to view my data in 3D, however I lose my symbology.

How can I keep my symbology and view it as 3D in Scene Viewer?

Here are a couple screen captures to help illustrate my issue.

Correct symbology, not in 3D

3D, but cannot change symbology


texts: How to clip point features by polygon grid?

texts:
I am trying to adapt the "Batch clip in ArcGIS 10" script I found on GIS.StackExchange to my needs. I would like to do this in Python as a learning experiment. My data:

Point data
Standard numbered grid

My goal is to output point_1, point_2, point_3, {...}, point_n

Help me think about the problem. Since I'm not going to be "clipping" polygon data, can I accomplish this with a batch select? I've been doing some reading about batch clips, but have not written anything useful yet.


texts: Add layers to GeoPackage using PyQGIS (QGIS 3)

texts:
How to add layers (vector or raster to an existing GeoPackage)?

While loading an existing layer of a GeoPackage is quite trivial:

path = '/home/matteo/geopackage.gpkg|layername={}'.format('admin')
iface.addVectorLayer(path, 'admin', 'ogr')

I couldn't find information about adding them via PyQGIS.

Is it possible? Will it also work for rasters?


texts: Distance between DB field and Point using PostGIS

texts:
I'm using Postgresql 9.0 with PostGIS 2.0 extension.
I've got the table and geographic field location:

=# SELECT * FROM geography_columns;
 f_table_catalog | f_table_schema |   f_table_name   | f_geography_column | coord_dimension | srid | type
-----------------+----------------+------------------+--------------------+-----------------+------+-------
          _dev   | public         |          objects | location           |               2 | 4326 | Point

=# \d objects
                                       Table "public.objects"
     Column     |            Type             |                           Modifiers
----------------+-----------------------------+---------------------------------------------------------------
 ...
 location       | geography(Point,4326)       |

And have got records with coordinates:

=# SELECT ST_X(location::geometry), ST_Y(location::geometry) FROM objects;
  st_x  | st_y
--------+-------
        |
 -73.01 | 15.01
 -73.02 | 15.02
 -73.03 | 15.03
(4 rows)

I want to grab records within distance in meters from Point(-70 15).

I've calculated this distances to each point (just for know it):

=# SELECT ST_Distance(ST_GeographyFromText('SRID=4326;POINT(-70 15)'), ST_GeographyFromText('SRID=4326;POINT(-70.01 15.01)'));
   st_distance
------------------
 1543.03968305103
(1 row)

=# SELECT ST_Distance(ST_GeographyFromText('SRID=4326;POINT(-70 15)'), ST_GeographyFromText('SRID=4326;POINT(-70.02 15.02)'));
   st_distance
------------------
 3086.04520065436
(1 row)

=# SELECT ST_Distance(ST_GeographyFromText('SRID=4326;POINT(-70 15)'), ST_GeographyFromText('SRID=4326;POINT(-70.03 15.03)'));
  st_distance
----------------
 4629.016529543

Now I see the closest point in 1543 meters and farthest point in 4629 meters.

I've used very similar query to example from here:

  Query and measurement functions use units of meters. So distance
  parameters should be expressed in meters, and return values should be
  expected in meters (or square meters for areas).

-- Show a distance query and note, London is outside the 1000km tolerance
  SELECT name FROM global_points WHERE ST_DWithin(location, ST_GeographyFromText('SRID=4326;POINT(-110 29)'), 1000000);

=# SELECT location FROM objects WHERE ST_DWithin(location, ST_GeographyFromText('SRID=4326;POINT(-70 15)'), 5000);
 location
----------
(0 rows)

and I didn't get locations within 5000 meters. More than that I cannot calculate distances in meters from location to Point:

=# SELECT ST_Distance(location, ST_GeographyFromText('SRID=4326;POINT(-70 15)')) FROM objects;
   st_distance
------------------

 323718.841976209
 324792.390871594
 325869.606072403
(4 rows)

There numbers are a lot bigger than numbers from distance in Point to Point example.

I'm guessing something wrong with my column because I cannot even grab points with same coordinates:

=# SELECT location FROM objects WHERE ST_DWithin(location, ST_GeographyFromText('SRID=4326;POINT(-70.01 15.01)'), 5000);
 location
----------
(0 rows)


texts: Defining a local CRS that is basically UTM31N rotated by 43.22 degrees

texts:
Working with some data in the North Sea and basically I'm a bit fed up of working at a slant of 43.22 degrees. I was thinking it must be possible to define an orthogonal local coordinate system as a .prj file for my area. However I have no idea how to do this...

The desired centre of such a grid would be at 2°0'0"E 55°0'0"N and rotated by -43.22 or +46.78 so that the boxes in the image below go from rotated to orthogonal.

I managed to solve this myself by defining a .prj using a tutorial found at:

https://geonet.esri.com/thread/11068

For anyone to check, I used the bottom left corner as my origin and set the rotation to -43.22/136.78 and ended up with:

PROJECTION["Local"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Scale_Factor",1.0],
PARAMETER["Azimuth",-43.22],
PARAMETER["Longitude_Of_Center",1.42624],
PARAMETER["Latitude_Of_Center",55.102476],UNIT["Meter",1.0]]


texts:
I am having an issue displaying satellite data on the correctly projection using cartopy and imshow. I am unsure of the initial data projection, but I am trying to put it into Plate Carree with a limited area extent using bbox. There are nans in the array and I tried to use np.ma.masked_invalid to mask.

https://scitools.org.uk/cartopy/docs/v0.15/examples/geostationary.html

I tried the following examples, but it did not resolve the issue.
[https://stackoverflow.com/questions/52883594/projecting-goes-16-geostationary-data-into-plate-carree-cartopy]2

What happens is the data is correct, but it is not plotting it into the correct projection with axes. I overlayed states and coastlines and they are all out of whack. This is using SCMI data for GOES-16 from NOAAPort.

Here is my code and image that shows the correct display (using different approach; last image on bottom) and an image showing the incorrect display (includes little box in top left corner and lighter colors) from this code... for comparison. It should be noted that contourf works fine. I tried to use pcolormesh, but it takes forever and using runs into a mem error.

edexServer='awips'
DataAccessLayer.changeEDEXHost(edexServer)
request = DataAccessLayer.newDataRequest()
request.setDatatype("satellite")
sector="ECONUS"
request.setLocationNames(sector)
request.setParameters("CH-13-10.35um")
cycles = DataAccessLayer.getAvailableTimes(request, True)
times = DataAccessLayer.getAvailableTimes(request)
fcstRun = DataAccessLayer.getForecastRun(cycles[-1], times)
response = DataAccessLayer.getGridData(request, [fcstRun[0]])
grid = response[0]
data = grid.getRawData()
imdata = np.ma.masked_invalid(np.atleast_2d(data))
#where_are_NaNs = np.isnan(data)
#data[where_are_NaNs] = 0
#data = ndimage.gaussian_filter(data, sigma=0.50, order=0, mode='nearest')
lons, lats = grid.getLatLonCoords()

bbox = [-130.0, -65.0, 20.0, 55.0]
fig = plt.figure(figsize=(6,6),dpi=200)
ax = plt.axes(projection=ccrs.PlateCarree())
#ax.set_extent(bbox, ccrs.PlateCarree())
coastlines = '/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/ne_10m_coastline.shp'
states = '/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/ne_10m_admin_1_states_provinces.shp'
globe = ccrs.Globe(semimajor_axis=6378137.0, semiminor_axis=6356752.31414)
crs = ccrs.Geostationary(central_longitude=-75,
                         satellite_height=35786000, globe=globe)
cpt=loadCPT('/home/awips/Downloads/IR4AVHRR6.cpt')
cpt_convert = LinearSegmentedColormap('cpt', cpt)
ir_norm, ir_cmap = colortables.get_with_range('ir_bd', 150, 357)
#cs = ax.pcolormesh(imdata,cmap='coolwarm', transform=ccrs.PlateCarree())
#cs = ax.contourf(lons,lats,data,cmap='Greys_r', transform=ccrs.PlateCarree())
#levels = np.arange(170, 357, 5)
#image = ax.imshow(data, origin='upper', extent=(-149,-53,14,51), cmap=cpt_convert, vmin=150, vmax=375)
image = ax.imshow(imdata, origin='upper',cmap=cpt_convert, vmin=150, vmax= 375, transform = ccrs.PlateCarree())


texts: Improving ArcPy script for ArcMap to only return rows of selected polygons instead of all polygons in layer

texts:
I am debugging a code provided to me that merges two polygons and alters the fields of "FID1' and "CALC_AREA."

I have realized that the current code searches every layer in the map document for fields containing "FID1" & "CALC_AREA", sorts them ascendingly, chooses the ID matching the largest "CALC_AREA" selected, and assigns the value to all IDs of merged polygons.

def populateMergeID():

mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd)[0]
for lyr in arcpy.mapping.ListLayers(mxd):

    # Describe cursor object for selected polygons sorted by area
    search = arcpy.SearchCursor(lyr, fields="FID1; CALC_AREA", sort_fields="CALC_AREA A")
    FID1 = ""
    area = 0

    for row in search:
        print "Current FID1!: "
        print row.FID1
    # Just want the first row FID1 value
        if row.CALC_AREA > area:
            area = row.CALC_AREA
            FID1 = row.FID1
# Corresponds to largest polygon in selection
    print "The container polygon has FID1 value!!: " + FID1

    del search

    with arcpy.da.UpdateCursor(lyr,"mergeid") as cursor:
        for row in cursor:
            row[0] = FID1
            cursor.updateRow(row)
            print row
    # del cursor

My question is that how would I write the script so that the search only looks for selected polygons instead of every polygon in the the map document.


texts: Cannot provide inputs for saga:difference using qgis_process-qgis (command line qgis)

texts:
I'm using the following command to find the difference between two Sentinel 2 rasters using saga:difference tool in Powershell (works fine in standard QGIS3):
./qgis_process-qgis run saga:difference --A=1.jp2 --B=2.jp2 --RESULT=test

However I get the following error:
****************
Warning: this algorithm contains known issues and the results may be unreliable!
****************

----------------
Inputs
----------------

A:      1.jp2
B:      2.jp2
RESULT: test

ERROR:  An error was encountered while checking parameter values
        Could not load source layer for A: 1.jp2 not found

Other tools such as native:buffer work fine.


texts: How to load html file with leaflet map using PyQGIS

texts:
I have a simple HTML file with Leaflet map layers. I want to load these layers in QGIS Map Canvas with the ability to choose one of the layers.

How can I load this HTML and show one of map layers with using python (QGIS plugin)?

This is an example of my HTML:

<!DOCTYPE html>

    <title>Leaflet Map Template</title>`
    <link rel="stylesheet"
          href="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" />


texts:
Does QGIS support NAD83 Montana State Plane Zone, US Foot?
I have vector data in a CSV file that I would like to display, but the coordinate system is MT83F.  I don't see a matching CRS.  Is there a user defined CRS that is appropriate?


texts: QGIS Coordinate Reference System Problem

texts:
I'm trying to load up a set of shapefiles in QGIS, to undertake a points in polygon analysis. I load up the points in OSGB 1936 / British National Grid (EPSG: 27700) (the same as the polygons), however they do not appear. I've enabled 'on the fly' transformations. See a screenshot below with no visible points layer:

I've opened another project with the points in WGS 84 and then saved them in the transformation required (EPSG: 27700). Still no luck.

I then opened them in Arcmap and they worked fine, see below:

I'm all out of ideas, are there any other suggestions for how to get this working in QGIS?

EDIT: See the requested metadata for the points layer below:

General
Storage type of this layer
ESRI Shapefile
Description of this provider
OGR data provider (compiled against GDAL/OGR library version 2.0.2, running against GDAL/OGR library version 2.0.2)
Source for this layer
C:/Users/ejo31/Desktop/cell.ID.locations/open.cell.id.27700.shp
Geometry type of the features in this layer
Point
The number of features in this layer
1332599
Capabilities of this layer
Add Features, Delete Features, Change Attribute Values, Add Attributes, Delete Attributes, Create Spatial Index, Fast Access to Features at ID, Change Geometries, Presimplify Geometries, Presimplify Geometries with Validity Check
Extents
In layer spatial reference system units
xMin,yMin -99.7747,-44.7261 : xMax,yMax 174.781,76.6116
Layer Spatial Reference System
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs


texts: Using xlsx data to clip shapefile in ArcGIS Desktop

texts:
I have an Excel table with just one column, with 1150 records, each of them a species name. I also have a shapefile with 24,000 records, each of them a species name and its geographic range. I need to clip out of the shapefile the 1150 records in my Excel table.

I thought I could just use ModelBuilder’s convenient row iterator to do this, but while the row iterator works fine the clip tool will not accept the result as the clip feature input. I am now trying to convert my Excel xlsx file into a shapefile, but since it doesn't have any geographic locations, it's just a list of names, this seems to be impossible.


texts: Get Graphics from KML Layer in ArcGIS JS API

texts: I just started working with the ArcGIS Javascript API and I was wondering how I would go about getting an array of Graphics from a KML layer. I do not know where to begin.


texts:
I am new to using MODIS data and was hoping to get some help starting out.  I have been able to download the hdf file but viewing it in R has been a problem for me.  What I have found so far, and tried, is using the gdalUtils package.  What I have done so far is this

gdal_translate('MOD11A1.A2012027.h11v04.005.2012028122822.hdf', 'test.tiff',
sd_index=1)

This gave me a tif image that I then put back into R using

data <- readTIFF('test.tiff')

This gave me a matrix in R that I could work with, which is what I ideally want.  The only issue I had was that the values in the matrix were all decimals and when I did an image plot the each axis went from 0 to 1 instead of the proper lat and lon.  the file I am using is LST in kelvin so these values are not correct. I know there is gdalwarp that does the reprojection but I am unsure how to use it.  The inputs for s_srs and t_srs I was unsure of.

Can anyone point me in the right direction as a beginner.  I have been reading the Steve Mosher blog so I have MRT and OSGeo4w installed plus all the packages he suggested.


texts: how to get a low level road network from openstreetmap data?

texts:
I've extracted the road network of a city from openstreetmaps (filetype: osm). One part of the road networks look like this:

However, there are too many details in this data.
I'd like to represent each road with just one linestring, just like this:

Is there any convenient (or maybe a tricky) way to reach that?
I've searched around our site, I found some map skeletonization/genratzation may help, but I'm not quite sure about that. I've also tried the "network toplogy cleaning" function in OpenJUMP, but it didn't work.
I would very happy to see a software/library can help me on this.

Or, is there any alternative way to extract rough(low level) road network of a city from other map source?


texts:
I have been having trouble loading 30m SRTM files into blender.
I download the SRTM (DEM) 30m files from the earthdata (NASA) website and I would like to prepare them in QGIS (making mosaic, clip, reprojection) for reading in the blender, but I have not been able to. How can I do this? If possible, with step-by-step illustrations, or even a video teaching how to configure an SRTM to open in blender.


texts: Extract date, time and GPS position linked to road crossing

texts:
I am using QGIS. I have a large dataset of approx a hundred animals and I want to determine the number of successful road crossings and extract the date, the time and the two GPS coordinates before and after the road crossing. Please see screenshot attached that shows the example for one animal.

I have intersected the vector line layer from the road network (in red) with the vector line of the animal movements that links its GPS locations and obtained a layer of points where the animal likely crossed the road (round white coloured dot in the screenshot).

When I create a buffer around the intersected points to include the GPS locations linked to each road crossing, I end up including points where the animal did not cross the road because they are included within the distance of the buffer.

How can I include only the GPS points linked to the data crossing? It would be one GPS location at each side of the road (in red).


texts: how to read row values instead of field values?

texts:
I use this to read the values of fields:

lyr = iface.activeLayer()

features = lyr.getFeatures()

for feat in features:
    attrs = feat.attributes()
    print attrs[1]

how can I read values of rows instead? and how can i loop through each row and print out the values like in a list for each row?


texts: CartoDB: How to generate automatic labels and then move those labels around the map

texts: I'm working with CartoDB and hoping to automatically generate labels that can then be moved around the map. I'm able to generate the labels, but it won't let me edit their placement.


texts: ArcGIS 10.2 with MatPlotLib Non-Blocking Plot

texts:
ArcGIS 10.2 ships with MatPlotLib 1.1.1 installed.  I am trying to display a non-blocking plot window.

Trials:

Launch an interactive plot in a Tk window from a toolbox.  This causes the plot to display and the toolbox dialog continues to run in the foreground.

    import matplotlib.pyplot as plt

    fig = plt.figure()
    plt.plot(range(10), range(10))
    plt.show()

Changing plt.show() to plt.show(block=False) should make the plot non-blocking.  This crashes ArcMap.

Attacking this from another direction - Create a python addin-in with a single button.  That button executes the above code.  Also crashes ArcMap.

Is it possible to get a non-blocking MatPlotLib plot via arcpy?


texts:
I'm trying to create a dashboard similar to the one done in https://colab.research.google.com/drive/1NeJCn51aSPPVNli2n8W3q_aUvBHoH-CJ?usp=sharing#scrollTo=Rs2m4uSDXGwX,
where a choropleth map is displayed representing values within a geographic unit. In my case, it's zip codes and I'm color-coding the map by smell reported cases for each zip code. I'm getting the following error message:

Error: undefined Error: "incorrect data check". CARTOframes is not
able to parse your local data because it is too large. Please, disable
the data compresion with encode_data=False in your Layer class.


texts:
I'm trying to drop the .shp from my list but when I do and print, I am only seeing the last one. When I print the list from ListFeatureClasses() I am seeing all of them. Any thoughts?

Here is the code:

import arcpy

workspace = r"C:\Users\ChapmaT1\CountyFiles"
inputSpatialReference = r"C:\Users\ChapmaT1\Script\roads.shp"

#get SR of inputSpatialReference
inputSR = arcpy.Describe(inputSpatialReference)
toSpatialRef = inputSR.SpatialReference
toSpatialRefName = toSpatialRef.Name

#Set Up Workspace
arcpy.env.workspace = workspace

#List of Shapefiles
shapefiles = arcpy.ListFeatureClasses()

#print shapefiles

#Good till Here########################################

for shape in shapefiles:
    shapeName = shape.replace('.shp', '')

print shapeName

The same can be seen using a non-ArcPy example (no import arcpy needed) of:

for shape in ["test1.shp","test2.shp","test3.shp"]:
    shapeName = shape.replace('.shp', '')

print shapeName


texts:
I am a beginner to python coding in GIS. I want to calculate the total number of counts in raster attribute table. In order to do that:

I would like to build a attribute table of raster data
use a loop through attribute table using ARCPY.SearchCursor()

Therefore, I wrote a code to build an attribute table. However, I am getting an error while running the code and not able to figure out the error in the code.

Code:

import arcpy
# Create the geoprocessing object
import arcgisscripting
gp = arcgisscripting.create()

#To overwrite output
arcpy.env.overwriteOutput = True

#Set environment settings
gp.workspace = "C:/Subhasis/Test/raster-read"

#Build Raster Attribute Table
gp.BuildRasterAttributeTable_management("sti-9", "Overwrite")

Error Message:
Traceback (most recent call last):
  File "C:/Subhasis/Test/raster-read/test.py", line 15, in <module>
    gp.BuildRasterAttributeTable_management("sti-9", "Overwrite")
ExecuteError: ERROR 000049: Failed to build attribute table
Failed to execute (BuildRasterAttributeTable).


texts:
I am using QGIS 2.18 and have plotted shapefile received from surveyor. data displays ok, now im trying to display coords of point, but the data is not in the attributes table.
the xyz values appear in something called a derived table and not in the feature attributes table. when i use identify feature.
the derived feature id is also different from the attributes table id number.
i tried copying the feature attributes, but did not give me the derived data where the coord info resides.
how can i either merge these tables or export the derived attributes table.


texts:
I have census data that uses the widely-established CSPro software.  It also comes with a Map viewer and the data files are .map and .mdf extensions.

Does anyone know how to convert these into an ArcGIS shapefile?  CSPro gives instructions on how to convert an ESRI shapefile to CSpro format but not the other way round.  It must be possible, but I do not know how.  I've tried OGR2GUI but without success - it does not recognize .map file extensions.


texts:
I have several scripts written in R language with package "sp" for Qgis. However I am stuck in the old versions (2.18) because whenever I use this same package in current versions 3.14 the script has an error. What alternatives could i do to work around this problem?
I didn't write the scripts, I only use them.
I've been researching and the "sp" package has been replaced by "sf" but I don't have the cognitive ability to rewrite all the scripts.
This is 1 of the 12 scripts I use in version 2.18.
When I try to run in version 3.14 I get this error.
##Determina_zonas_de_manejo=group
##Raster_de_Entrada=raster
##Poligono_mascara=vector
##Quantidade_de_Celulas_para_Remover_das_Bordas=number 0
##Raster_Corte_por_Mascara=output raster

library(raster)
library(sp)
require(rgdal)
library(rgeos)

raster_in  <- Raster_de_Entrada
pol_mask <- Poligono_mascara
rm_bord <- Quantidade_de_Celulas_para_Remover_das_Bordas

#raster_in <- brick('C:/Dados/ManejoPrecisao/UGD_geodata/ndvi.tif')[[1]]
#pol_file <- 'C:/Dados/ManejoPrecisao/UGD_geodata/talhao01.shp'
#pol_mask <- readOGR(dirname(pol_file), strsplit(basename(pol_file), '[.]')[[1]][1] )
#rm_bord = 1

cel_tam <- res(raster_in)[1]
rm_bord_tam <- rm_bord * cel_tam

rm_bord_buf <- gBuffer(pol_mask, width= -rm_bord_tam)

if(rm_bord == 0){
  crop_in <- crop(raster_in, extent(pol_mask), snap = 'out')
  raster_out <- mask(crop_in, pol_mask)
  Raster_Corte_por_Mascara = raster_out
} else{
  crop_in <- crop(raster_in, extent(rm_bord_buf), snap = 'out')
  raster_out <- mask(crop_in, rm_bord_buf)
  Raster_Corte_por_Mascara = raster_out
}