texts: Linux QGIS 2.18.4 GRASS7 tools not working

texts:
GRASS tools are showing & GRASS 7 shows as activated in the Providers menu, but I get the "Missing dependency, GRASS isn't configured properly" message. How do I configure it properly?

GRASS 7.2 standalone shows in the main menu but doesn't run - a window flashes up momentarily.

This is a recent install of 2.18 after removing 2.16.
I did this to install:

add-apt-repository ppa:ubuntugis/ubuntugis-unstable
then
apt-get update # apt-get install qgis apt-get update # apt-get install qgis apt-get update # apt-get install qgis
sudo apt-get install qgis python-qgis qgis-plugin-grass

I've read threads on here but Cannot run Grass algorithms in QGIS 2.18 on Ubuntu 16.04 may be close although I haven't mixed sources - how do I find what version of GRASS - both standalone & plugins is (partially?) installed?

In a previous post, AndreJ stated:

  You ran into gdal-abi-2-1-2 but it is not installable because that is only available for xenial. Trusty is still at gdal-abi-2-1-0. You might need to upgrade GDAL before installing QGIS.

I've overcome the xenial problem (ie using trusty), is the gdal version still a problem?  How do I find out?

FWIW, SAGA is installed and working.


texts: How is attribute data in DBF file tied to shapefile location data in SHP file?

texts:
I opened a DBF file in Excel and I'm looking at attribute data from a shapefile export.

How is the row of attribute data in the DBF file connected to it's corresponding locational data in the SHP file?

I don't see anything that looks like a key.


texts: Any tool to view .mvt files content in human readable form

texts: I have tried with protobufeditor but didn't succeed to open it. I need something similar to check for errors in existing files, and binary data is too hard to read when file is large. I have downloaded vector_tile.proto file and compiled it with google libs.


texts:
I try to rotate a custom Marker icon inside Mapbox gl-js. All the commented attempts not working.
Are we absolutely need to go with geojson sources to rotate an icon and a custom one that is not able to fit on Maki svg standard! ?

$.get(mapyim, function (resp) {
  resp.features.forEach(function(marker) {
  var el = document.createElement('div');
  el.className = 'myMarkerClassCss';
  // el.style.transform = 'rotate(' + marker.properties.angle + 'deg)';
  el.addEventListener('click', function() {
    window.alert(marker.properties.title + '/' + marker.properties.date);
  });

  // add marker to map
  new mapboxgl.Marker(el, {
    // rotate: marker.properties.angle,
    offset: [
      -26 / 2,
      -33 / 2
    ]
  })
  .setLngLat(marker.geometry.coordinates)
  .addTo(map);
 });
});


texts: Showing feature count for layout extent on legend using QGIS?

texts:
I'm using QGIS 3.6.1 on Win10.

I have a print layout with a map of an area within my total work data set.

I add an item to the print layout legend with the summation button on for "Show feature count for each class of vector layer." I expect this to show the count of the features visible within the layout map extents. Instead, it shows the total count for the entire data set.

Is there a way to just count the data within the map layout extents and show this on the legend?


texts: Mosaic from two satellite images using python

texts:
To generate mosaic of the images I used gdal_merge as shown below.

python gdal_merge.py -o out.tif -n 0 -a_nodata 0 in1.tif in2.tif

I find a clear demarcation between input 1 and input 2 image as shown in the image, but I am expecting the final output with no demarcation or smooth surface in the output. How to rectify it?


texts: Leaflet - draggable marker and coordinates display in a field form

texts:
I have to do a draggable marker and its coordinates should be displayed in fields. It will be a part of a contact form in PHP. I created a draggable marker, help me what to do now.

var marker = L.marker(new L.LatLng(53.471, 18.744), {
draggable: true
}).addTo(map);

http://jsfiddle.net/xTh5U/

Here is example in Google Maps API, I need the same in Leaflet.


texts: Buffering Bounding Boxes and/or Polygons QGIS

texts:
I just created a Bounding Boxes layer using a layer of oddly shaped polygons as the input.

The Bounding Boxes layer has a bunch of symmetrical rectangular polygons on it now, all oriented parallel with the X axis.

I need to buffer all of these rectangular polygons now.
When you buffer a line using the buffer tool, and select Flat end style, the tool doesn't buffer the ends of the line.

The buffer result is flush with the end of the input line.

I need to do the same sort of thing with these rectangular polygons.

I want the buffer only to apply to the North and South sides of the rectangular polygons and not the East and West sides.

So the rectangular polygons will get taller when the buffer is applied, but not wider.

Is there a tool that will do this?

When I use the buffer tool it always applies the buffer to all sides of the polygons no matter what the settings are.


texts: Creating Generic Transverse-Mercator projection in ArcPy with ArcMap

texts:
I'm attempting to create a Transverse-Mercator projection in arcpy. However I keep running into an error trying to set the latitudeOfOrigin, even when the documentation explicitly says this is a required attribute. I am using ArcMap 10.5.
I've tried creating the projection from scratch and tried simply modifying the parameters on an existing spatial reference. In both cases I get the following error when trying to set latitudeOfOrigin:

NameError: The attribute 'latitudeOfOrigin' is not supported on this
instance of SpatialReference.

What is the proper way to create this projection?
#create from scratch
def test_spatial2(self, lat, lon):
    sr = arcpy.SpatialReference()
    sr.projectionCode = 43006
    sr.projectionName = 'Transverse_Mercator'
    sr.type = "Projected"
    sr.linearUnitCode = 9001 #meters
    sr.gcs = arcpy.SpatialReference(4326) #WGS 84
    sr.scaleFactor = 1
    sr.centralMeridian = lon
    sr.centralMeridianInDegrees = lon
    sr.latitudeOfOrigin = lat
    sr.falseEasting = 0
    sr.falseNorthing = 0
    sr.create()

#create from UTM
def test_spatial3(self, lat, lon):
    sr = arcpy.SpatialReference(32610) #WGS84 UTM North Zone 10
    sr.centralMeridian = lon
    sr.latitudeOfOrigin = lat
    sr.falseEasting = 0
    sr.falseNorthing = 0
    sr.create()


texts: Convert Parallel Polylines (edges of pavement) to Polygons

texts:
I would like to find a way to convert an Edge of Pavement layer into polygons that would represent the total paved area. It's a large shapefile, covering some 100 sq miles or so, so doing anything manually is out of the question.

I've tried the Feature to Polygon tool but got the result below. The blue lines are the edges of pavement and the tan represents the polygons that were created. Not only did the tool not fill in the spaces between the lines, but it left large gaps throughout the rest of the area as well (in white).

I also tried downloading and using the Split Polygons by Lines tool found here, placing a large square polygon beneath all of the lines. Unfortunately this tool has been processing for hours and doesn't appear to be working. Perhaps I have too may lines processing at once and I need to split it into smaller segments.

Just wondering if road polygons are something anyone else has created. I ultimately am looking to make a layer of impervious surface by combining these polygons with a buildings layer.


texts: query to get predicted location after traveling on pgrouting route for some hours, postgis 2, postgresql 9.3

texts: I am using Postgis 2.0 on Postgresql 9.3 on Ubuntu 14.04 x64. I am trying to find an efficient way to fulfill this routing request. So I am able to get the shortest route using the source and target locations. However once I have that route, I would like to know where a car or person will be x hours after starting from the source location. I know that the cost function has the travel time in hours, so I should be able to find the longitude and latitude location after traveling for some number of hours. However, I am not familiar enough with pgrouting to know the most efficient way to run this kind of query. Any suggestions?


texts: Improving PostGIS performance for spatial intersection

texts:
I have a fairly large table with line geometries, representing roads (about 900k rows) and a table with polygons representing land use (about 8k rows). I want to know what type of land each road passes through. For most it is only one, but some pass through two different polygons.
I've written the following statements and tested them on a subset of the road data (around 4.5k), and they work (take around 33 seconds to create the view). I've tried running it with the full table, but my computer crashed after 2 hours.
Is there any way to optimize the query, so it runs a bit faster?
CREATE MATERIALIZED VIEW ways_lc AS
    SELECT w.osm_id, string_agg(l.label_modified, ', ') AS landcover FROM waysdk w
    JOIN land_cover l ON ST_Intersects(w.geometry, l.geometry) GROUP BY osm_id;

UPDATE waysdk w SET landscape = landcover
    FROM ways_lc lc WHERE w.osm_id = lc.osm_id;

There are indexes and spatial indexes on both tables. The reason I'm creating a materialized view is to also create indexes on the view, but it didn't make any difference from performance.
I have tried looking into PostgreSQL/PostGIS performance issues, but have mainly encountered solutions for very large and complex geometries, which isn't the case here.
Running EXPLAIN for the query returns this (which I to be honest have a hard time interpreting):
HashAggregate  (cost=15696875.24..16296475.21 rows=958763 width=40) (never executed)
    Group Key: w.osm_id
  ->  Gather  (cost=1000.15..12379366.37 rows=21489936 width=27) (never executed)
          Workers Planned: 2
        ->  Nested Loop  (cost=0.15..10229372.77 rows=8954140 width=27) (never executed)
              ->  Seq Scan on waysdk w  (cost=0.00..168900.85 rows=399485 width=142) (never executed)
              ->  Index Scan using idx_land_cover_geometry on land_cover l  (cost=0.15..25.17 rows=1 width=954) (never executed)
                      Index Cond: (geometry && w.geometry)
                      Filter: st_intersects(w.geometry, geometry)


texts: Edit multiple vertices of different features simultaneously

texts:
I am using ArcGIS 10.2.2 and would like to edit multiple vertices of different features simultaneously.  For example, I have a man hole with three laterals - is it possible to select all four features and move the end points of the laterals and manhole node at the same time.  See sketch below - move manhole and lateral end points to blue "+".  I have tried using topology but I am only able to get it to move the manhole and one line segment of the three.

I need to use this ability to adjust thousands of such feature combinations to newer high accuracy imagery I just received.  Moving each of the multiple vertices one at a time will be painfully slow.


texts: How to monitor changes in layers panel for QGIS?

texts:
I'm writing Python script in QGIS 2.18 to alert the user when multiple layers with different projections are added to the canvas. My codes are shown below.

from qgis.gui import QgsMessageBar

input_layers = []
layers = iface.legendInterface().layers()
for layer in layers:
  if layer.type() == QgsMapLayer.RasterLayer:
     input_layers.append(layer)

pcs_list = []
for input in input_layers:
  pcs_list.append(input)
  pcs_str = list(set([pcs.crs().authid() for pcs in pcs_list]))

text = "You have multiple layers with different CRS in this project: " + str(pcs_str).translate(None, "[]u")

if len(pcs_str) >1:
    iface.messageBar().pushMessage("Alert",text, level=QgsMessageBar.WARNING,duration=20)
else:
  pass

This works fine when I run the script via the Python console. However I'd like this alert to pop-up automatically when new layers (with different projections) are added to the layers panel. I'm aware that layer objects in the layers panel can be accessed through the the layer tree root. Is there any method to monitor changes in the layer tree root?


texts: What geospatial table and indexing functions are suitable for a borehole database design for large volume of serialised constant increment data?

texts:
Pretty common use case I suspect, but want to make sure I optimise from the get go. PostGIS I'm suspecting is the weapon of choice, and one way I'll be interrogating data is QGIS. Large part of data set originally sourced from LAS files.
Raw Data comes as typically 3 tables:

collars/wellhead - point data of well head xyz ,
wellpath - linestring data, dip/azimuth surveys at regular intervals, puts well/borehole in 3d
space
geophysics - 1d data, regular intervals along welpath like M values, up to say 15 parameters per depth step. Typically all 20 params at same intervals.

Size is an issue (well, for me). Given most of the data is by definition at regular intervals, say 1cm, I thought there must be some standard functions that could make the storage retrieval efficient. I'm aware of things like ST_LineLocatePoint

collars - 3000 xyz points
wellpath - 30 dip/azi points per collar, regular 10m intervals
geophysics - 50000 rows per collar, regular 1cm intervals (CWS Las 1.0), 20 single precision observations per intervals.

A typical csv containing the above weighs in at 10gb. Nice to have appropriately indexed, hence the question. Its the regular 1cm interval data that is the bulk that I'd like to have indexed appropriately.


texts:
Does anyone know a way to save my username for PostGIS in QGIS, so I'll always be logged in without the need to enter my username and password everytime?

Anytime I start new project in QGIS, I need to log in and enter my username and my password and I wonder if there is any way that it will always be logged in and I'll not have to insert my username and password.


texts: Get Geometry object from Envelope in NetTopologySuite

texts: I've built some code to take multiple geometries and merge them into one envelope, but now I need to convert that envelope into a polygon/IGeometry shape that represents that area.  Is there a way to do this?


texts:
I was recently introduced to the Avenza PDF app for iPad and I've been using it in the company I work for to help our field crews find their way around, and record waypoints and tracks. The app takes georeferenced PDF's and uses the iPad GPS to show a field person exactly where they are. It works very well and the field guys love it. The only problem is the iPad isn't synchronised with our Office media, which is all PC and Microsoft products, plus our field crews are now carrying an iPad plus a laptop which they need for reporting duties. If we were able to purchase the Surface then we could do away with the iPad and the laptop and keep the Surface.

So I've been asked to research the MS Surface. Does anyone out there know if the Surface has decent (if any) GPS capabilities? Or whether there is an app that would enable me to download custom made maps into the tablet? I've found Microsoft MapPoint 2013, but I have no experience with this program.

Basically I'm looking for an Avenza app equivalent for the Surface.


texts: How To Use Ampersand (&) Symbol For a Label In Symbology 10.2?

texts:
I am using unique values to symbolize my data. I am changing the labels as needed so they are reflected in my table of contents. I want to label one category "Oil & Gas". However when I use the ampersand (&) in the label, it becomes an underscore "_".

Any thoughts on how to keep the ampersand?


texts: Automatically populate attribute table based on feature location within another feature?

texts:
I have two datasets. One of all the suburbs in a particular state. The other of all the housing lots across the state. I want to add a field to the lot feature class' Attribute table which can then use some tool to automatically populate which suburb each lot falls within.

Is there a simple process to do this?


texts:
Is there any alternative to Google's S2 geometry library? It looks it was updated 4 years back and hardly any documentation.

I need a library to divide cities into zones of the area I need. S2 looks good solution, but its documentation is very poor. I have even no clue how to use it. I am looking at its Java library which will be compatible with the project I am doing.

If no alternative, any pointers how to start on S2? Like, how do I create cells with IDs of the resolution I need? How do I use them to search for the proximity given a point??


texts: Does ogr2ogr from the Docker osgeo/gdal support PDF layers?

texts:
I'm trying to use the Docker osgeo/gdal image to extract vector data from a GeoPDF, i.e. https://www.floridastateparks.org/sites/default/files/media/file/SFHPSP_TrailMapBike_11-27-19.pdf. It's compiled with the poppler PDF backend.
gdalinfo -mdd layers shows plenty of layers, but ogr2ogr doesn't recognize any of those layers, e.g. ogr2ogr -f "ESRI Shapefile" other.shp SFHPSP_TrailMapBike_11-27-19.pdf --config OGR_PDF_READ_NON_STRUCTURED YES Layers.Bike_Trails or --config GDAL_PDF_LAYERS ALL
Is there an option to get ogr2ogr to recognize these layers, or is something else (like a different driver or a better PDF) needed?


texts: How to speed up topogeometry layer in postgis topology?

texts:
I have converted e simple polygon layer to e a topogeometry layer in postgis topology.

And now I have 3 main layer, Edge_data (lines),nodes(points) and the polygons topogeometry (.topo).

I can open these layers in Qgis and everything is ok beside the topogeometry layer (polygon), it takes to much time to render and it is not practical.

Is this normal for this type of geometry or it can be done something to speed up things, the number of polygons is 25000.


texts: clip aerial map from Add Basemap Imagery

texts:
I wish to clip an aerial image from the "Imagery" of Add Basemap in Arcmap 10.3 and save as a sid file.
How can I do that? Can someone help me please.

Example, I am looking to clip only the state of Minnesota from the "Imagery" of add basemap.


texts: QGIS 3 geometry/topology checking problem

texts:
I have a polygonal shapefile which has one overlap, where one polygon overlaps another polygon.

Vector -> GeometryChecker (Rule: must have no overlap) finds this error:

When I use the new QGIS 3.81 Vector -> Check geometries (with check for gap enabled) No error is shown:

Now I check with Vector -> Geometry Tools -> Check  Validity , also No error is shown:

Question: I am misusing Vector -> Check geometries and Vector -> Geometry Tools -> Check  Validity or is there a bug????

Test-data are here: https://www.dropbox.com/s/yj981h8b948mz00/stuff.zip?dl=0


texts: Exporting feature class to SQLite database using ArcPy with ArcGIS Pro?

texts:
I am trying to export 5 feature classes from Geodatabase to SQLite database using arcpy script. I have ESRI script that creates empty SQLite Database all I need now is just to update the script to make it copy my features classes.

# -*- #################
"""
Tool name: Create SQLite Database
Source: CreateSQLiteDatabase.py
Author: ESRI

Create a SQLite database of the selected flavor: ST_Geometry, SpatiaLite, GeoPackage.
"""

import arcpy
import os
import sys

def create_sqlite_database(in_sqlite_location, geometry_storage):
    """ create a new SQLite database. """

    # Strip the file extension from the location.
    in_sqlite_database = os.path.splitext(in_sqlite_location)[0]

    # Verify that the target folder exists.
    if not(arcpy.Exists(os.path.dirname(in_sqlite_database))):
        arcpy.AddIDMessage("ERROR", 695)
        sys.exit (1)

    if arcpy.Exists(in_sqlite_location):
        if arcpy.env.overwriteOutput == False:
            arcpy.AddIDMessage("ERROR", 159, in_sqlite_location)
            if __name__ == '__main__':
                print (u'ERROR ' + arcpy.GetIDMessage(159))
                sys.exit(1)
            else:
                raise arcpy.ExecuteError (arcpy.GetIDMessage(159))

    # Create the SQLite dataset
    arcpy.gp.CreateSQLiteDB(in_sqlite_location ,geometry_storage)

if __name__ == "__main__":
    # Get the parameters from the script tool
    create_sqlite_database(arcpy.GetParameterAsText(0),
                   arcpy.GetParameterAsText(1))

The problem is that I cannot use feature class to feature class because in basic licence I cannot copy feature from geodatabase to sqlite database.

I am using ArcGIS Pro 1.3 Basic.


texts: Using where clause for arcgis analysis in ArcPy?

texts:
I am trying to write a python code for select analysis.

FIPS: column in feature nameofinputfeature whose values needs to be compared with OFIPS

OFIPS: variable whose value is known each time the code is run.

I am having a problem in syntax of the {where clause}

my code is:

arcpy.Select_analysis("nameofinputfeature", "C:\\a.shp", '"FIPS" = %OFIPS%')


texts:
I am trying to create a postcode search using a dialog box but if I press "cancel", it automatically shows ALL of the addresses within the area and doesn't just close the box. Not sure why!

def pc_search(self):
#Get the postcode to search for
qid = QInputDialog()
title = "Enter Postcode"
label = "Postcode:"
mode = QLineEdit.Normal
default = ""

pc, ok = QInputDialog.getText(qid, title, label, mode, default)
print pc

# Open the address data table and search for the postcode
uri = QgsDataSourceURI()
uri.setConnection("ictdsgdev5", "5432", "gistest", "postgres", "gistest1")
uri.setDataSource("public", "OSADDBASE_LEICS", "geom", "POSTCODE like '" + pc.upper() + "%" + "'")
print uri.sql()
vlayer = QgsVectorLayer(uri.uri(), "Address", "postgres")
print vlayer.isValid()

if vlayer.isValid():
    # Add the address layer to the list of layers
    QgsMapLayerRegistry.instance().addMapLayer(vlayer)

    # Select the found addresses and zoom to them
    selectList=[]

    for feature in vlayer.getFeatures():
        geom = feature.geometry()
        selectList.append(feature.id())

    vlayer.setSelectedFeatures(selectList)

else:
    #MessageBox displaying no address found
    QMessageBox.information(None, "Information", "No address(es) found for " + pc.upper())


texts: How to turn all lines into one direction?

texts:
How can I turn all features of my polyline (postgis and/or shapefile) into one direction - it doesn't matter what direction. I know how to switch the direction of  lines but I have no idea how to bring all lines (I have thousands of them) into one single direction.

EDIT:

A line can have only two directions, the one or the other like 0 or 1. I want them all to have the direction 0 OR 1. For example the left line in the image is in direction 1 the middle line is in direction 0 and the right one is in direction 1. Now I want the middle line to direction 1 OR the right and left line to direction 0. So that all lines have the same direction 1 OR 0.

geom;given_direction;wanted_direction
line;1;1
line;0;1
line;1;1
line;1;1
line;0;1
line;1;1


texts:
This is my code, and it generates this error. I'm wondering how to fix it.
Can't convert a computed Geometry to GeoJSON. Use evaluate() instead.
    in createNDVIComposite(), line 30
    in <global>, line 57
    in <global>, line 59

var roi = roi.geometry();
//==================================//
// LANDSAT NDVI COMPOSITE CREATION //
//==================================//

// NDVI Functions
/*These 2 functions are used to calculate the NDVI for Landsat 5/7 and Landsat 8, respectively.*/
var getNDVI57 = function(image) {
  var ndvi = image.normalizedDifference(['B4', 'B3']);
  return(ndvi);
};
var getNDVI8 = function(image){
  var ndvi = image.normalizedDifference(['B5', 'B4']);
  return(ndvi);

};

// Main function
var createNDVIComposite = function(){
  var yearrangeStart = 1999;
  var yearrangeStop = 2018;
  var imgList = [];
  for(var loopYear = yearrangeStart; loopYear <= yearrangeStop; loopYear +=1){
    var start = ee.Date.fromYMD(loopYear, 1, 1);
    var end = ee.Date.fromYMD(loopYear, 12, 31);

    //Landsat 8
    var l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
      .filterBounds(roi)
      .filterDate(start, end)
      .filterMetadata('CLOUD_COVER', 'less_than', 80)
      .map(getNDVI8);

    //Landsat7
    var l7 = ee.ImageCollection("LANDSAT/LE07/C01/T1_SR")
      .filterBounds(roi)
      .filterDate(start, end)
      .filterMetadata('CLOUD_COVER', 'less_than', 80)
      .map(getNDVI57);

    //Landsat5
    var l5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR")
      .filterBounds(roi)
      .filterDate(start, end)
      .filterMetadata('CLOUD_COVER', 'less_than', 80)
      .map(getNDVI57);

    //Merge collections
    var mergedCollection = ee.ImageCollection(l8.merge(l7).merge(l5))
                            .reduce(ee.Reducer.median())
                            .reduce(ee.Reducer.percentile([90]));

    var finalOutput = mergedCollection.clip(roi).rename(loopYear.toString());

    imgList.push(finalOutput);
    //Generate filename for export
    var filename = ("Annual_NDVI_"+roi+"_").concat(loopYear.toString());
  }

  // Visualization parameters.
  //var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
  var args = {
    crs:'EPSG:3857',
    dimensions: '300',
    region: roi,
    min: -1,
    max: 1,
    //palatte: [ndviParams],
    palette: ['blue', 'white','green',' red'],
    framesPerSecond: 2,
  };

  var yearImgCol = ee.ImageCollection.fromImages(imgList);

  var thumb = ui.Thumbnail({
    // Specifying a collection for "image" animates the sequence of images.
    image: yearImgCol,
    params: args,
    style: {
      position: 'bottom-right',
      width: '320px'
    }});
  Map.add(thumb);

}

//print(ui.Thumbnail(Class.select('classificacao'), videoArgs));

//-----------------------------------//
// MAIN FUNCTION EXECUTION //
//-----------------------------------//
var comp = createNDVIComposite();


texts: Issue with circular patterns showing up after orthographic projection

texts:
For a project I'm working on, I'm trying to project a GOES-16 image to an orthographic projection centred at the continental United States. In a later stage, I found these odd circular patterns appear in a result derived from an image I had reprojected to this orthographic projection. I have been able to trace this back to the projection, and can reproduce the patterns by computing the number of original GOES-16 pixels that are part of each pixel in the final orthographic reproduction.
Why is this? And if possible, how do I resolve it (besides choosing a different projection)?
I have included a (hopefully) working example below.
import numpy as np
import cartopy.crs as ccrs
import pyproj
import pandas as pd
import matplotlib.pyplot as plt

# Extent in orthographic projection coordinates
EXTENT = [-2373970.2054802035, 3638399.906890716, -2388159.0372323524, 1620933.528091551]

# Orthographic projection
p_ortho = pyproj.Proj('+proj=ortho +ellps=GRS80 +lat_0=39.8283 +lon_0=-98.5795 +sweep=x +no_defs')

# GOES-16 nominal satellite height
h = 35786.0234375

# GOES-16 projection
proj_goes = pyproj.Proj(proj='geos', h=h*1000, lon_0=-75.2, sweep='x')

# GOES-16 scanning angles
x_goes = np.arange(-0.151844, 0.151844+5.5998564e-05, 5.5998564e-05)
y_goes = np.arange(0.151844, -0.151844-5.5998564e-05, -5.5998564e-05)

# Coordinates for lines of constant x and y scan angle
xx_goes, yy_goes = np.meshgrid(x_goes[::100], y_goes)
lon_x, lat_x = proj_goes(xx_goes*h*1000, yy_goes*h*1000, inverse=True)

xx_goes, yy_goes = np.meshgrid(x_goes, y_goes[::100])
lon_y, lat_y = proj_goes(xx_goes*h*1000, yy_goes*h*1000, inverse=True)

# Create x-y grid to reproject
xx_goes, yy_goes = np.meshgrid(x_goes, y_goes)

# Project scanning angles to corresponding longitude latitude
lon_goes, lat_goes = proj_goes(xx_goes*h*1000, yy_goes*h*1000, inverse=True)

# Filter for CONUS domain (roughly)
lon = lon_goes[(lon_goes > -136)*(lon_goes < -45)*(lat_goes>10)*(lat_goes <55)]
lat = lat_goes[(lon_goes > -136)*(lon_goes < -45)*(lat_goes>10)*(lat_goes <55)]

# Now project longitudes, latitudes to orthographic projection
xx_p, yy_p = p_ortho(lon, lat)

# Use dataframe for calculation
df = pd.DataFrame({'x': xx_p, 'y' : yy_p})

# Min and max coordinates in orthographic projection
[x_min, x_max, y_min, y_max] = EXTENT

# Number of bins
x_bins = 300
y_bins = 200

# Calculate bin size
x_size = (x_max-x_min)/x_bins
y_size = (y_max-y_min)/y_bins

# Calculate bin number for xx_p, yyp
df['x_bin'] = np.floor(((df['x'] - x_min)/x_size)).astype(int)
df['y_bin'] = np.floor(((df['y'] - y_min)/y_size)).astype(int)

# Filter because I was lazy in defining the CONUS domain (?)
df = df[(df.x_bin >= 0) & (df.x_bin <=x_bins-1) & (df.y_bin >= 0) & (df.y_bin <=y_bins-1)]

# Reshape results into 2D array for display
res = np.zeros((y_bins, x_bins))
grouped = df.groupby(['x_bin', 'y_bin']).count().reset_index()
res[y_bins-1-grouped.y_bin.values, grouped.x_bin.values] = grouped.x.values

# Show results
fig = plt.figure(figsize=(10, 5), dpi=300)
proj = ccrs.Orthographic(central_latitude=39.8283, central_longitude=-98.5795)
ax = fig.add_axes([0, 0, 1, 1], projection=proj)
ax.set_extent(EXTENT, proj)
ax.set_aspect('auto')
cf = ax.imshow(res, origin="upper", extent=EXTENT)

# Lines of constant scan angle
for i in range(lat_x.shape[1]):

    ax.plot(lon_x[:,i], lat_x[:,i], linewidth=0.25, c='w', transform=ccrs.PlateCarree())

    ax.plot(lon_y[i,:], lat_y[i,:], linewidth=0.25, c='w', transform=ccrs.PlateCarree())

ax.coastlines()
fig.colorbar(cf, label="Number of GOES-16 ABI pixels per 100 orthographic projection pixels")


texts: How to display vectors in arcgis json format in openlayers?

texts:
For this project I am working on, I would like to display an ArcGIS Online FeatureService (JSON) in a webmap that I have built using OpenLayers.

I am able to read the JSON from the service, but I can't seem to 'convert' into a vector layer.

My code is the following:

var data = ($.ajax({
    url: '...',
    success: function (data) {
    console.log(data);
    }

    }));

   var geojson_format = new OpenLayers.Format.GeoJSON();
   var vector_layer = new OpenLayers.Layer.Vector();
   map.addLayer(vector_layer);
   vector_layer.addFeatures(geojson_format.read(data));

This unfortunately produces the following errors:

TypeError: g is null (FF)

TypeError: 'null' is not an object (evaluating 'g.layer=this') (Safari)

Uncaught TypeError: Cannot set property 'layer' of null (Chrome)

Does anyone know how to fix this?


texts: Remove polygons having less than 4 coordinates in R

texts:
I have a multiPolygon shapefile of buildings in a particular city. I read it in with readOGR and and transform it to my desired planar projection with the following code:

buildings_data<-readOGR("buildings.shp")
buildings<- spTransform(buildings_data, CRS("+proj=utm +zone=45 +a=6377299.151
                                        +b=6356098.145120132
                                        +towgs84=295,736,257,0,0,0,0 +units=m +no_defs "))

It gives me the following output along with transforming:

There were 19 warnings (use warnings() to see them)

All the warnings are:

In Polygon(coords = crds) : less than 4 coordinates in polygon

So how do I remove these polygons with less than 4 coordinates?


texts: Any restrictions on the use of ArcGIS REST API?

texts:
I'm building a PHP script which ends up with an array of latitude, longitude pairs (around 2,500), and I want to be able to reverse geocode them into full adresses, then recording them on my database.

I ran a search on Reverse Geocoding and have narrowed my options to Google and ESRI. I'm leaning towards ESRI because they have a REST API, whereas Google works via Javascript, which would give me much more trouble to implement.

But before I start coding I was curious whether the ArcGIS REST API has any restrictions at all: is it free or do I need to subscribe? Do I need an API key? If so, where do I get it? Are there restrictions regarding requests per minute / hour / month? I've looked around their docs, but can't seem to find anything non-technical.


texts: Performing centrality analysis with sDNA using ArcGIS Pro

texts:
I am trying to perform centrality analysis with sDNA with ArcGIS Pro (ver. 2.7.1). However, I have failed to run the Preparation & Integral Analysis in sDNA with the same reason.
What is the problem and What should I do?


texts: Finding overlapping area of vector and raster layers in QGIS

texts:
I'm working in QGIS with 3 different layers:

A raster layer indicating 0 for non-rice and 1 for rice area
A raster layer indication 0, 1, 2 for different hazard levels
A vector layer containing the different municipalities. It shows a map with the municipal boundaries and in the attribute table the attribute NAME_2 provides the municipality name.

What I eventually need is for each municipality the total rice area, and the rice area situated in each hazard level. Thus 4 values for all municipalities.

I'm unsure how to do this as there are vector and raster formats (converting them to the same format didn't give the desired result) and as I will need the calculation for each municipality individually.


texts: Displaying WMTS tiles over Google Map (iOS app) using capabilites XML as reference?

texts:
I have written an iOS app which can display a number of different layers on a Google Map. I am trying to add another layer which works fine on ArcGIS but I can't get it to work in my app.

There are no API keys or logins or anything required by the service.

Here is the capabilities reference XML.

Here it is working on ArcGIS Online:

Based on the information in the capabilities XML I have compiled this template URL:

https://gisbasemap.ecan.govt.nz/arcgis/rest/services/Imagery/Canterbury_Imagery_1990_1994/MapServer/WMTS/tile/1.0.0/Imagery_Canterbury_Imagery_1990_1994/default/default028mm/{z}/{x}/{y}

And here is a URL at a location and zoom that I would expect to return a tile:

https://gisbasemap.ecan.govt.nz/arcgis/rest/services/Imagery/Canterbury_Imagery_1990_1994/MapServer/WMTS/tile/1.0.0/Imagery_Canterbury_Imagery_1990_1994/default/default028mm/18/256785/166329

If I can get this to return a tile in a web browser I know I can get it to work in my app.

Have I misinterpreted the capabilities XML or is this something that is not possible in this case (ie. access only from ArcGIS Online)?


texts: Calculating average time per unit in Operations Dashboard for ArcGIS?

texts:
In the Operations Dashboard, I have a table with each item a number of units and a total time of how long worked.
Number of units and total time can make in visible.
But now I want to know how much time is consumed per unit.
example,
There are a total of 63 units and the total time is 56:30 hours.
So would the number calculated 0:53 hours and want to have visible in the Operations Dashboard.


texts: Converting state name abbreviations to full names using ArcGIS Field Calculator?

texts:
I have a data table with columns that includes city names (ex: San Jose) and state name abbreviations (ex: CA). Is there a way to define the state names as their full names (i.e., CA = California) and then combine city name and state name (ex: San Jose, California)?

I have approximately 120 records. So what I think I have to do is define each state as their full names (ex: NY = New York, FL = Florida, and so on) first. I guess I could create a new field and add these to it, and then combine this new field with full state names with city names? Does that make sense or is there another way to do it? Need to achieve this using Python.

Something like the following, but it doesn't seem right:

def stateNames(): # Add this pre-logic script in the code block
  CA = 'California'
  FL = 'Florida'
  GA = 'Georgia'
  HI = 'Hawaii'
  MN = 'Minnesota'
  NY = 'New York'
  # and so on...
return stateNames()

stateNames + !City_Name! # Add this in the box below pre-logic script


texts: Controlling default value of arcpy.AddField_management

texts:
I have a large collection of shapefiles that use text fields to hold numerical values.  I have created new fields, and populated them with the numerical equivalent of the text field.  When adding the new field using:

arcpy.AddField_management(fc, numField,"DOUBLE")

all rows are assigned a zero value.  My problem is that I have both zero values and nulls in my text field (nulls are represented by a single space), and they mean something different, but when I transfer the value of the text field into the numerical field, both zero and null appear as zero.
I've tried forcing the numerical field to null when I populate it using the else clause:

with arcpy.da.UpdateCursor(fc, [strField, numField]) as cursor:
    for row in cursor:
        tstField = re.findall("[0-9]", row[0])
        if len(tstField)>0:
            row[1] = round(float(row[0]), 1)
            cursor.updateRow(row)
        else:
            row[1] = ""

but this still leaves me with zeroes for nulls.

I know that the Double field type can contain nulls, so what am I missing/doing wrong?

I'm using ArcGIS 10.1

(I realize that the way I'm checking for nulls is clunky, but looking for spaces - even when copied and pasted from the table - gave no joy.  Using re let me check for numeric values)


texts:
How can we open a 3D MAX File(3D Building model with exterior and interior elements)in Autodesk Infrastructure/LandXplorer or ESRI CityEngine?

How can we georeference a 3D vector building model?


texts:
Is there a possibility to find lower and/or upper case letters in QGIS?

I know that it is possible to convert all lower to upper: upper('hello WOrld') → 'HELLO WORLD' and all upper case to lower: lower('HELLO World') → 'hello world', but I can't find just a 'true' or 'false' output.

I have to find all upper and all lower cases to give them a different styling.


texts: Exporting multiplied raster layers as one in QGIS

texts:
I want to use QGIS to create a map from a DEM. I want color and hillshade, so I did the following:

Duplicated my DEM
For the style of the copied layer, I selected hillshade, and the multiply option
For the style of the original layer, I selected Singleband Pseudocolor, chose my ramp, classified and chose the multiply option as well

This gives me the visual output I want. How do I export this as a single layer, keeping the color and hillshade?
I want to export the result as a Leaflet map, but using the qgis2web plugin only shows the top layer (singleband pseudocolor).
What am I missing?


texts: Appending problematic polygon from File Geodatabase to Oracle ArcSDE using ArcPy?

texts:
I have a File Geodatabase that has some polygons in it. I am using append to copy them to an Oracle SDE 10.2 database. However there are some polygons that when imported have a null value for shape in the SDE. An example is shown below. Why is this happening, why is is valid in the file and not in Oracle, and more importantly, how do I get the polygon in to Oracle?

The red polygon is from the File Geodatabase, the blue is one previously imported in to Oracle. Below is a close up of the hole / intersection


texts: Difference between video captured with Google Earth and Google Earth Pro

texts: Is there any difference between video which is captured with Google Earth and Google Earth Pro? Maybe some kind of tag showing 'pro'?


texts:
I am using pgrouting on a postgis database created through osm2pgrouting. It performs very good on a limited dataset (3.5k ways, all shortest path A* searches < 20 ms).

However since I have imported a bigger bounding box (122k ways) from europe.osm the performance went down a lot (a shortest path costs around 900ms).

I would think that using A* most of those edges will never be visited as they are out of the way.

What I have done so far in an attempt to improve the speed:

Put an index on the geometry column (no noticeable effect)
Increased my memory from 8GB to 16GB
Change the postgresql memory settings (shared_buffers, effective_cache_size) from (128MB, 128MB) to (1GB, 2GB) (no noticeable effect)

I have a feeling that most of the work is being done in the C Boost library where the graph is being made so optimizing postgresql will not give me much better results. As I do minor changes to the set of rows I select for A* for every search I am a bit afraid that the boost library cannot cache my graph and has to rebuild all the 122k edges every time (even though it will only use a very limited subset every query). And I have no idea how much is spent doing that compared to the actual shortest path search.

Does any of you use pgrouting on a 122k or greater OSM dataset? What performance should I expect? What settings affect the performance most?


texts:
I have a shape file in ETRS89 CRS created with ArcGIS 10.6, prj file says:

GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

If I the create a shape file in R, with the sf library, using the same projection, epsg = 4258:

Gasolineras <- st_transform(GasolinerasWGS84, crs = 4258)
sf::st_write(Gasolineras, 'Gasolineras.shp')

and then I read it bakc into ArcGIS, it complains, because they are not in the same CRS.

But they are!

Looking at the prj file created by R, it looks a bit different than the one created by ArcGIS, but they differ only in the labels, not in the parameters:

GEOGCS["GRS 1980(IUGG, 1980)",DATUM["D_unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

It says datum "D_unknown", so it appears that ArcGIS understand that the proyections are different.

Is they any way to avoid this, and force R to write particular labels in the prj file of a shape file?


texts: Using ouput data variables in subsequent calculation within ArcGIS ModelBuilder?

texts:
Using ArcGIS Desktop 10.3.1, I have two simple polygon shapefiles, each with just a single feature.

I'm new at ModelBuilder, and have 2 tasks:

calculate the acreage within each of the two shapefiles, and then
calculate their relative size; that is, simply  divide one shapefile acreage by the other.

Here's a screenshot of my model:

I successfully completed step 1 by using the Python expression !shape.area@acres!  within each Calculate Field tool.  The correct acreage is added to the just-added field in each shapefile.

But the model returns the syntax error 000539 at step 2, which confuses me.  The ESRI documentation says that the output from each of my two Calculate Field tools in Step 1 is a variable, represented by the name of the green output oval.  So I would expect that using this expression in the Calculate Value tool would correctly give me the relative acreage: %pss_out% / %wetlandA_out%

Here's a screenshot of my Calculate Value tool:

Once I get this problem solved, I plan to populate a separate table field with the output from the Calculate Value tool.

In summary, how do I divide the output from one of the two Calculate Field tools in Step 1 by the other?


texts: Map Size not updating on map.udpateSize()

texts:
I am Using OpenLayers to publish a map. But I can't figure out the way to make the map size dynamic according to screen size. As the map remains same on every screen irrespective of the screen size. I also tried map.updatesize(). But that is not working too.
My code is
layer = new ol.layer.Tile({
    title: 'Basemap',
    baseLayer: true,
    visible: true,

    source: new ol.source.TileWMS({
        url: 'http://mlinfomaps.in/geoserver/wms',
        params: {
            'LAYERS': 'MSSDS_WS:MSSDS_BASEMAP',
            'TILED': true
        },
        serverType: 'geoserver',

        // Countries have transparency, so do not fade tiles:
        transition: 0,
        crossOrigin: "anonymous"
    })
});

var mousePositionControl = new ol.control.MousePosition({
    coordinateFormat: ol.coordinate.createStringXY(4),
    projection: 'EPSG:4326',

    target: document.getElementById('mouse-position'),
    //undefinedHTML: '&nbsp;'
});

var projection = new ol.proj.Projection({

    extent: [7970828.0670139585, 1681517.1260521673, 9127778.927138386, 2541280.8202038296]
})

function scaleControl() {
    control = new ol.control.ScaleLine({
        units: 'metric',
        bar: true,
        steps: 4,
        text: true,
        minWidth: 100
    });
    return control;
}

var sourceMeasure = new ol.source.Vector({
    crossOrigin: "anonymous",

});
var vectorMeasure = new ol.layer.Vector({

    source: sourceMeasure,
    style: new ol.style.Style({
        fill: new ol.style.Fill({
            color: 'rgba(255, 255, 255, 0.2)'
        }),
        stroke: new ol.style.Stroke({
            color: '#df1c29',
            width: 2
        }),
        image: new ol.style.Circle({
            radius: 6,
            fill: new ol.style.Fill({
                color: '#df1c29'
            })
        })
    }),
    title: 'Measure Overlay',
});

// The Map
//var overlay = new ol.Overlay
//    ({
//        element: container,
//        autoPan: true,
//        autoPanAnimation: {
//            duration: 250
//        }
//    });

var view = new ol.View({
    center: ol.proj.transform([76.7997, 18.6298], 'EPSG:4326', 'EPSG:3857'),
    zoom: 7,
    extent: projection.getExtent(),
    //minZoom:7,
    maxZoom: 12
});

var map = new ol.Map({
    target: 'map',
    controls: ol.control.defaults().extend([mousePositionControl, scaleControl()]),
    view: view,
    //fit: view.fit(),
    //overlays: [overlay],
    layers: [layer]

});
const extent = projection.getExtent()
map.getView().fit(extent);
map.updateSize();

window.onresize = function () {
    setTimeout(function () {
        map.updateSize();
    }, 200);
}

The map on different screen looks like this:
On my screen it looks like

On my colleague's screen it looks like

Can you help me get through this?


texts: How to select intersecting lines in the same shapefile

texts:
In QGIS 2.16.3, I have a series of bus routes (with different attributes), all in the same shapefile. For simplicity, they all run every 30 minutes. They often share sections of route - for example, bus #1 and #2 use the same road out of the city centre. On this section, therefore, buses run every 15 minutes.

I want to select these roads where two or more buses run (i.e. where lines in the same shapefile intersect). Along roads with two or more buses, the lines overlap perfectly.

The best I've got so far is to reduce transparency, so that roads with more than one bus appear darker than those without - but it isn't what I'm looking for.


texts: Creating field which holds stretch of length in ArcGIS Desktop?

texts:
I need to prepare a layout in arcmap, in which i have to show the table showing road stretch and width.
Example - Stretch   -   Width
         0 to 5 km  -    10m,
         5 to 12 km -    07m
Like this....
I have segment wise length in my attributes, which is showing like a below image.

For example road HHD68 starts with HHD68A at 0m and ended at 4019m, HHD68B started at 4019m and ended at 4053m, and so on up to HHD68M.


texts:
I have some geotagged photos saved on a network drive which I want to import into QGIS 3.2. I have installed the import photos plugin.  I have gone through the import photos process and selected the input folder location and output folder location.  But when I click on the path to the photo it keeps trying to take me through internet explorer to open the link (which it obviously can't open)instead of taking me to my photos on the network drive.

Anyone know how to solve this?


texts: Is it possible to extract just the East and West horizon angles for a point in ArcGIS?

texts: I have a group of observer points that I need horizon angles for (just East 0 deg and West 180 deg).  I am thinking that I should be able to somehow extract these from the Viewshed, however the process is unclear.  Does anyone know how this could be accomplished in viewshed or another application?


texts:
Not quite sure if I should even be doing this.  What I am wanting is I will be dynamically generating my SLD file and it looks like when you update the SLD in the geoserver admin it does a reload.  So I tried to do a reload using the rest API and curl and it does not appear to work.
Here is my Curl command:

curl -uadmin:password-XPOST http://localhost:8080/geoserver/rest/reload

If there is another way to clear everything so my SLD reloads that'd be awesome as well.  Just needing to get this working and am not sure why it isn't.


texts: How to add a dummy band to selected images in an image collection?

texts:
I'm using the GEE Python API. The Hansen & Al. glad alert dataset is currently broken and the 20 last images are missing the 'conf19' band. In my process I need to use it whatever the image so I'd like to add a dummy 0 value band to each image that do not contains the 'conf19' band.
I tried the following :
all_alerts  = ee.ImageCollection('projects/glad/alert/UpdResult')

# add the missing conf2019 to all images
def addConf2019(image):
    #retreive bands name
    names = image.bandNames()

    #add conf19 if missing
    if not names.contains('conf19'):
        conf19 = image.expression(
            'conf19 = 0 * A',
            {'A' : image.select('conf20')}
        )

        image = image.addBands(conf19)

    return image

alerts_cor = all_alerts.map(addConf2019)

but when I display the bands names of the imageCollection :
band_names = alerts_cor.map(lambda image:
    image.set('bandNames', image.bandNames())
).aggregate_array('bandNames').getInfo()

band_names

not a single one changed. Do you see what I am doing wrong ?


texts: Converting multipoint geometry to ONE single point using ArcGIS Desktop?

texts:
I have an ArcGIS file geodatabase with two featureClasses.

One is a point featureClass which represents the centroids of the cells of a regular grid.

The second represents the corners of each grid cell. This is a multipoint of up to four points per location, because of neighbouring grid cells.

I need to append the centroids to the corners (so that five points can later be associated (by 'nearest') with point data in each grid cell).

But ArcGIS will not permit the point geometry to be appended to the multipoint geometry.

What I do not want is to get an exploded multipoint featureClass, with up to four points per location. Perhaps this could be an interim step, with points 2,3 and 4 subsequently deleted.

What I do want is for each 'corner' to be represented by a single point.


texts: Line to point conversion with start point and end point detection?


texts: QGIS Field Calculator Error for String values

texts:
I am attempting to re-code some text data in the attribute table using the QGIS Field Calculator tool, but I keep getting an annoying error message:
"An Error occured while evaluating the calculation string: "No Root Node! Parsing Failed?"

FYI, I am using the expressions bar, this is happening while in editing mode.


texts: Downloading Sentinel-2 L2A with Google Cloud BigQuery?

texts:
I'm trying to download Sentinel-2 data in Level 2A from Google Cloud Storage.

Level 1C works perfect with the following (already discussed) method:
Download script for L1C

Now I want to do the same with L2A data but I do not know how to change the existing function query_sentinel.

Could anybody help me?


texts: Programmatically press the "select all" button in WFS capabilities of QGIS project properties with Python

texts: I would like to ask if someone knows a way to press the "select all" button in Project>Properties>QGIS Server>WFS capabilities using python, or if is there a specific command to do the same thing.


texts: PyQGIS chaining consecutive processing algorithms timing Issue

texts:
I am trying to run the merge vector layers algorithm for multiple shapefiles. Then run the join attributes by location between a merged polygon shapefile and merged multiline shapefile. All files are placed in the same directory.
The problem I am having is that once the merge vector layers algorithm is completed, the join attributes by location algorithm cannot find the merged vector layer files and triggers an error message:

List out of range for 'INPUT' in the join attributes by location algorithm.

Running the same script twice, I am able to execute the join attribute by location algorithm without errors, but running the script once causes an issue.
I tried using QgsTask and QgsTaskManager to control the timing of both processing algorithms waiting on a complete signal for the task. However, I still run into the issue. I also tried time.sleep() to try to delay the time between the 2 algorithms, but that did not work as well.
How can I get a complete signal from 1 processing algorithm to trigger the execution of the next processing algorithm? Or if that's not the issue, what the heck is going on?
  out_path = os.path.join(loc,shape_type[0])
        parameters = {'LAYERS': shape_type[1],
                      'CRS': 'EPSG:4326',
                      'OUTPUT': out_path}

        #run merging algrithm n times depending on number of shape file subtypes
        merge_layers =processing.run("qgis:mergevectorlayers", parameters)

 #join use predicate 'crosses' =6,
    #take attributes of first located feature only one feature to one feature transform = 1
    parameters = {'INPUT': SHAPEFILE_SUBTYPES[0][1][0],
                    'JOIN': SHAPEFILE_SUBTYPES[1][1][0],
                    'PREDICATE': 6,
                    'METHOD': 1,
                  'OUTPUT': out_path}

    #run spatial inner join algorithm
    processing.run("qgis:joinattributesbylocation", parameters)


texts:
I'm trying to setup Postgres 9.5 on my local machine x64 based win 8.

I've downloaded postgresql-9.5.1-1-windows.exe and also downloaded postgis-bundle-pg95x64-setup-2.2.2-1.exe

Postgres runs fine, but when trying to install postgis at the end of the instalation (when installer is attempting to install sample database) I get this error:

ERROR:  could not load library "C:/PostgreSQLx86/9.5/lib/postgis-2.2.dll": %1 is not a valid Win32 application.

I've tried with the x64 and x32 versions also as you might notice i tried to install it on a non-spaces path PostgreSQLx86.

Still no luck.

Anyone can help me out ?


texts: Changing font of legend for gxp OpenGeo?

texts:
I already changed fonts in all *.css files but a font of content of legend is still default. I think that I made changes in all style files.

How can I change a font of content of legend?


texts: How to integrate a Python Script into Processing Modeler in QGIS

texts:
I have a model in QGIS and I've generated a Python Script to calculate the matching between two strings (Levenshtein Distance with the FuzzyWuzzy Package) for one of the OUTPUTS.

##############################################################################
from fuzzywuzzy import fuzz

from fuzzywuzzy import process

from PyQt5.QtCore import QVariant

fn = 'W:\BENUTZER\OUTPUT_SP_2019_ToCheck.shp'

vLayer = QgsVectorLayer(fn, '', 'ogr')

vLayer.dataProvider().addAttributes([QgsField('Matching', QVariant.Int)])

vLayer.updateFields()

vLayer.startEditing()

for f in vLayer.getFeatures():

    f['Matching'] = fuzz.ratio(f['TITLE'], f['pName'])

    vLayer.updateFeature(f)

vLayer.commitChanges()

##############################################################################
Now with this calculated distance (stored in the field "Matching"), I have to continue with the workflow generation. My question is, there's a way to include this Python code into my Model, or how can I generate the value in the Graphical Modeleler obtained with the Python script?
This is the result that I obtained with the script

Thanks for your help!


texts: Convert Polygon to MultiPolygon using features attributes

texts:
I'm trying to convert SHP files to GeoJson using ogr2ogr and that works well with Polygon.

Now I'm struggling to convert multiple Polygon geometry features to a single one (MultiPolygon) for each CODE_COM.
Actually, the file looks like this:

{
"type": "FeatureCollection",
"name": "FEUILLE",
"features": [
{ "type": "Feature", "properties": { "FEUILLE": 1, "SECTION": "0A", "CODE_DEP": "01", "NOM_COM": "L'Abergement-Clémenciat", "CODE_COM": "001", "COM_ABS": "000", "ECHELLE": "2500", "EDITION": "2019-04-01", "CODE_ARR": "000" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 4.9307031, 46.173755 ], [ 4.9308624, 46.1737449 ], [ 4.9307031, 46.173755 ] ] ] } },
{ "type": "Feature", "properties": { "FEUILLE": 1, "SECTION": "0A", "CODE_DEP": "01", "NOM_COM": "Andert-et-Condon", "CODE_COM": "009", "COM_ABS": "000", "ECHELLE": "1250", "EDITION": "2019-04-01", "CODE_ARR": "000" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 5.6565458, 45.7883894 ], [ 5.6565725, 45.7883458 ],
...

I guess I need to join the data using an SQL query, the problem is that I don't really know how to do that.
According to the output, it seems obvious to me that the features should be grouped by the "CODE_COM" property, that helps referencing the data.

Here is the ogr2ogr command I'm using to get this result:

ogr2ogr -f "GeoJSON" -lco RFC7946=YES -t_srs "EPSG:4326"

I've tried adding an SQL query like that but something is missing:

ogr2ogr -f "GeoJSON" -lco RFC7946=YES -t_srs "EPSG:4326" -dialect sqlite -sql "select st_union(geometry) as geometry from FEUILLE"

EDIT:
I finally managed to make it work changing the SQL query like this:

select st_union(geometry), FEUILLE, SECTION, CODE_DEP, NOM_COM, CODE_COM, COM_ABS, ECHELLE, EDITION, CODE_ARR from FEUILLE GROUP BY CODE_COM

But I now have a big issue: the Polygons aren't converted as MultiPolygons features, so there are a lot of missing lines.

Here is the result without the SQL query (what's expected but with features grouped as MultiPolygons):

And here is the same result with the SQL query:


texts: Accessing system tables using ArcObjects?

texts:
I have a case to get the maximum of year value from the table schema in GeoDatabase.

More detail: I have a table named SUPPLYMASTER it contains the columns like Supply2000, Supply2001, Supply2002, Supply2003, ... up to Supply2015 and it can be extended in future.

I can write the logic in simple SQL with the system tables as:

SELECT MAX(REPLACE([name], 'Supply', '')) AS [MaxYear]
FROM sys.columns
WHERE object_id = object_id('SUPPLYMASTER') AND [name] LIKE 'Supply%'

It filters the column names starting with Supply from the SUPPLYMASTER table, the replace will remove the Supply text from the filtered column names and Max will return the maximum number from the filtered values. So I will get the result as 2015.

I tried the same in ArcObjects, but when trying to access the system table, it throws error.

I tried to access the system table like the below:

Config settings as:

<appSettings>
    <add key="Database" value="DatabaseName"/>
    <add key="Server" value="ServerName"/>
    <add key="UserName" value="UserName"/>
    <add key="Password" value="Password"/>
    <add key="Version" value="DBO.DEFAULT"/>
</appSettings>

OpenTable code as:

IPropertySet connectionProperties = new PropertySetClass();
propertySet.SetProperty("Database", System.Configuration.ConfigurationManager.AppSettings["Database"]);
propertySet.SetProperty("INSTANCE", "sde:sqlserver:" + System.Configuration.ConfigurationManager.AppSettings["Server"]);
propertySet.SetProperty("USER", System.Configuration.ConfigurationManager.AppSettings["UserName"]);
propertySet.SetProperty("PASSWORD", System.Configuration.ConfigurationManager.AppSettings["Password"]);
propertySet.SetProperty("VERSION", System.Configuration.ConfigurationManager.AppSettings["Version"]);

Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SdeWorkspaceFactory");
IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);
IWorkspace workspace = workspaceFactory.Open(connectionProperties, 0);

IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspace;

ITable sysTable = featureWorkspace.OpenTable("sys.columns");

IQueryFilter sysTableFilter = new QueryFilterClass();
sysTableFilter.SubFields = "name";
sysTableFilter.WhereClause = "object_id = object_id('SUPPLYMASTER') AND [name] LIKE 'Supply%'";

ICursor cursor = sysTable.Search(sysTableFilter, true);
IRow row = null;

while ((row = cursor.NextRow()) != null)
{
  // .... data processing
}

It throws the error in the line ITable sysTable = featureWorkspace.OpenTable("sys.columns");

  DBMS table not found

How can I access the system tables using ArcObjects?

Note: I tried with INFORMATION_SCHEMA.COLUMNS also, it throws the same error.


texts: QGIS Convert bathymetry shapefile from fillable to strokeable

texts:
In Java I would like to stroke the bathymetry depth lines instead of filling the area. I'm starting with the Natural Earth Bathymetric Data:  http://www.naturalearthdata.com/downloads/10m-physical-vectors/. When I apply a stroke with these files I get extra lines where the data is split into features.

Working in QGIS with the ne_10m_bathymetry_K_200 file I've tried dissolve, Polygons to lines, and in Grass tried to follow the steps outlined in this post:
Dissolve polygons in QGIS (GEOS) with tolerance

No matter what I do I'm left with a stubborn line just North of Australia and a line at the right edge of the screen.  If I zoom in on the line North of Australia it does look like there is a gap between these polygons.

How can I remove those last two lines?


texts: How to convert line simple features to topological network?

texts: using Gdal/ OGR is there a way to convert a shp, kml, or PostGIS line layer into a network of nodes and segments so it can be used in packages such as networkx?


texts: CRS for Lat/Long points and QuickMapServices

texts:
I have data from Google maps with latitudes and longitude recorded like this:

lat
long
location

35.05
135.76
Kyoto

44.24
-123.72
Corvallis

The coordinates are from Google maps. The only CRS which displays all of the points is WGS 84 EPSG: 3857; other CRS display the points in wrong locations, only the first point, or no point at all. When I use the EPSG: 3857, however, no QuickMapServices maps work. The most they will display is a blue background.


texts:
i've done the union first and then rasterized it. But now i want instead of that, to rasterize each one and then merge it as one. How do i do that? I mean  without the union first off.

 on<-union(union(extent(shp1),extent(shp2)),extent(shp3))
 r<-raster(on,res=100)
 r1<-rasterize(shp1,r,field="vat")
 r2<-rasterize(shp2,r,field="vat")
 r3<-rasterize(shp3,r,field="vat")
 raMerge <- merge(r3, r2, r1, fun=sum, na.rm=T)

i turned it to :

ra1<-raster(shp1,res=100)
ra2<-raster(shp2,res=100)
ra3<-raster(shp3,res=100)
r1<-rasterize(shp1,ra1,field="vat")
r2<-rasterize(shp2,ra2,field="vat")
r3<-rasterize(shp3,ra3,field="vat")
merge<-r1&r2&r3

but says

 Error in r1 & r2 :
 Cannot compare Rasters that have different BasicRaster attributes. See compare()


texts: Iterating through tables in map document and replacing workspace paths using ArcPy?

texts:
I am using Python 2.7 and ArcGIS 10.4.1 to create a script that walks through folders and scans for mxds. I use a set to get unique data sources and copy the data to a new folder keeping the same data format. I would then like to access the map document and repath all layers and tables to the copied data. I'm struggling to get the tables repathing code to work. In the below example, I have a fGDB table in the mxd and have copied the fGDB that it sits in to a new folder. I would thus like to repath to the table in the copied fGDB.

import arcpy, os
mxdrep = arcpy.mapping.MapDocument("CURRENT")
dfsrep = arcpy.mapping.ListDataFrames(mxdrep)

ignoreStringList = ["Vendor", "vendor", ".sde"]
#Target folder containing the copied fGDB
dbFolder = os.path.join(projFolderPath, "Databases")
for dfrep in dfsrep:
   tablesRep = arcpy.mapping.ListTableViews(mxdrep, "", dfrep)
   tabRepDataSrc = tableRep.dataSource
   #Excludes data sources containing list of sub strings
   if not any(folders in tabRepDataSrc for folders in ignoreStringList):
      #Create the table view for each table
      tableRepView = arcpy.mapping.TableView(tabRepDataSrc)
      # Get the directory level to replace
      tabPath = os.path.dirname(tabRepDataSrc)
      tabRepPath = os.path.dirname(tabPath)
      #Repath the tables
      tableRepView.findAndReplaceWorkspacePath(tabRepPath, dbFolder)

This gives a rather cryptic error;

Runtime error
Traceback (most recent call last):
  File "<string>", line 9, in <module>
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\_mapping.py", line 1284, in findAndReplaceWorkspacePath
    return convertArcObjectToPythonObject(self._arc_object.findAndReplaceWorkspacePath(*gp_fixargs((find_workspace_path, replace_workspace_path, validate), True)))
ValueError: Layer: Unexpected error


texts:
Trying to apply a map.SetFilter over points with continuous values, grouped by categories like:

0-8000 orange

8000-16000 green.

16000-30000 yellow.

30000-40000 blue.

Bigger than 40000 purple

As example, I want to visualize only values between 8000-16000 and >35000, that is, orange and purple

Jsfiddle: https://jsfiddle.net/perikut/x1cknef8/6/

Using 'any' gives me all colors. Using 'all' doesn't give me any color. Did not find any similar example around, maybe the formulation of the filter is not correct?

map.setFilter('hospitals',['any',[">", ["get", "npv_h"], 8000],["<", ["get", "npv_h"], 16000],[">", ["get", "npv_h"], 40000]])

Also, the first SetTimeout, with the filter below, only partially works. You can still see orange points and I checked that they have a value < 8000

map.setFilter('hospitals',['any',[">", ["get", "npv_h"], 8000]])

thanks in advance


texts: Is there a way to improve the render quality of ESRI basemaps when exporting a map to PDF or JPG?

texts:
I've had this problem with ESRI base maps before but not this bad. The problem is that when I use an ESRI basemap as my backdrop, it looks fine in my MXD but when I export the map to PDF, JPG, or whatever format, the quality of the basemap in the export file is a lot worse than what I see on-screen in ArcGIS. Here's the most recent, egregious example:

What I see in ArcMap:

What it looks like in my exported PDF:

I've tried several things to solve this problem, and the result you see above is the best I can do:

Turning up the export quality on the PDF to 600dpi or higher
Exporting to JPG, or to TIF with no compression at all
Clearing my cache on the basemap layers I'm using
Printing to the Adobe PDF print driver rather than using Export to PDF

Some people apparently have this problem due to a slow internet connection, but since I can see the tiles in my MXD and they look fine, I don't think that's my problem.

Is there something I'm missing?


texts: Overlapping polygons with buffers and need to keep associated data

texts:
I have a point file that has been buffered by 400m. I have overlayed it on top of a polygon file and I need to keep all data that overlaps for each point. What is the best way to do this?
I have unioned the 400m buffer (pink in the image) and then the blue selected parts are what I would need to keep for the top two points - but not the bits that overlap for each of the points...I need to do this for 1000's of records.


texts:
I'm new to GIS but I've already a signifiant experience in development. I work on new project. We're trying to use the power of Qgis to create a GIS tool for our commercial team.

QGis is istalled from OSGeo4W (3.8), Python (3.7) and PyQt (from pip) are installed too.

Correctly set PYTHONPATH :

set PYTHONPATH=C:\OSGeo4W64\apps\qgis\python

When i try to run my python script (copied from the official documentation "PyQGIS Developper Cookbook"), we've got an error :

from qgis.core import *

# Supply path to qgis install location
QgsApplication.setPrefixPath("C:\OSGeo4W64\apps\qgis", True)

# Create a reference to the QgsApplication.  Setting the
# second argument to False disables the GUI.
qgs = QgsApplication([], False)

# Load providers
qgs.initQgis()

# Write your code here to load some layers, use processing
# algorithms, etc.

# Finally, exitQgis() is called to remove the
# provider and layer registries from memory

qgs.exitQgis()

ImportError: DLL load failed: %1 n’est pas une application Win32 valide.


texts: I want to pass a feature service as an input(Input parameter type -GPFeatureRecordSetLayer) to a geoprocessing service. can anyone tell me a sample code? I read all the relevant questions available, but I didn't find an answer in which feature service is used as an input.


texts: Why does local ArcGIS API for JavaScript over https trigger info request over http?

texts:
I got the esri js api 3.14 deployed on a server and an application that is reachable over https.

Now I got a proxy script that delivers the map requests.

When you look at the requests of the page you can see that there is ONE request (an info request to the rest interace) that goes directly to the internalServer.

The Chrome Browser says:

Firefox says:

I added the proxy script rule like this:

urlUtils.addProxyRule({
                urlPrefix: "http://srvutdemo102/arcgis/rest/services/utjsc_de/Grundkarte",
                proxyUrl: "ESRI_proxy.jsp"
            });

Is it possible to tell the js-api to use the proxy script for this "initial" info request?


texts: Calculating number of days between first and last frost in Google Earth Engine?

texts:
I want to do two things.

First, I want to map the number of nights that have temperatures below 0 in one year. This code seems to do this.

Second, I want to calculate the number of days between the first and last nights with sub-zero temps. I managed to get the first one in the series but the output is a format I don't know: 994118400000.

How can I get the last night the series and measure the span of time between the first and last dates in the series?

var NightSA = ee.ImageCollection('MODIS/006/MOD11A1').select('LST_Night_1km');

var count2001 = NightSA.filterDate('2001-07-01','2002-06-30').map(function(NightSA)
{return NightSA.lt(0).set('system:time_start',
NightSA.get('system:time_start'))});

var first2001 = count2001.first().get('system:time_start');
print (first2001);

var nightSurfaceTempVis = {
  min: 0,
  max: 365,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};

Map.addLayer(count2001.count(),nightSurfaceTempVis,'Nights < 0');


texts: Why does the ArcMap exported KMZ always show my local path information when displayed in Google Maps API?

texts:
I used ArcMap (arcpy) to export shapefiles into KMZ and then display them via Google Maps API v3. It looks fine initally:

But when I hover my mouse above the KMZ and click it, it will alway show my local path information:

It is very annoying. Is there a solution to avoid this? I just want the KMZ to show the shape...


texts: Unable to open shapefile with GRASS GIS?

texts:
Is GRASS GIS able to work with a shapefile (*.shp)?

I can´t open my shapefile (328000Kb) in that program.

I want to dissolve a shapefile in GRASS because on QGIS the process is very slow.


texts: Extracting tree cover values from Hansen to feature collection using Google Earth Engine?

texts:
I would like to extract mean pixel value of treecover 2000 from Hansen from each geometry (polygon) of different shapefiles that I have import into Google Earth Engine (var table_comp). I create for this reason a  for loop.

This is briefly my script:

var gfc =ee.Image('UMD/hansen/global_forest_change_2018_v1_6');
var treecover =gfc.select(['treecover2000'])

for (var py=0; py<19; py++){
  //var py=9;
  var year=py+2000

  var table_comp = load feature collection

  var tree_mean_value = treecover.reduceRegions({
    collection: table_comp,
    reducer: ee.Reducer.mean(),
    scale:30});

I would like to ask some clarification about the option scale of reduceRegions, in the help I read "A nominal scale in meters of the projection to work in" but I believe to not understand very well. I put 30 because I think this as spatial resolution of my starting image.

Can anybody can help me?


texts: Updating layer datasource in multiple map documents using ArcPy gives SyntaxError?

texts:
I would like to update the data source of a specific layer in multiple map documents using ArcPy.

However, I'm getting a <type 'exceptions.SyntaxError'>: invalid syntax error on if lyr.datasetName == inputlayer in the code below.

import arcpy, string, os
import arcpy.mapping as MAP

try:
    MXDList = string.split(arcpy.GetParameterAsText(0), ";") #ArcMap document(s)
    inputlayer = arcpy.GetParameterAsText(1) # featureclass
    newwork = arcpy.GetParameterAsText(2) # Workspace

    for MXDPath in MXDList:
        MXD = MAP.MapDocument(MXDPath)
        for dataframe in arcpy.mapping.ListDataFrames(MXD):
            lyrs = arcpy.mapping.ListLayers(dataframe)
            for lyr in lyrs:
                if lyr.datasetName == inputlayer
                    lyr.replaceDataSource(newwork, "FILEGDB_WORKSPACE", "Center")
    MXD.save()

What causes the SyntaxError here?

I was attempting to reference this forum thread.


texts: how to reproject an image programatically in c# using gdal

texts:
I have to reproject an image programmatically in c# using gdal library.
For the moment i'm a little bit lost in gdal and i don't really now which methods i have to use.
I think i will have to calculate tie points, set reference systems and then reproject but i don't now how to manage this.

Can somebody help me with exemples or tutorials?

Best regards,
Christophe.


texts:
After a working break of a couple of weeks, I start do work on my last project. But I can not edit my postgis-layers anymore! Its possible that I did several system-upgrades, and maybe there was also a postgis-update, and I updated to QGIS 2.12 but I can not exactly tell, what could be the reason.

I have the right ownership of the table and all Privileges (arwdDxt*)

As you can see in the editing-bar - I can not make any changes to the layer


texts: The TIN interpolation tool does not work within the processing models in Qgis 3.2

texts:
I want to use the TIN interpolation tool within a model that I am creating in the process toolbox; and so far I have not been able to use it. it does not let me connect the vector layers with the TIN tool. What should this problem be?

I try to connect the layer and can not ...

It does not allow to make the connection between the input parameters of the model with the TIN interpolation algorithm, because I can not select the layers to interpolate ...


texts: Calculate NDWI of each image in a Imagecollection?

texts:
var L8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
.filterDate('2016-06-01', '2016-07-31') // you can change date here
.filter(ee.Filter.lt("CLOUD_COVER", 0.75))
.filterBounds(collection[i])
.map(maskLandsatclouds)
.sort('DATE_ACQUIRED')
.select(L8bands)

 var count = L8.size();
 for (var j=0;j< 10 ;j++){
 // Get the number of images.
 print(j);

Correct if wrong; the code below is taking the median of all available images in a collection.

 var mosaic_L8 =L8.median().clip(collection[i]); // here we are taking the
 median at each pixel in the collection
 Map.addLayer(mosaic_L8, L8vis, "mosaic_L8")

...and here, it is calculating NDWI of an image:

 // Create an NDWI image, define visualization parameters and display.
 var ndwi = mosaic_L8.normalizedDifference(['B3', 'B5']);

 var ndwiViz = {min: 0.5, max: 1, palette: ['44c9f1', '1637f1']};
 //Map.addLayer(ndwi, ndwiViz, 'NDWI', false);

 // Mask the non-watery parts of the image, where NDWI < 0.10.
 var ndwiMasked = ndwi.updateMask(ndwi.gte(0.10));
 Map.addLayer(ndwiMasked, ndwiViz,'NDWI masked '+ ee.Number(j).getInfo());

 Map.centerObject(ndwiMasked);

I have enclosed this code in a for-loop to perform the operation of NDWI on each image in a collection, but am unable to access individual images. It should be something like L8.normalizedDifference(['B3', 'B5']);
but this is not a correct command.


texts: pyqgis create layer from a list of points

texts:
I am working with pyqgis in QGIS3 and i want to create a point layer from an intersection beetween a point and a polygon layer. I can identify wich points intersects my polygon layer but i can't use them in a new layer. I think i have to store the result in a new point layer with a list, but it not seems to be the good argument.

How can i "populate" a point layer with a list like ?

[<QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>, <QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>, <QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>, <QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>]

Here is my code

# Points
myPoints = "C:\pyqgis\mypoints.shp"
layerPoints = QgsVectorLayer(myPoints, '', 'ogr')

# Polygon
myPolygon = "C:\pyqgis\mypolygon.shp"
layerPolygon = QgsVectorLayer(myPolygon, '', 'ogr')

# Intersection
featureList = []
for i in layerPoints.getFeatures():
        for p in layerPolygon.getFeatures():
            if i.geometry().intersects(p.geometry()):
                featureList.append(geom.asWkt())

# Create layer with result
intersection = QgsFeature()
intersection.setGeometry(QgsGeometry.fromPointXY(featureList))

# Add it to the map
result = iface.addVectorLayer(intersection, "result", "ogr")


texts: Batch export map series pages to GIF files using ArcPy in ArcGIS Pro

texts:
I have created a map series on ArcGIS Pro 2.5.1 displaying distribution data for different bird species over a UK county. Each page displays a different species. The index feature is a point feature class which is just a list of the different species, each with the same coordinates to centre the map. I want to export each page as a GIF file.
This Python code snippet from this help page seemed to be for the same general task:
https://pro.arcgis.com/en/pro-app/arcpy/mapping/mapseries-class.htm
relpath = os.path.dirname(sys.argv[0])

p = arcpy.mp.ArcGISProject(relpath + "\\MapSeries\\US_States.aprx")
l = p.listLayouts()[0]
if not l.mapSeries is None:
  ms = l.mapSeries
  if ms.enabled:
    for pageNum in range(1, ms.pageCount + 1):
      ms.currentPageNumber = pageNum
      print("Exporting {0}".format(ms.pageRow.STATE_NAME))
      pageName = ms.pageRow.STATE_NAME
      l.exportToPNG(os.path.join(relpath, "Output", f"Ex2_{ms.pageRow.STATE_NAME}.png"))

I have tried to modify this to fit my own workflow:
relpath = os.path.dirname("C:\Yusef\Pembrokeshire Avifauna\Breeding 2000s gif")
p = arcpy.mp.ArcGISProject("C:\Yusef\Pembrokeshire Avifauna\Pembrokeshire Avifauna_working.aprx")
l = p.listLayouts()[0]
if not l.mapSeries is None:
  ms = l.mapSeries
  if ms.enabled:
    for pageNum in range(1, ms.pageCount + 1):
      ms.currentPageNumber = pageNum
      print("Exporting {0}".format(ms.pageRow.Species))
      pageName = ms.pageRow.Species
      l.exportToGIF(os.path.join(relpath, "Output", f"Ex2_{ms.pageRow.Species}.gif"))

For which I get the following error message, which I'm not sure how to interpet:
Traceback (most recent call last):
  File "<string>", line 12, in <module>
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\_mp.py", line 1029, in exportToGIF
    return convertArcObjectToPythonObject(self._arc_object.exportToGIF(*gp_fixargs((out_gif, resolution, gif_color_mode, clip_to_elements), True)))
OSError: C:\Yusef\Pembrokeshire Avifauna\Output\Ex2_Barn Owl.gif

Can anyone tell me what I'm doing wrong?


texts: Clip and ship by a queried layer extent and by user defined extent

texts: I am trying add a functionality to the clip and ship script. I would also like the users to extract data by the extent of a queried feature layer. For instance, if County = "Alameda" then extract data that falls only within the Alameda county extent. Else, allow the users to hand draw an selection area. If any of have clues or ideas I would love to hear about it.


texts: Field survey and Google Earth crop points

texts:
I am trying to apply machine learning to classify and distinguish crop types in my AOI, using optical satellite images. I don't have access to a pre-classified land use land cover reference map. And since this is a supervised classification, ground truth data is mandatory for a successful and accurate result.
So the next best solution is to organize a field campaign to collect data (GPS data, crop types, sewing dates, current state: growth stage of the plant...). I am not familiar with such a task, hence I wanted to get some insights on how to proceed since this is not an easy task and it is costly.
I already checked the previously asked questions here with similar tags, but I haven't found anything related or answers my question which is how to organize and prepare the field campaign. What I think is that I need first to visualize my AOI on Google Earth, check the land parcels with a good level of "greeness" and road accessibility, then check NDVI spectral signatures to somewhat infer the number of crops existing there and finally, choose locations (mark points) based on NDVI ensuring a somewhat equal number of locations for each possible crop type.


texts: What events are fired when a feature is copy and pasted?

texts:
I am trying to detect when a feature has been copied and pasted.  The OnCreateFeature event fires, but I cannot tell if the feature object was copy/pasted.  Is there an event or other method to differentiate pasted features from newly created features?

EDIT: Two possible approaches I've thought of, but not sure of efficacy or downsides.

In the OnCreateFeature event, I can read the 'new' feature's attributes, and compare them to the DefaultValue of each field.  If every value is default, it's probably a new feature and not a copied feature.
Get the geometry of the 'new' feature, and compare it to the geometry of every selected feature.  If there is a match, the feature is a copy/paste.


texts: Looking for all points in a square around a given point in PostGIS

texts:
For a given point (x,y), I would like to find all the points in the neighborhood (defined by myRadius, a number, in meters).

For now, I'm using a circle around each point, but some reason, I'd like to draw a square around each point instead of a circle (square in the orientation of x-y coordinates).

Command I use for the circle (distance):

SELECT * FROM table_of_points WHERE (ST_Transform(the_geom,21781) && ST_expand(ST_SetSRID(ST_GeomFromText('Point(x y)'), 21781), myRadius )

For the square, can I use:

SELECT * FROM table_of_points WHERE (ST_Transform(the_geom,21781) && ST_MakeEnvelope(x-myRadius, y-myRadius, x+myRadius, y+myRadius, 21781 ) )

It works, but I'm just wondering if it's the fastest way of doing this.


texts: Why does multi-threaded rendering not use all the CPU potential in QGIS 2.4?

texts: I have just downloaded QGIS 2.4 which is great and ran some tests. I have enabled the multi-threaded rendering but it seems that the program do not use all the potential of my CPU. Attached pic shows CPU usage when QGIS is rendering (top) and some other multi-threaded program I am using. As You can see the "other" program takes all the cores to a 100% while QGIS only took half of them. The rendering time in this case was 78586 ms, so it's not like it couldn't be any faster. Max cores to use is set to 8. Is this normal or am I missing something?


texts: Moving labels after converting to annotation in ArcGIS for Desktop?

texts:
I'm trying to move some of my labels individually around some of these points so that they don't overlap with each other or crowd each other out. When I converted them, I chose to store the annotations in the map and used a reference scale of 1:100:000.

But for some reason I can't seem to move them after I ran the conversion. When I use the select elements tool, I can't select the labels. It just drags a box over the map like normal; I can't pick up the labels and move them. I've attached a picture to describe what I mean (I've got Gas_export selected, but Standalone ramps is the layer I converted to annotation).

I've got the data frame activated (there's only one), I'm not editing, and I checked to make sure I'm still using the same reference scale. I'm sure it's something totally obvious that I've missed, but I'm really stuck.


texts: Trying to turn a GRID Stack 7.x to GeoTIFF

texts:
I have a set of imagery that contains the following folders:

aerial
aerialc1
aerialc2
aerialc3
info

When I load the folder (or folders) into ArcMap, they all display perfectly. When I try and export to GeoTIFF, the resulting imagery is always black and unusable. I've tried Project Raster as well since I wanted to change the projection type. Nothings seems to be working. I also attempted to load the folders into Global Mapper, but that didn't work either.

I just need to get the imagery into a GeoTIFF format. Its about 770 GB in total.


texts: Loading KML file with ESRI javascript API

texts:
I am trying to load a kml file using the ESRI javascript API v. 3.5.  What I did with v. 2.6 no longer works.  How do you add KML data with the latest esri javascript api?  I am testing on my local machine, which used to work fine.

        <link rel="stylesheet" href="http://serverapi.arcgisonline.com/jsapi/arcgis/3.5/js/dojo/dijit/themes/claro/claro.css">
        <link rel="stylesheet" href="http://serverapi.arcgisonline.com/jsapi/arcgis/3.5/js/esri/css/esri.css">

In the console I get an error

  Origin null is not allowed by Access-Control-Allow-Origin.

This seemed to work fine with the older api why the difference?  How can we test locally since kml files have to be a public url without running a server?

EDIT

I ended up installing IIS on Windows 8.  I then used localhost.  In Chrome it would not work still gave me an error of

  Origin ...localhost is not allowed by Access-Control-Allow-Origin

So I tried IE again and it worked with localhost once I enabled intranet settings.  However, I did discover if you use the computer name ie.. http:// computername instead of localhost it will work in chrome as well.


texts: where i can add style in this code leaflet

texts:
var urlborne = 'http://localhost:8080/geoserver/Gestion-proprete/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Gestion-proprete:borne&maxFeatures=2000&outputFormat=text/javascript&format_options=callback:getJson&srsName=EPSG:4326';

    var Borne = new L.GeoJSON();

$.ajax({
    jsonp: false,
    url: urlborne,
    type: 'GET',
    dataType: 'jsonp',
    jsonpCallback: 'getJson',
    success: function (data) {
        console.log(data);
        Borne.addData(data);
        map.fitBounds(Borne.getBounds());
    }
});


texts:
I've followed instructions in tileserver-gl github page and successfully install tileservre-gl, I mean npm command finished without error and all packages were downloaded.

zurich_switzerland file is also downloaded.

But when according to instructions I run tileserver-gl zurich_switzerland.mbtiles it says:

command not found.

What's the problem?


texts:
I need to extract together the values of a raster tiff (radar SAR information) with a layer of several polygons (each polygon is a class, resulting from a classification of 5 soil types).
I want to know all the values contained in each class to export the attribute table.
I have tried the option to extract by polygon but it is only for a single polygon . I also tried the option to extract by layer of points, but it gives me an error because the raster has a lot of information. Also by mask...
Any other ideas?


texts: Automatically changing field entry in related table when entering it in feature class?

texts:
I am working to try and develop a code or process that I can run that will allow me to automatically change an XY field in a related table when I change an associated XY field. Let me explain our situation a little first. We have a system where we are trying to maintain utility data in an SDE, and while doing this data is being entered by technicians in the field and brought back in. When data is entered in the field, there is a FACILITY_ID associated with each item, such as hydrants, values etc. This data is pushed to each cooresponding our feature class. However, when multiple actions are taken on the same feature, we are storing each of those changes in a different table while only reflecting the last action in the feature class. These two have a relationship class. So, there is only one entry for each FACILITY_ID in a feature class, but there may be multiple entries in the 'related' table for each FACILITY_ID.

Our problem is what to do if we move a hydrant/value in ArcMap and have to change the XY coordinate. We don't mind recalculating the X and Y values in the feature class table, but we would like a process that could recognize that change in the feature class XY and automatically change all of the X and Y values for the associated related table FACILITY_ID without having to go in and edit them by hand. Is there a simple way to do this? We've been looking at ArcObjects or ArcPy, but does a command like this already exist in ArcMap?

I should clarify, the 'related table' is more just a table that has a simple relationship class connected to the feature class. The second table in a non-spatial one.


texts: Percent cover from overlapping shapefiles

texts:
Apologies if this has been answered elsewhere. I have tried for several days now to come up with a solution to this problem, but haven't found an exact answer that applies to my current situation (or at least one that I could understand enough to implement myself).

I am trying to extract the percent of particular drainage areas that are covered by different surficial geologies (i.e. sand, silt, bedrock, etc.). I have shape files for almost 100 catchment areas (some overlap each other due to being sub-watersheds of larger river systems). I also have a shape file of surficial geology.

Below are some images to help illustrate:
The first is what the catchment areas look like. The second is the geology layer with 3 catchments displayed over top (note the nestedness).Finally, I am looking to produce a table of percentages with catchments against geological features.

Right now I am only using QGIS, but am open to trying other programs as well. I have used R in the past, but am not very good at python.

EDIT: the result of using union on the geology and catchment layers:


texts:
I am new to QGIS and I want to overlay an old hand drawn forest map (.png file) over a British OS map (.tif file) of the same forest, so I can digitize the compartments in the old forest map.

I have dragged both files into a new project, and they are in the window, but their respective positions are incredibly far away from each other, they aren't "on top" of one another. I have to click "zoom to layer" to find each layer.

I searched through this website and google, and the solution which always comes up is to make sure they have the same coordinate system. If I right click on the properties of both my layers, they both say "EPSG:7405 - OSGB 1936 / British National Grid - Projected" in the CRS section - so I don't know what else to do.  Can anyone help? Thanks in advance.


texts: Looping ImageCollection and FeatureCollection together in Google Earth Engine

texts:
I am looking for a way to filter images for certain dates per feature in the feature collection. I want a function that would take each individual feature from the feature collection and then check from the image collection if there is an S1 image available for that region for a particular date. The dates are also listed in a column in the feature collection for each feature. So for feature1 is there an image for 2016-04-17 in the S1 collection, if yes select that feature. The feature collection table contains a column called 'SrcImgDate'. So I just want to filter out a feature collection with features that have an image for the srcimgdate next to them.

Expected output: FeatureCollection with all features that intersect with an image at the given srcimgdate and location

This is the imageCollection with a general date range

var Sentinel1 = S1.filter(ee.Filter.eq('transmitterReceiverPolarisation', pol))
                  .filterMetadata('instrumentMode', 'equals', 'IW')
                  .filterDate('2016-04-01','2016-08-30' )
                  .select(['VH','VV'])
                  .filterMetadata('resolution_meters', 'equals' , 10);

Here is the link to the FeatureCollection. You can limit it to 1000, as it has 11000 feature and would give an error to use such a large dataset.
https://code.earthengine.google.com/?asset=users/rawailnaeem/CA


texts: Geoprocessing Service to return Multidimensional Raster from a NetCDF file

texts:
I have tried to make this work for weeks with no results. I am trying to create a model in ArcMap that uses the Multidimension Tool "Make NetCDF Raster Layer" and return a raster to which I can apply a Time Slider with the Silverlight API once its result is published as a Geoprocessing Service (ArcGis 10.1).
If I create such model and run it in ArcMap, and then publish the result, then I can successfully call the service as an ArcGISDynamicMapServiceLayer and link a Time Slider to it on my web app. But these type of NetCDF files contain more than one variable and for obvious reasons I want to be able to select any of such variables from my web app and get a different raster (with its time dimension) for each one of them.
So it occurred to me that I could publish a model that gets the variable by input parameter and automatically generates the raster I need. But when I run this Geoprocessing Service and assign the resulting raster to an ArcGISDynamicMapServiceLayer within the Silverlight environment I get a Raster that I can display but it doesn't have a time dimension/extent.

Any suggestions?


texts: Seeking Easy Guide to ESA satellite data?

texts:
Can anybody provide user friendly guidance to ESA datasets?

What is available?

When will new stuff be available?

I struggle  to work out what might be interesting amongst the many websites and their somewhat verbose format. For example this question indicates that there is probably some high resolution imagery available or coming soon somewhere.... When I try to find anything I find page after page after page about data access but very little information about the type of data itself.

The whole Copernicus Project looks fascinating!

I understand that it is all very new and many things aren't functioning yet nevertheless I would love to begin preparing ideas.

Can anybody point me in the right direction?


texts: error message when clipping raster with vector layer

texts:
When I try to crop a raster using QGIS Clipper tool, I get the following error message: ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.. Why?

The code QGIS generated to perform the clip was as follows: $ gdalwarp -q -cutline <clip_mask>.shp -crop_to_cutline -of GTiff <input_raster>.tif <output_raster>.tif

The same error message is generated if I try to run the gdalwarp command by itself from the commandline.


texts:
I want to spatial join two tables with the same polygons, so I only want to extract where the polygons does not merge with the other.

SELECT
       e.SGid as 'table_name2'
From [table]
    inner join table2 e with (index = s123_idx, nolock) ON e.SHAPE.STWithin(j.Shape.STCentroid()) = 1
         where (SUBSTRING(j.SG_ID,1,8) + SUBSTRING(j.SG_ID,12,13)) <> e.SGid


texts: What is the best way to handle CRSs in Geopandas so that both versions =0.7 are handled?

texts:
I am developing a python package called RivGraph. I've recently run into an issue related to geopandas' decision to overhaul CRS representation, as described here and here. The short of it is that in geopandas versions <0.7, crs's were stored as dictionaries of proj4 strings, e.g. {'init':'epsg:4326'}. In geopandas 0.7 and moving forward, crs's are stored as pyproj objects, which are richer and more robust. Ok, great.

I share RivGraph as a conda-installable, and direct people to install it using conda to ensure all the dependencies are worked out. What I've discovered is that, depending on the user's platform, sometimes geopandas <0.7 is installed, and sometimes geopandas 0.7 is installed. So that leads me to ask this question:

What is the best way to make my package compatible with both of these versions of geopandas? I guess I realize that I can check the geopandas version in my package, and write if statements that use one syntax vs. the other, but this seems tedious and ugly. I can't quite figure out if there is a single syntax that is compatible across both versions of geopandas when it comes to crs management.


texts: Creating summary text layer like realestate.com.au using CartoDB?

texts:
I've been playing around with the cartodb platform and haven't been able to find out how to create a data summary or view of the averages for given location. I see that the team at realestate.com.au managed to do it. But I'm not sure how.

Any ideas?