texts: Query for WAY that connects two INTERSECTIONS

texts:
I'm trying the retrieve the "WAY" that passes two given "INTERSECTIONS".

What I do is query for the NODE that intersects the INTERSECTION using the locations of this intersections (which is returned by OSRM). I use the locations as a bounding box in overpass:
    nodes(north, west,south,east);
    way(bn);
    out;

The result is there, but all ways attached to the two nodes are returned. How do I filter so only the WAY that passes both NODEs is returned?


texts:
When I am converting a raster using the polygonize tool, I get this error:

ERROR 4: 'C:/Users/Vitor/Desktop/Cartografia EE_Portimo/COS RASTER/cos_raster_reeditada.tif' does not exist in the file system,
and is not recognised as a supported dataset name.

I tried with other rasters and keep getting the same error.

Does anyone know what it is about or how I can solve it?


texts: Importing Python Library Inside QGIS Plugin for all OS using Subprocess

texts:
I am writing a QGIS python plugin which needs pandas library to function. I have searched the GIS.SE and there are many similar questions but none talks about using subprocess and that if it will work on all OS. People have suggested using pip.main() but the pip.main() has been deprecated. I have developed the following method of importing/installing pandas.
I need to know if this will work on all OS. I am on Windows and this works fine on windows.
try:
    try:
        import pandas as pd
    except ImportError:
        import pathlib as pl
        import subprocess
        import sys

        qgis_Path = pl.Path(sys.executable)
        qgis_python_path = (qgis_Path.parent / "python3.exe").as_posix()

        subprocess.check_call(
            [qgis_python_path, "-m", "pip", "install", "--user", "pandas"]
        )
        import pandas as pd

        feedback.pushInfo(
            "python library pandas was installed for qgis python"
        )
except:
    feedback.reportError(
        "Failed to import pandas. Tried installing pandas but failed.\nPlease manually install pandas for the python that comes with your QGIS or run this algorithm with Skip Output option for Report as HTML."
)


texts: How to get users click a map (OSM) and collect location data plus parameters?

texts:
I wonder if there is a simple service out there on the web that allows to collect location data from users easily. A click on a OSM map should give the location, then a form to collect a comment, value, category etc.
Think http://www.fixmystreet.org on a much simpler level, Google Forms with added map, kind of.

Idea: Collect incidents where cyclists experience dangerous situations, generate heat map where most threats are.


texts: Blank Popup when clicking on the map Leaflet-WMS

texts:
I created an interactive map using Leaflet, when I click on my layers I can get data in popup, my problem is that when I click on the map where there are no features it shows a blank popup. I want to block this popup when it's blank.

This my code:

map.addEventListener('click', onMapClick);
popup = new L.Popup({maxWidth: 1000});
function onMapClick(e) {
var latlngStr = '(' + e.latlng.lat.toFixed(3) + ', ' +         e.latlng.lng.toFixed(3) + ')';
var BBOX =         map.getBounds()._southWest.lng+","+map.getBounds()._southWest.lat+","+map.getBounds()._northEast.lng+","
+map.getBounds()._northEast.lat;
var WIDTH= map.getSize().x;
var HEIGHT = map.getSize().y;
var X = map.layerPointToContainerPoint(e.layerPoint).x;
var Y = map.layerPointToContainerPoint(e.layerPoint).y;
var layerObj = {};
layerObj["test_layer:EV_Lot_ALHOUDA_wgs"] = alhouda6_layer;
layerObj["test_layer:Imm_Coll_Lot_ALHOUDA_wgs84"] =alhouda7_layer;
layerObj["test_layer:Villa_Lot_ALHOUDA_wgs84"] = alhouda8_layer;
layerObj["test_layer:Eco_Lot_ALHOUDA_wgs84"] =alhouda5_layer;
layerObj["test_layer:Equip_Lot_ALHOUDA_wgs84"] = alhouda4_layer;
layerObj["test_layer:Rec_Lot_ALHOUDA_wgs84"] =alhouda3_layer;
layerObj["test_layer:Parking_Lot_ALHOUDA_wgs84"] = alhouda2_layer;
layerObj["test_layer:P200Log_Lot_ALHOUDA_wgs84"] =alhouda1_layer;

layerObj["test_layer:R3_Lot_HFOUNTY_wgs84"] =founty1_layer;
layerObj["test_layer:R5_Lot_HFOUNTY_wgs84"] =founty2_layer;
layerObj["test_layer:R9_Lot_HFOUNTY"] =founty3_layer;
layerObj["test_layer:R15_Lot_HFOUNTY_wgs84"] =founty4_layer;
layerObj["test_layer:R1_Lot_HFOUNTY_wgs84"] =founty5_layer;
layerObj["test_layer:R2_Lot_HFOUNTY_wgs84"] =founty6_layer;
layerObj["test_layer:Equip_Lot_HFOUNTY_wgs84"] =founty7_layer;
 queryLayers = [];
for (key in layerObj) {
if(map.hasLayer(layerObj[key])) {
 queryLayers.push(key);
 }
}
var urlLayers = queryLayers.join();

var URL = 'http://localhost:8082/geoserver/test_layer/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&LAYERS=test_layer:Imm_Coll_Lot_ALHOUDA_wgs84,test_layer:Villa_Lot_ALHOUDA_wgs84,test_layer:Eco_Lot_ALHOUDA_wgs84,test_layer:Equip_Lot_ALHOUDA_wgs84,test_layer:Rec_Lot_ALHOUDA_wgs84,test_layer:EV_Lot_ALHOUDA_wgs,test_layer:P200Log_Lot_ALHOUDA_wgs84,test_layer:Parking_Lot_ALHOUDA_wgs84,R3_Lot_HFOUNTY_wgs84,R5_Lot_HFOUNTY_wgs84,R9_Lot_HFOUNTY,R15_Lot_HFOUNTY_wgs84,R1_Lot_HFOUNTY_wgs84,R2_Lot_HFOUNTY_wgs84,Equip_Lot_HFOUNTY_wgs84&QUERY_LAYERS='+urlLayers+'&propertyName=Zone_Urb,Type_Urb&STYLES=&BBOX='+BBOX+'&FEATURE_COUNT=5&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&FORMAT=image%2Fpng&INFO_FORMAT=text%2fhtml&SRS=EPSG%3A4326&X='+X+'&Y='+Y;
 popup.setLatLng(e.latlng);
popup.setContent("<iframe src='"+URL+"' width='400' height='100' frameborder='0'></iframe>");
if (urlLayers)

map.openPopup(popup);


texts: In OpenLayers 3, how do I draw an extent as a LinearRect?

texts:
In OpenLayers 2, this was an easy feat. I was able to simply instantiate a LinearRect given the coordinates of my BoundingBox, and then add the LinearRect as a feature to a vector layer I had added to the map. In OpenLayers 3, I can't seem to find any easy way to do this using extents. Say I have the current extent (it is transformed to match OSM, since the back-end for this application serves EPSG:4326):

# Colorado Springs
extent = ol.proj.transformExtent([aoi.west, aoi.south, aoi.east, aoi.north], "EPSG:4326", "EPSG:900913")

I am trying to add this extent to the map using the following code:

lr = new ol.geom.LinearRing([
  [aoi.west, aoi.north],
  [aoi.east, aoi.north],
  [aoi.west, aoi.south],
  [aoi.east, aoi.south]
]).transform("EPSG:4326", "EPSG:900913")

layer = new ol.layer.Vector
  source: new ol.source.Vector()
  style: new ol.style.Style
    fill: new ol.style.Stroke
      color: '#ffcc33'
      width: 2
    stroke: new ol.style.Stroke
      color: '#ffcc33'
      width: 2
    image: new ol.style.Circle
      radius: 7
      fill: new ol.style.Fill
        color: '#ffcc33'
  features: [
    new ol.Feature
      color: '#BADA55'
      where: 'inner'
      geometry: lr
  ]

map.addLayer(layer)

There are no errors thrown – the only problem is that the shape is not displayed on the map (although the base OSM layer is). Checking the maps layers (console.log map.getLayers().getArray()) reveals that the vector layer is making its way into the map.

For reference, the aoi object looks like this:

aoi = { west: -104.8, south: 39.8, east: -104.7, north: 39.9 }


texts: Markers display with leaflet and javascript

texts:
I have a problem displaying data relating to schools on a div 'info' when click the marker with leaflet and javascript:

here's the error:

Uncaught TypeError: layer.setStyle is not a function


texts: How are bound variables specified in an ArcGIS REST SOE?

texts:
In the documentation for using handler functions with ArcGIS Server Rest SOE's, they don't explain how the keys in the boundVariables collection are determined.  When I search through the entire solution for "customLayersID" I don't see anything.  How would I know to use "customLayersID" as a key?

private byte[] CustomLayer(NameValueCollection boundVariables, string outputFormat,
    string requestProperties, out string responseProperties)
{
    responseProperties = null;

    //The LayerID.
    int layerID = Convert.ToInt32(boundVariables["customLayersID"]);

    //Execute.
    CustomLayerInfo layerInfo = GetLayerInfo(layerID);

    string json = layerInfo.ToJsonObject().ToJson();

    return Encoding.UTF8.GetBytes(json);
}


texts: Is it possible to use a heat map as a selection mask of elements?

texts:
I have a layer with thousands of polygons representing buildings, as is normal, buildings are concentrated in large population centers, but my map will be used mostly in rural areas. So I do not need so many buildings in big cities, but It is essential that all buildings that are far from the cities appear.

What would be the best tool or method to achieve this?

I look for some type of tool or method of work that can eliminate "random" and in a spatially distributed way the small buildings, preferably those near / within the large population centers.

UPDATE-> I think that using a heatmap like "filter" to select features maybe will be a solution, but I don't know if is it possible.

NOTE: I convert polygon layer to centroids layer and I got a HeatMap.

You can get example shapefile here-> builds.shp


texts: Subsetting Sentinel 2 image in Google Earth Engine

texts:
I am new to Google Earth Engine and I am working with Sentinel-2 imagery for landcover classification.

I need to limit the Sentinel 2 images to a specific study area i.e. Area of Interest. What specific commands do I need to achieve this?


texts: ArcHydro: Flow Accumulation 863988034 Cells

texts:
I've been using Arc Hydro for the last 7 years to determine catchment characteristics. I'm currently setting up a new model to determine structure locations along a proposed railway alignment. This is the largest study area I've ever had to process and currently its failing when I try to generate the Flow Accumulation Grid. My Flow Direction Grid for the study area has 863 398 034 cells.

I'm looking for advice in how to successfully generate my Flow Accumulation Grid. If ArcGIS can't handle such large DEM's will GRASS be able to process such a large DEM and generate a ArcGIS GRID as an output that I could use for further processing of catchment characteristics.


texts: Find tree whose closest to the greatest variety of other trees using sf and R

texts:
Given a large spatial database of trees (600,000) find the tree that within radius r encircles the greatest variety of types of trees.

Looking at an answer by andyteucher I've simulated some tree data below using random trees with types A to I. I've then used buffer to draw a circle around each tree and intersect that with all the trees. Not sure what's next or if this direction is appropriate.

library(sf)
library(tmap)
library(tidyverse)
library(leaflet)
library(stringi)

mean_lat <- 43.9
sd_lat <- 0.1
mean_long <- -79.384293
sd_long <-  0.1

set.seed(42)
trees <- data.frame(lat = rnorm(500, mean = mean_lat, sd = sd_lat),
                      long = rnorm(500, mean = mean_long, sd = sd_long))

trees <- mutate(trees, tree_type = stri_rand_strings(500, 1, pattern = "[A-I]"))

radius = 500

# Convert to sf, set the crs to EPSG:4326 (lat/long),
# and transform CRS to Lambert conformal conic projection - North America
tree_sf <- st_as_sf(trees, coords = c("long", "lat"), crs = 4326) %>%
  st_transform(102009)

# Buffer circles by 100m
tree_circles <- st_buffer(tree_sf, dist = radius)
enclosed_trees <- st_intersection(tree_sf, tree_circles)

# Transform back to 4326 to make Leaflet happy
tree_circles <- st_transform(tree_circles, 4326)
enclosed_trees <- st_transform(enclosed_trees,4326)

leaflet() %>%
  addTiles() %>%
  addPolygons(data = tree_circles) %>%
  addCircles(data = enclosed_trees)


texts: Projection error, geotransform not valid

texts:
I am trying to extract data from gridded weather datasets such as GLDAS from pixels in which weather stations are located.  I created an ArcGIS shapefile  of weather station points.  These weatherstation point coordinates were in GCS_North_American_1983.

My GLDAS raster data has a strange geographic coordinate system (see below).  To properly extract the pixel values of my raster underneath the points in my shapefile I need to reproject the shapefile to be in the same GCS as the GLDAS rasters. When I attempt to project the weather station points to the GLDAS GCS I get the following error: "000151 : The geotransformation(s) is not valid."

I am using the Project tool in ArcMap.  I select the input coordinate system to match the data to be reprojected (GCS_North_American_1983) and I select the output coordinate system to match my GLDAS raster.  The spatial reference for the raster data is as follows:

  GCS_Coordinate_System_imported_from_GRIB_file
  Authority: Custom

  Angular Unit: Degree (0.0174532925199433)
  Prime Meridian: Greenwich (0.0)
  Datum: D_unknown
   Spheroid: Sphere
     Semimajor Axis: 6367470.0
     Semiminor Axis: 6367470.0
     Inverse Flattening: 0.0

I did not create that GCS, that's what the raster came with.  Projecting to the above GCS fails yet I have been able to project the same weather station points to this subtly different output coordinate system from a similar gridded weather dataset called NLDAS:

  GCS_Coordinate_System_imported_from_GRIB_file
  Authority: Custom

  Angular Unit: Degree (0.0174532925199433)
  Prime Meridian: Greenwich (0.0)
  Datum: D_unknown
   Spheroid: Sphere
     Semimajor Axis: 6371200.0
     Semiminor Axis: 6371200.0
     Inverse Flattening: 0.0

I can't seem to figure out why I get an error for the first GCS projection and not the second.


texts: How can I add labels to charts in Mapserver?

texts:
I have created a mapfile with a piechart and I want to add labels to the various pieces of the pie.

The layer in the mapfile is defined as follows:

LAYER
DATA "Test"
NAME "Test"
METADATA
    WMS_TITLE "Test"
    WMS_ABSTRACT "Test"
    WMS_SRS "EPSG:28992"
    WFS_TITLE "Test"
    RESULT_FIELDS "Test"
END
SIZEUNITS PIXELS
STATUS ON
TOLERANCE 10
TOLERANCEUNITS PIXELS
DUMP TRUE
TYPE CHART
PROCESSING "CHART_TYPE=PIE"
PROCESSING "CHART_SIZE=30"
UNITS METERS
MAXSCALEDENOM 150000
LABELITEM "item"
CLASS
    NAME "object1"
    STYLE
        SIZE [size]
        COLOR 44 162 95
        OUTLINECOLOR 0 0 0
    END
    LABEL
        COLOR 0 0 0
        OUTLINECOLOR 225 225 0
        #MINDISTANCE 500
        MINFEATURESIZE AUTO
        #OFFSET 0 0
        #ANGLE AUTO
        TYPE TRUETYPE
        FONT arial
        SIZE 2
        #ANTIALIAS TRUE
        #POSITION CC
        #PARTIALS FALSE
    END
END
CLASS
    NAME "object2"
    STYLE
        SIZE [size]
        COLOR 69 130 189
        OUTLINECOLOR 0 0 0
    END
END
END

The mapfile is working, I see the piecharts. However, the labelling is not working.

Does anyone know how I can make this work? Or is this not possible?


texts: Is it possible to have R write out .sbn and .sbx files?

texts:
I'm working on an situation where a client wants, in addition to the .shp & .prj files that I can produce with e.g. rgdal::writeOGR, the .sbn & .sbx files as well.

My understanding is this is an ARCGIS proprietary file type (?), as explicated here, which is designed to facilitate lookups/queries on larger shapefiles (which indeed I am dealing with).

Is this truly a proprietary file type, or is there some way to produce these indices using R tools?


texts: I am having an (.mdb) file .I want to export this file to (.dgn) format .How can i do this in Qgis ? Is there any plugin available ?


texts: Creating simple x y grid with no projection for non geographical spatial analysis in QGIS?

texts:
I want to analyse a part of the body for which I have high resolution images that have Non-geographical coordinates on an x y grid.

How can I create this as a project in qgis without creating crs based problems?

I guess I make my own but I'm not sure of the syntax to use.


texts:
I am trying to find the length of a polyline inside a polygon using this query:
DECLARE @p geometry, @l geometry;

SELECT @l = SHAPE FROM [Vz_DSM_Stationing].[dbo].[Route_Sections] WHERE OBJECTID = 3;
SELECT @p = SHAPE FROM [Vz_DSM_Extra].[dbo].[PERMIT_BOUNDARIES] WHERE OBJECTID = 8;

SELECT @p.STIntersection(@l).STLength();

This query retrieves the shapes correctly, but selects NULL at the final statement.


texts:
Does anyone out there know the geometry storage encoding for ESRI's Runtime Geodatabase (.sqlite) format? aka mobile geodatabase.

There are 'GDB_ tables', 'st_ tables', 'RTREE indexes' and other fairly standard geodatabase bits, but I'm having trouble identifying the binary encoding of the SHAPE column in the feature tables. It's not WKB, it appears to have a very standard header and a consistent 'footer' character (I think).

These are XY points stored in 4326.


texts: Using Data Driven Pages (DDP) to create set of maps with same extent but different layers?

texts:
Is it possible to utilise the Data Driven Pages (DDP) tool in ArcMap or ArcPy to create a set of maps with the same extent but showing different layers?

To clarify I regularly make large sets of maps for projects all in the same location often with the same baselayers but modified to focus on different aspects for the relevant section of a report such as monitoring locations, land-use, key features etc. which all require separate .mxds to be created and maintained.

A DDP type system would allow me to utilise the dynamic text and ensure a standard template / output throughout, it would also make life a lot easier when updating a general figure number or comment as a report is written in one .mxd rather than multiple for a client.

I'm using ArcMap 10.1, but having read through the ArcGIS DDP online help I assume this would have to be done with ArcPy as I've been unable to figure out how to do it within ArcMap. I'm aware that QGIS has the Atlas Plugin, although I'm not as familiar with its capabilities but would happily use it as an alternative.


texts: Truly lossless TopoJSON / keep all coordinates in TopoJSON

texts:
I'm displaying GPX traces on a Leaflet map. I'm using TopoJSON to drastically reduce file size, which works pretty well.

My problem is the following : I want to be able to display altitude and timestamp for each point of the trace. In order to do that, I need the TopoJSON cli to keep ALL of the original coordinates, to match their indexes with indexes of the timestamps array.

But so far I'm unable to do that. On my GPX file converted to GeoJSON, I got 4062 coordinates, and the GeoJSON file converted to TopoJSON has one arc with 3913 vertices.

Using this :

togeojson trace.gpx | topojson --no-quantization -p > trace.topojson

Gives :

bounds: -18.033666 27.732765 -17.981946 27.772151 (spherical)
topology: 2 arcs, 3915 points
prune: retained 2 / 2 arcs (100%)

So, is it possible to keep the same amount of points, so that I could map precisely coordinates in the TopoJSON file and timestamps?


texts: How to generate bounding boxes for Spacenet tif files using geojson vector data

texts:
I have .tif raster files from the SpaceNet dataset, along with corresponding object descriptions in geojson files. I want to be able to generate Darknet style bounding box information ([object_class, x_center, y_center, width, height]) from this data. Note that the class is same for all objects in one geojson file.

I am unable to find any suitable tutorial or guide to perform this transformation. The spacenet utilities guide on GitHub is old and doesn't work. Could anybody please help me with this transformation? I am not sure if rasterio can be used to do this.

Download sample .geojson file here: https://ufile.io/drcjt
Download sample .tif file here: https://ufile.io/wsq92

Edit: Please note that the Darknet style bounding box ([object_class, x_center, y_center, width, height]) assumes a rectangle in output, rather than a polygon.


texts: Is there any library for converting projection transformation in php?

texts:
I would like to transform data from lat/long ETRS89 to y/x UTM ETRS89 through php.

Any idea?


texts: Google Earth Engine, how to distinguish between rivers/streams and ponds/lakes in a water mask

texts:
I am trying to pull out every lake and pond in a region of SW Alaska. I have successfully identified all water features using the following code:

// study area
var rect = ee.Geometry.Rectangle([-168, 58, -155, 64], 'EPSG:4326', false); //xMin, yMin, xMax, yMax
Map.addLayer(rect,{},"rectangle", false);

// Function to mask clouds using the Sentinel-2 QA band.
function maskS2clouds(image) {
  var qa = image.select('QA60');
  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = Math.pow(2, 10);
  var cirrusBitMask = Math.pow(2, 11);
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));
  // Return the masked and scaled data.
  return image.updateMask(mask).divide(10000);
}

// Bring in the Sentinel-2 collection and create MNDWI from a median of the filtered collection
var s2 = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2016-06-01', '2018-10-15')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .filter(ee.Filter.dayOfYear(153, 274)) // summer images only
                  .map(maskS2clouds)
                  .median()
                  .select("B2","B3","B4","B5","B6","B7","B8","B11");
var mndwi = s2.normalizedDifference(['B3', 'B11']).rename('mndwi');
Map.addLayer(mndwi.clip(rect), {min:-1, max:1, palette:['brown','black','blue']}, 'Modified Normalized Difference Water Index MNDWI', false);

// Change this threshold to allow more or less water in mask
var ndwiThreshold = 0.4;
var water = mndwi.gt(ndwiThreshold);
Map.addLayer(water.clip(rect), {min:0, max:1, palette:['black', 'blue']}, 'water', true);
Map.centerObject(rect,9);

but how could I separate lakes/ponds from streams and rivers? I've been trying to use aspect ration but haven't gotten that to work yet.


texts: Updating the Entry in geometry_columns for Spatial Table Created from a Query

texts:
I created a spatial table from a query. Now, I notice in geometry_columns, the following entry for the table:

 f_table_catalog | f_table_schema |    f_table_name    | f_geometry_column | coord_dimension |  srid  |     type
      my_db      | public         | my_table           | geom              |               2 |      0 | GEOMETRY

I want the srid to be 4326 (this is the srid of the original table) and the type to be Linestring (currently Geometry). Running the following,

SELECT Populate_Geometry_Columns('public.my_table'::regclass, false);

updated the srid in geometry_columns to 4326. However the type is still Geometry. How can I make sure that the geom column in my_table is properly recognized in the geometry_columns and its type is Linestring?


texts: Error when installing Nominatim on Ubuntu

texts:
Due to a large amount of geocoding tasks I want to install Nominatim on my local machine following the official manual. I installed the required software without problems, however when installing Nominatim itself as follows

./configure
make

I encounter the following error

middle-pgsql.c:1660:39: error: format specifies type 'long' but the argument has
  type 'unsigned int' [-Werror,-Wformat]
                        size * 8, sizeof(osmid_t) * 8);
                                  ^~~~~~~~~~~~~~~~~~~
1 error generated.
make[2]: *** [middle-pgsql.o] Error 1
make[2]: Leaving directory `/home/jharvard/Nominatim/osm2pgsql'
make[1]: *** [all] Error 2
make[1]: Leaving directory `/home/jharvard/Nominatim/osm2pgsql'
make: *** [all-recursive] Error 1

Though the configuration file run smoothly, it showed the following warning that I don't know how to deal with. Can anybody help? Is that likely to be the problem?

checking for libxml - version >= 2.0.0... no
*** Could not run libxml test program, checking why...
*** The test program compiled, but did not run. This usually means
*** that the run-time linker is not finding LIBXML or finding the wrong
*** version of LIBXML. If it is not finding LIBXML, you'll need to set your
*** LD_LIBRARY_PATH environment variable, or edit /etc/ld.so.conf to point
*** to the installed location  Also, make sure you have run ldconfig if that
*** is required on your system
***
*** If you have an old version installed, it is best to remove it, although
*** you may also be able to get things to work by modifying LD_LIBRARY_PATH

I have a Mac OS 10.10.2 Yosemite on which I am running Ubuntu 14.04.1 LTS (Release 14.04) using Virtual Box 4.3.22 (the Ubuntu distribution I am using is distributed for an edx course).


texts: Can QGIS merge CSV files with different colums?

texts:
I am trying to merge a month's worth of daily CSV files without using ArcGIS. Each file has ~25,000 entries so this is a high data load project.

Sample data can be found here. Note that we are using the Prerun product.

Arc can handle it (using simple merge geoprocess tool), but we have to make this open source as it will be used to teach students who don't have access to ESRI software.

We tried a simple c# code but ran into the problem that some of the days have extra columns, which snafu'd our code.

Before we spend the time to make our code much more robust to handle the disparity in the products, I wanted to ask if there was a way to merge tables (CSV files) within QGIS? This would be the simplest methodology for the people we're trying to benefit.


texts: Extract layer attributes and addresses QGIS

texts:
I have congressional district layer and a csv file with addresses. I'd like to extract all addresses and their corresponding congressional district. How do I do that using Qgis?

The CSV File contains 3 columns including X,Y coordinates. When I import the CSV file with addresses the CRS = WGS84 (EPGS 4326).

The Congressional District layer CRS = +proj=lcc +lat_1=34.03333333333333 +lat_2=35.46666666666667 +lat_0=33.5 +lon_0=-118 +x_0=2000000 +y_0=500000.0000000001 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

I am using QGIS version 2.0.1.

@Morgan Ellingham, after following your steps, I ran into the following error prompt:


texts: How to Remove Dark Bands Where Scenes Overlap in Mosaic Landsat TM5 Image?

texts: I have an image mosaic of 5 Landsat TM5 scenes. Where the scenes overlap there is a dark band, only on the left and right sides, but the top and bottom are seamless. What is the best way to correct this in Erdas 2011?


texts: Viewing ArcPy code from results of process run

texts: How do I copy the ArcPy code from a process I have run in ArcMap where I have run a process with specific environment settings etc (density of, distance to etc) so that I can run the same code iteratively on a number of other features?


texts: Installing 64 bit Python on ArcGIS Desktop?

texts:
I am trying to have Python 64-bit with ArcGIS so that I can use my entire RAM when doing data analysis. Apparently ArcGIS 10.2.2 Desktop comes just with 'Python 32-bit'.

Should I install a fresh Python 64-bit and import arcpy from there?

Is that going to work?

I need to use arcpy and scipy libraries for my data analysis and I need to call both of these libraries in my Python code.

Is there a better way to run 'Python 64-bit' with ArcGIS?


texts: How to embed ArcGIS Server Map in personal Google or Blogger website?

texts:
I have created my own free Google website, and I want to embed a web map on it that I have created with ArcGIS Server.

Can I do this, and if so, how?


texts: How to find the POI (points of interest) Within a Given Radius?

texts:
I have a table with x,y coordinates with the_geom column. Now I need to find the nodes (Points Of Interest) present in a given radius. Kindly help me in knowing how to find it.

Thanks.


texts: Changing order of appearance for feature class with multiple overlapping polygons in ArcGIS Pro?

texts:
I am using ArcGIS Pro.

I have a feature class with many overlapping polygons.   I am changing the transparency of each of the symbols individually so that you can see through each polygon to the underlying features,  much like it would appear if they were 2 different feature classes.

How do I make the underlying polygon appear on top of, instead of underneath, the larger polygon?

Screenshot example embedded below.


texts: While importing *.tif in PostGIS db using raster2pgsql command i am facing issue "ERROR: invalid memory alloc request size 1073741824"

texts:
I just need to import a .tif file into my PostGIS db. the TIFF file I have is 430 MB. While using this command following error is thrown:

 raster2pgsql -I -C -c *.ndvi.tif public.gistest1 | psql -U XXXX -h
 localhost -d gisdb -p 5432

Password for user cropin_suhal: Processing 1/1: *.ndvi.tif

BEGIN
CREATE TABLE
ERROR:  invalid memory alloc request size 1073741824
ERROR:  current transaction is aborted, commands ignored
until end of     transaction block
ERROR:  current transaction is
aborted, commands ignored until end of     transaction block
ERROR:
current transaction is aborted, commands ignored until end of
transaction block    ROLLBACK

The work around I did is to increase max_wal_size (integer)=10GB in postgresql.conf file but of no use.


texts:
What is the best way to join following tables? I have 55 records in the second table and I would like to just "drag and drop" the data to the first one.

Shapefile table

Database table


texts: dynamic visualization of spatiotemporal data on a grid

texts:
I am handling spatiotemporal data and it is a quite new area for me.

Basically I have data related to the position of cars at a given time in Montreal. I would like to dynamically visualize the density of cars in a grid.

I am using R with some appropriate libraries such as plotKML, spacetime, sp, ..., which are useful to link R with Google Earth. I have already manually computed the density of cars and the data looks like that:

                   Time   Latitude   Longitude Dens
793 2016-02-23 23:20:56 45.4963903 -73.5596835 257
183 2016-02-23 04:48:38 45.5471917 -73.5492618 57
141 2016-02-23 04:39:58 45.5471535 -73.5491820 157
258 2016-02-23 05:04:08 45.5357973 -73.5593413 13
846 2016-02-23 23:31:56 45.5031632 -73.5791742 420
236 2016-02-23 04:59:33 45.5471562 -73.5492615 269
128 2016-02-23 04:37:18 45.5470988 -73.5492218 11
886 2016-02-23 23:40:11 45.5106078 -73.5786010 50
524 2016-02-23 22:25:34 45.4961470 -73.5591993 189
741 2016-02-23 23:10:21 45.4961615 -73.5591900 57

Now I should concretely visualize it in some way. I was thinking to use a grid, and plotting bubbles for each cell, with radius proportional to the density, or nothing if in that cell at that given time the density is 0.

I think most of the work can be done by properly using the plotKML function, since in some case it is possible to "intelligently" display bubbles (as shown in the wiki tutorial of plotKML here). For example, if we get rid of the Time variable, one (static) solution is the following:

library(sp)
library(plotKML)
library(spacetime)

# from raw data in df to "clean" data
mydf = createmydf(df)

spdf <- SpatialPointsDataFrame(
  coords = mydf[,c("Longitude","Latitude")],
  data = data.frame(mydf[,"Dens"]),
  proj4string = CRS("+proj=longlat +datum=WGS84"))

# static plot:
plotKML(spdf, LabelScale = .1, colour_scale=rep("red", 2))

And the output for this is this (static) picture:

Now the question is: is it possible to do the same but in a dynamic way (i.e. considering the temporal dimension)?.

I tried to use the STFDF function. For this I computed a full space-time grid with n spatial locations and m times. This represents the input of STFDF, together with n x m observations (Densities in my case). This is the code:

# gridSpace: n x 2 matrix with Latitude and Longitude positions
# gridTime: a sequence of times, a vector of length m
# Dens: vector of length n x m
sp <- SpatialPoints(gridSpace[c("Longitude","Latitude")])
proj4string(sp) <- CRS("+proj=longlat +datum=WGS84")
mydf.st = STFDF(sp, gridTime, data.frame(Dens))

plotKML(mydf.st)

In this case the plot is fine but there are just points of fixed dimension and changing the color with the time. Does somebody know a way to change this representation with bubbles of dynamic radius or bars? I think it should be related to the plotKML function. Or should I switch to another time function such as STIDF?


texts:
I'm trying to calculate the VCI data for my entire State.

I'm using MODIS's 16 day NDVI Average. From this post, I see that I would need to calculate the Min & Max NDVI values for each pixel based on several year's data.

My state is covered by 4 different MODIS Scenes, and hence downloading and processing this area for multiple years would be a long and resource intensive process.

Is there any Dataset which gives the Min & Max NDVI Values based on long-term data?


texts: Error occurring when converting to dB during the python-based classification

texts:
I am trying to perform an object-based classification using a python-based script with rsgislib. Although, when converting to dB, an error is appearing. Any ideas why is this happening?
 Create and Calculate dB columns
ratDataset = gdal.Open(clumpsImg, gdal.GA_Update)
HHPow = rat.readColumn(ratDataset, 'HHPow')
HVPow = rat.readColumn(ratDataset, 'HVPow')
HHHVPow = rat.readColumn(ratDataset, 'HHHVPow')

HHdB = numpy.where(HHPow>0, 10*numpy.log10(HHPow), 0)
HVdB = numpy.where(HVPow>0, 10*numpy.log10(HVPow), 0)
HHHVdB = numpy.where(HHHVPow>0, 10*numpy.log10(HHHVPow), 0)

rat.writeColumn(ratDataset, "HHdB", HHdB)
rat.writeColumn(ratDataset, "HVdB", HVdB)
rat.writeColumn(ratDataset, "HHHVdB", HHHVdB)
ratDataset = None

And I keep getting the following error:
Writing Stats (Mean, ) to Output RAT
/data/3_to_5_copy/16_PerformObjClass.py:63: RuntimeWarning: divide by zero encountered in log10
HHdB = numpy.where(HHPow>0, 10numpy.log10(HHPow), 0)
/data/3_to_5_copy/16_PerformObjClass.py:64: RuntimeWarning: divide by zero encountered in log10
HVdB = numpy.where(HVPow>0, 10numpy.log10(HVPow), 0)
/data/3_to_5_copy/16_PerformObjClass.py:65: RuntimeWarning: divide by zero encountered in log10
HHHVdB = numpy.where(HHHVPow>0, 10*numpy.log10(HHHVPow), 0)
Creating output image using input image


texts: ArcGIS Online - Legend for raster/image (in web map, web map app, via Image service or other)

texts:
Is there a way to have an AGO map, or mapping application, with rasters than be symbolised in the legend?

We currently use Image Services which load fine but cannot be listed in the legend - so the symbology becomes a little useless.

Is there a way to attach legend properties to the image service?

If we were to switch to a WMS would that help?


texts: Can you auto-populate a field in Survey123 from another field by only pulling certain characters?

texts:
I have a field where the user inputs text Aa1Bb2Cc3

Am I able to auto populate another field by only populating the field with Bb or 3.

Auto populate the field with characters [4:6] or just [9]

How would I do this?


texts: What's the meaning of 'fb' at the end of PyQgis Algorithm ID

texts:
I'm following this PYQGIS tutorial, the algorithm ID of zonal statistics in my qgis is "native:zonalstatisticsfb".
The algorithm doesn't work unless I delete the fb at the end of the ID.
So, I'm wondering what the fb at the end of the ID means.


texts: QGIS: Making changes to vector layers stored in PostGIS

texts:
I have a vector layer stored in PostGIS. I split several lines in the layer and simplified others by removing vertices.

When I tried to save the changes to the layer I got postgres more postgres errors about duplicate gids.

Is this a result splitting features and if so how do I get around the issue?

Russell


texts: How to find forget password in QGIS DB manager?

texts:
I am new to QGIS. I forget the DB manager password

.

Is there a way to fix it?


texts: Filtering features outside features of another layer using QGIS geometry generator

texts:
I am trying to filter features in one layer outside features in another layer using geometry generator.
I have a buildings layer and an area layer and now I want to only show the buildings outside the areas.
The buildings layer is a polygon layer and the area layer is also a polygon layer consisting of several areas (features).
Since it is a many feature area layer I think I need to use an aggregate function e.g. collect() and/or aggregate() from the QGIS expression builder.
I trying the following expression on the geometry generator on the buildings layer...
difference(
    $geometry, collect(geometry
                            (get_feature('layer_id','attribute','value')
                                )
                            )
                       )

...which ends up in an syntax error message.


texts: ERROR 000824 while trying to run IsNull Method outside ArcGIS?

texts:
I am trying to perform IsNull function to a raster using PyScripter. This is the code I use:

import arcpy
import arcinfo
InRaster = "C:\\data\InRaster"
IsNullRa = "C:\\data\IsNullRa"
arcpy.gp.IsNull_sa(InRaster, IsNullRa)

In ArcGIS Python window, it works OK. However, in a different editor such as PyScripter I get this error:

ERROR 000824: The tool is not licensed. Failed to execute (IsNull).

I have ArcGIS 10.1 and do have spatial analyst. I tried installing  Microsoft Visual C++ 2008 Redistributable Package (x64) from here but it did not work. Do you know how I can run IsNull outside of ArcGIS? Thanks.


texts: I wonder to know if -3000 in this MOD13Q1 product represents missing value? I saw the product which describes the fill value as -3000 and No data value as N/A. I am not sure whether these two values all represent missing values or the fill value of -3000 represents the other thing.


texts: Fixed length line segment in QGIS: why does helping circle not appear always?

texts:
I am working with fixed line lengths. Advanced Digitizing Panel is activated and I set the line length. Usually there is a the blue "helping" circle whose radius is the value entered in the distance text box (see picture).

But sometimes there is no circle. Why is that?

I would like to always see the blue circle.


texts:
I have a raster of estimated soil moisture values and 30 points of field measurements.

I wanted to correlate the measured and estimated values using scatterplot but the result in not good. There might have been errors while field measurements and my map has a high resolution of 10cm, so I decided to calculate an average of 10 pixels in each direction for each point so I have an average value of 1m for each point which hopefully will increase the result accuracy.

How can I calculate just an average of specific pixels in ArcMap 10? I know about zonal statistics but here I dont have another raster dataset. I have to chose the pixels I need around each point.


texts: How to store Google Maps DrawingManager Circle overlay in PostgreSQL Postgis GEOMETRY(GEOMETRY, 4326)

texts: What is the correct way to store a Google Maps DrawingManager circle overlay? Should I just store the center as a POINT and the radius in a separate column or store a polygon using ST_Buffer()?


texts: Dealing with overlapping features in OpenLayers

texts:
If users create a polygon feature that completely overlaps another one, like this …

… then there's no way you can select the inner feature, i.e. the one that's overlapped by the other. I'm currently using a simple SelectFeature like this:

var select = new OpenLayers.Control.SelectFeature(
  [this.annotationLayer, this.controlPointsLayer], {
    hover: true,
    renderIntent: "temporary"
    }
);
this.map.addControl(select);
select.activate();

But obviously, I can only get it to select the outer feature. What is the typical way of dealing with this? How can I make it possible for the user to select the inner feature as well?

The simplest solution would probably be to make the outer feature only selectable within the area that's left when subtracting the inner feature – but I have no idea how to go about doing that.


texts: GeoNetwork: How to reset the admin password?

texts:
I am unable to log on with the admin account to a certain GeoNetwork
instance, but I do not recall ever changing the password. How can I reset (or
change) the admin password without logging on to the web interface?


texts: Testing and training points can not be visualized in Google Earth Engine

texts:
I learned supervised classification through the example of Google Earth Engine, but in this example, I found that I cannot use Map.addLayer () to visualize the training data and test data. When I use this code:

Map.addLayer (trainingPartition, {'color': '# ff0000'}, 'train');
Map.addLayer (testingPartition, {'color': '# 3b8b00'}, 'test');

Nothing appears on the screen.
Does anyone know why this is the case?
Here are  Key codes:

// Optionally, do some accuracy assessment.  Fist, add a column of
// random uniforms to the training dataset.
var withRandom = training.randomColumn('random');
// We want to reserve some of the data for testing, to avoid overfitting the model.
var split = 0.7;  // Roughly 70% training, 30% testing.
var trainingPartition = withRandom.filter(ee.Filter.lt('random', split));
var testingPartition = withRandom.filter(ee.Filter.gte('random', split));
Map.addLayer(trainingPartition,{'color':'#ff0000'},'train')
Map.addLayer(testingPartition,{'color':'#3b8b00'},'test')

Here is link: https://code.earthengine.google.com/1c2f985baf131e88b3dd96fbbab9b588


texts:
My county ( Pima, Arizona ) uses a number called “average cross slope” which they define as

  ACS = (I * L)/A

where

I = contour interval
L = length of contour lines in the area
A = area

Given a unit square with corners (0,0) (1,1), a contour interval of 1 and two contour lines (0,0)(0,1) and (1,0)(1,1), I get
L=2,I=1,A=1 and ACS = (I*L)/A = 1*2/1 = 2

Using the same unit square, a contour interval of 0.5 and three interval lines (0,0)(0,1), (0.5,0)(0.5,1) and (1,0)(1,1) I get

ACS = (I*L)/A = (0.5*3)/1 = 1.5

Since these two examples represent the same topography, I would expect the result to be the same.

Also, the county refers to average cross slope as a percentage. Intuitively, since the topography is a 45 degree angle running left to right, I would think the average slope would be 1. (slope = rise/run), or 100%

So I have three different answers for the same topography, 2, 1.5 and 1.

I’m not a GIS person, I’m a software guy and owner/builder.

Is this a common equation in GIS?


texts: Clustering of polygons on map using Leaflet?

texts:
There is a plugin for Leaflet to cluster the markers.The result looks like so. Result:

If I like something to do to show clusters of circles and polygons?

For example, when you hover over a cluster it shows the ground on which he considers the number. I would like to have this mesh of polygons is always on display and she was the color of the color circle. Something in concept. Wanted result:

Any idea experts this system?

Is it even possible to do this?


texts: Adding text in integer field for label using QGIS?

texts:
I am labeling mile points along several different lines in QGIS 3.2. For the labels I am using a field that is just integers so the labels are 1, 2, 3, etc but I would like to add some text so that that for one line it'd say, for example A1, A2, etc.

Is there a function to add the text characters in either the field/a new field or in the label styling window?


texts: I am trying to generate random points on a route from one point to another. What I have till now is, given a source and a destination, I can get a route using Google's API. But how do I get the coordinates of random points on the route?


texts: ARCGIS Calculate Value Tool - Always returns same value

texts:
I want to check if an environmental shapefile was generated from a previous process.

I utilized the Calculate Value tool in Modelbuilder but the problem is that no matter if this Wetlands layer exists or not, it always returns "true". I've tested the code in IDLE and it runs with no issues.

Expression

  x("%Wetlands%", "%Workspace%")

Code Block

import arcpy
def x(Wetlands, Workspace):
 arcpy.env.workspace = "C:/Omitted/Environmental.gdb"
 if arcpy.Exists("Wetlands"):
   return "true"
 else:
   return "false"

Data type

  Boolean


texts:
I have an area, example 4000m x 4000m. So, 4x4km. The origin can be placed anywhere, maybe in center, maybe on the top-left, nevermind.

Now, I want to fill this area with some stuff from OSM. I read a lot about UTM, zones, etc. UTM quadrant is divided into 6x10 100km squares, this squares are divided into 10x10 10x10km squares, this squares are divided to 10x10 1x1km squares and so on.

Getting back to my 4kmx4km square. How do I convert coordinates to place objects in my square and the position responds to reality?

To make it harder, let's say that I want to put on the square objects from squares from different UTM quadrants from different zones, for example: most southern-east 1x1 from 34V, most southern-west from 35V, most northern-west from 35U and most northern-east from 35U.

Additional question: Is it possible to create one continuous maps (builded from smaller pieces), which covers the region composed of several UTM zones (for example, from 29W to 34S).

How to translate UTM coordinates properly?

EDIT
Ok, I implement Albers, Lambert and Equidistant projections for some tests. Input args: a = 6378206.4 m, e = 0.00676866, e2 = 0.0822719 Standard parallels: 55.0, 49.0 Origin: 52.205746, 22.514281

And results:

Albers : x = -489.71509; y = -2291.5

Lambert : x = -489.15250; y = -1653.5

Equidistant: x = -489.71533; y = -2289.5

I checked the algorithms and are ok, the results agree with the sample data. What is the reason that the result of the Lambert projection stands out from the rest? Does this mean that Lambert is better to applied to larger areas than the other two?


texts:
I'm using this ESRI leaflet plugin Example to load a tiledMapLayer.

<script>
  var map = L.map('map').setView([30.70, -81.47], 12);

  L.esri.tiledMapLayer({
    url: "https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer"
  }).addTo(map);
</script>

available here: https://esri.github.io/esri-leaflet/examples/tile-layer-1.html

This works when you have only one titled map service on the URL or just number 0 which it defaults to, but I would like refer to the other titled services on this url.

World Imagery (0)
Low Resolution 15m Imagery (1)
High Resolution 60cm Imagery (2)
High Resolution 30cm Imagery (3)
Citations (4)

I was expecting to be able to write code like this but no luck:

<script>
  var map = L.map('map').setView([30.70, -81.47], 12);

  L.esri.tiledMapLayer({
    url: "https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/2"
  }).addTo(map);
</script>

I've appended the /2 at the end of the service expecting to load High Resolution 60cm Imagery (2) from the service but I get no returns.


texts: Appending rows to empty FeatureClass with ArcPy

texts:
My goal is to delete all rows in an existing feature class, then append new rows from a new shapefile... everything works except for the append. The table-view portion as the input is giving me the trouble. I am working in ArcGIS Desktop 10.5.1

I got all the code working fine except for the append, it produces this error:

Traceback (most recent call last):
  File "C:\Users\jhead\Desktop\UpdateCursor_Append.py", line 30, in <module>
    arcpy.Append_management(sourceView, emptyFC, 'TEST')
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 4271, in Append
    raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Table View.
ERROR 000840: The value is not a Raster Layer.
Failed to execute (Append).

My code is below:

emptyFC = r"Database Connections\TEST ENVIRONMENT.sde\TEST_ENVIRONMENT.DBO.featureclass"

sourceSHP = r"path\to\file.shp"

sourceDBF = r"path\to\output\folder"
print 'Converting new parcels to DBASE'
arcpy.TableToDBASE_conversion(sourceSHP,sourceDBF)

print ' Converting was successfull!'

print 'creating temporary table view in memory...'
sourceDBF2 = r"path\to\file.dbf"
sourceView = arcpy.MakeTableView_management(sourceDBF2, "sourceView")

print ' Appending new rows to existing Feature Class...\n'

arcpy.Append_management(sourceView, emptyFC, 'TEST')

According to the documentation I should be able to use a shapefile as an input data set, but using the shapefile also produces the exact same error as the table view.

http://desktop.arcgis.com/en/arcmap/latest/tools/data-management-toolbox/append.htm


texts: Can't close custom popup at fixed position

texts:
I created this custom popup at a fixed position. When clicking a marker a popup opens, I'm able to close the popup and open another one clicking a different marker, then I'm not able to close the popup. The popup wrapper remains open and only changes it's content when clicked on a different marker? The eventListener works partially. Need that piece of code that does the job.

var redFlag = L.icon({
iconUrl: 'images/mapmarker2.png',
iconSize: [34, 34],
iconAnchor: [17,34]
});

const myLayer = L.geoJSON(art, {
    pointToLayer: function (feature, latlng) {
          return L.marker(latlng, {icon: redFlag});

},
onEachFeature: function ( feature, layer) {
    layer.on('click', function(e){

        var getWrap = document.getElementById('wrapper');
        var wrap = getWrap.appendChild(document.createElement('div'));
        wrap.className = 'wrapper';
        wrap.innerHTML =
         `<div class="close">X</div>`+
         `<div class="popUpContent" style="background-color:#e8f4ff">` +
         `

`+
         `
${feature.properties.name}
`+
         `<span class="puTitle">"${feature.properties.title}"</span><br>` +
            `<div class="extra3">${feature.properties.extra}</div></div>`+
            `</div>`;

    if(!feature.properties.title){

        wrap.innerHTML =
            `<div class="close">X</div>`+
            `<div class="popUpContent" style="background-color:#e8f4ff">` +
            `

`+
            `
${feature.properties.name}
`+
            `<div class="extra3">${feature.properties.extra}</div></div>`+
            `</div>`;

            document.querySelector('.close').addEventListener( 'click',
            closePopup);
            function closePopup(e){

              if(event.target === 'close'){
                  document.querySelector(".wrapper").style.display='none'
              }else
     if(document.querySelector(".wrapper").style.display='block'){
                  document.querySelector(".wrapper").style.display='none';
                 }
              }
           }
       });
   }

});

mymap.addLayer(myLayer)


texts: Mapping several cities with Locator from ArcGIS Online?

texts:
I'm trying to map several points in the map, all of them are cities, I used the example from the example in ArcGIS Geocode an Address, and I locate it well, but I just can map ONE at a time.

In this section :

public var lugares:Array = new Array("GUADALAJARA", "ECATEPEC"); //Cities from Mexico
var parameters:AddressToLocationsParameters = new AddressToLocationsParameters();
//parameters such as 'SingleLine' are dependent upon the locator service used.
//parameters.address = { SingleLine: onelineaddress };
parameters.addresses = lugares ;
// Use outFields to get back extra information
// The exact fields available depends on the specific locator service used.
parameters.outFields = [ "Loc_name" ];
locator.addressesToLocations(parameters, new AsyncResponder(onResult, onFault));

when I changed the instructions to plural parameters.addresses = lugares and locator.addressesToLocations It displays an error that a Token must be added and 499 error code,I added token:Object = null but nothing, that error I googled it and it says that the locator service may have permissions... Does anyone has another link? I can't find it, or another solution to map cities? There is going to be more than 2 cities, (just 4 my test I used 2, but it will be maybe hundreds obtained from DB).


texts: QuickOSM TimeOut with even the smallest area in QGIS 3

texts: QuickOSM worked beautifully in QGIS 2 but since I changed to QGIS 3 it never works. It always "timeout" for anything even in the smallest Canvas of a few houses. When I Query for "highway" it does Timeout. I tried with QGIS 3.2, 3.4, 3.8. Multiple servers. Always the same. What am I missing?


texts: Search Sentinel-2 images using python api

texts:
Based on the documentation here, I wrote the following script:

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
from datetime import date
api = SentinelAPI('user name', 'password', 'https://scihub.copernicus.eu/dhus')
polygon = r"C:\fields\filed1.geojson"
footprint = geojson_to_wkt(read_geojson(polygon))

products = api.query(footprint,
                     date = ('20170801','20170830'),
                     platformname = 'Sentinel-2',
                     cloudcoverpercentage = (0,30))

products_df = api.to_dataframe(products)

I got the expected result. I'm trying to do the same thing but the search needs to be based upon x,y coordinates.

I read here that I just need to provide the x,y coordinate (as a tuple) to the footprint so I tried:

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
from datetime import date

api = SentinelAPI('user name', 'password', 'https://scihub.copernicus.eu/dhus')
footprint = (41.9000, 12.5000)

products = api.query(footprint,
                     date = ('20170801','20170830'),
                     platformname = 'Sentinel-2',
                     cloudcoverpercentage = (0,30))

products_df = api.to_dataframe(products)

and I get an error:

SentinelAPIError                          Traceback (most recent call last)
<ipython-input-56-98f8dc400863> in <module>()
      8                      date = ('20170801','20170830'),
      9                      platformname = 'Sentinel-2',
---> 10                      cloudcoverpercentage = (0,30))
     11
     12 products_df = api.to_dataframe(products)

C:\Users\Ran\Anaconda2\lib\site-packages\sentinelsat\sentinel.pyc in query(self, area, date, raw, area_relation, order_by, limit, offset, **keywords)
    132                           order_by, limit, offset, query)
    133         formatted_order_by = _format_order_by(order_by)
--> 134         response, count = self._load_query(query, formatted_order_by, limit, offset)
    135         self.logger.info("Found %s products", count)
    136         return _parse_opensearch_response(response)

C:\Users\Ran\Anaconda2\lib\site-packages\sentinelsat\sentinel.pyc in _load_query(self, query, order_by, limit, offset)
    240
    241     def _load_query(self, query, order_by=None, limit=None, offset=0):
--> 242         products, count = self._load_subquery(query, order_by, limit, offset)
    243
    244         # repeat query until all results have been loaded

C:\Users\Ran\Anaconda2\lib\site-packages\sentinelsat\sentinel.pyc in _load_subquery(self, query, order_by, limit, offset)
    282                 # returned as the total results value when the query string was incorrect.
    283                 raise SentinelAPIError(
--> 284                     'Invalid query string. Check the parameters and format.', response)
    285             total_results = int(json_feed['opensearch:totalResults'])
    286         except (ValueError, KeyError):

SentinelAPIError: HTTP status 200 OK: Invalid query string. Check the parameters and format

How can I perform that?


texts: In cartodb, how can I insert a new row in a table via API?

texts:
We've discovered cartodb and we want to use it in our project, but in the API documentation we've only seen the sql.execute() method (docs) in order to insert a now row in a table.

Is this the way to insert points dynamically in our map?


texts: Using AJAX to run PostGIS SQL statements and display results via GeoServer

texts:
I am learning Web GIS, and I can't find a good example that gives a big picture of a web interface of running spatial analysis using PostGIS.

What I am looking for is a graphical/GIS version of pgAdmin running on the web. Given a SQL statement, the web interface calls PostGIS, possibly generating a result table, which then gets rendered in a graphical/cartographic form rather than in a tabular form as in pgAdmin.

I looked around, and most question/answers here is regarding how to display an existing table in PostGIS. But there is nothing that covers the complete loop. That is, start query, store the result in the database, display it, and maybe using the result table for a new query.

I have a little bit of familiarity with GeoServer and its RESTful API, and very basic knowledge of JavScript/jQuery/HTML5.

Can some suggest a viable route for achieving this query-store-display loop?

Note: I am less interested in the visualization part and getting GeoJson's because they cannot be directly used by PostGIS for further analysis.

From @John's answer, GeoJSON seems to be a very good option and I think there is good reason why GeoJSON is so popular. On the other hand, I was hoping to find a more "server-side" solution that relies primarily on PostGIS doing the spatial analysis and GeoServer doing the rendering/cartography of the PostGIS tables/views. Hopefully, this will involve minimal client side programming (other than issuing SQL commands and GeoServer API calls).

I can do this manually with pgAdmin together with the GeoServer administration interface or RESTful API curling. But I am not sure how to do the same programmatically from the Web with e.g. AJAX and REST APIs. I have tried to use WPS, but without much success. For example, I have tried to perform simple raster map algebra, but the WPS in GeoServer takes two input rasters from server and only gives you an output image for download, not another raster dataset for further analysis. Here, the "loop" is broken.

Besides, WPS functions vary from one package to another. If there is something wrong in a WPS, it's less easy to get it fixed unless the upstream team is really active.

So I am thinking to not rely on WPS, and instead rely on PostGIS to do all vector and raster analysis, as I think Postgres is more tested and the relational/map algebra therein is complete. And I was hoping to use AJAX as a "command center" to issue SQL statements to Postgres (maybe with some server side PHP proxy??), and to render some of the resulting tables to the web via GeoServer (or another map server) (e.g. by publishing/updating the result tables/views).


texts: Google Geocoder Lookup get postal code by country and city

texts:
I'm currently using Google Geocoder API to determine location data for my web service.
It's important for me to get the postal code for a given country and city, but it's look like that Google is limiting the result.

Here's the API call:

http://maps.googleapis.com/maps/api/geocode/json?address=deutschland,+saarlouis&language=de&sensor=false

The result is like that:

{
   "results" : [
      {
         "address_components" : [
            {
               "long_name" : "Saarlouis",
               "short_name" : "Saarlouis",
               "types" : [ "locality", "political" ]
            },
            {
               "long_name" : "Landkreis Saarlouis",
               "short_name" : "Landkreis Saarlouis",
               "types" : [ "administrative_area_level_3", "political" ]
            },
            {
               "long_name" : "Saarland",
               "short_name" : "SL",
               "types" : [ "administrative_area_level_1", "political" ]
            },
            {
               "long_name" : "Deutschland",
               "short_name" : "DE",
               "types" : [ "country", "political" ]
            }
         ],
         "formatted_address" : "Saarlouis, Deutschland",
         "geometry" : {
            "bounds" : {
               "northeast" : {
                  "lat" : 49.36187570,
                  "lng" : 6.815549499999999
               },
               "southwest" : {
                  "lat" : 49.26046930,
                  "lng" : 6.67501510
               }
            },
            "location" : {
               "lat" : 49.31346060,
               "lng" : 6.752286499999999
            },
            "location_type" : "APPROXIMATE",
            "viewport" : {
               "northeast" : {
                  "lat" : 49.36187570,
                  "lng" : 6.815549499999999
               },
               "southwest" : {
                  "lat" : 49.26046930,
                  "lng" : 6.67501510
               }
            }
         },
         "types" : [ "locality", "political" ]
      }
   ],
   "status" : "OK"
}


texts:
I have a raster file of 1600 km2 with cell size 2x2 m that I need to polygonize and the script takes a really long time to finish (>25 minutes, but I don't know how long it effectively takes, since I had to kill the task)

How can I make smaller tiles of the large raster in QGIS? Let's say we make 16 tiles, so 4x4 tiles. What is the best method for that in QGIS 3.3+?


texts: Using table joins with ST_DWithin to optimise process?

texts:
I am trying to select all segments (among millions of them) that are less than X metres from a few thousands polygons obtained from a raster dataset (meaning they've got quite a complex geometry - many vertices) OR less than X metres from another set of thousands polygons (less complex geometry). I am attempting to do so using ST_DWithin and a double LEFT JOIN, but it's proving to be very slow (more than 48 hours running and it hasn't finished yet).

Projected Coordinate System (British National Grid)

All 3 tables have a spatial index

This is the query I'm trying with:

CREATE TABLE result AS
SELECT t1.segments_id , t2.polyg1_id , t3.polyg2_id
FROM segments t1
LEFT JOIN polyg1 t2
ON ST_DWithin(t1.geom, t2.geom, X)
LEFT JOIN polyg2 t3
ON ST_DWithin(t1.geom, t3.geom, X)

Once this query is finished I would just need to discard all records from result that have no polyg1_id or polyg2_id values.

Is there any better (more efficient) ways to achieve what I'm trying to obtain?


texts: How to alter the presentation of FORMA Raw Output FIRMS in Earth Engine?

texts:
I'm attempting to view the FORMA Raw Output FIRMS dataset in Google Earth Engine. The thumbnail for the dataset shows it looking mostly transparent with green areas where the deforestation has occured - this is the image I would like to view and export.

However, when I run the code from the Earth Engine Snippet, it results in mostly black with the areas that were green in white instead.

What would I need to change about the code to result in the first image, and then export it into a file that I could view in Google Earth Pro?


texts: Reshaping polygon shapefile using ArcGIS for Desktop?

texts:
I have a polygon shapefile of my country, but the borders of it are not accurate (according to Google Earth maps borders).

I searched the web and I could not find an accurate shapefile for the borders. I drew the borders on Google Earth as polylines, then I imported them into ArcGIS 10.1 for Desktop.

Now, the old shapefiles of my country are very good, except for the borders line.

How do I cut and reshape the old polygon according the new polyline?

Some places need extending while some others need to shrink.


texts: How can I use Python to save several MXDs within the same script?

texts:
I have a program (720 lines of code) that runs multiple tasks on Arcmap .mxd files. Such as update layers, update text, add layers, etc.

Within my program I have multiple "mxd.save()" instances. All of which occur within for loops.

For example:

for item in mxds:
    mxd = arcpy.mapping.MapDocument(item)
    df=arcpy.mapping.ListDataFrames(mxd,"Project Area")[0]
    print "HOWDY!"
    mxd.save()

Question is can I just have at the very end of my program have one for loop that saves all the mxds I was working on, while still retaining all the changes?

UPDATE:

Problem I think am having is this. I have one of my for loops that adds a layer to the map. If I do not mxd.save() after that layer is added. My program later when it attempts to find the layer that was added cannot find it. I guess it's not held in memory.

UPDATE2:

THANKS! I was doing something similar, but wrong none the less. Each time I had a for loop I would create a new .MapDocument object which I realized (thanks to you) was unnecessary. I now define the mapdocument object once at the start of my program and work from there. Do this has also greatly increased the speed of my program.


texts:
I'm writing a QGIS plugin using PyQGIS and I have a weird problem.

When executing following line of code:

processing.runandload("qgis:createpointsalonglines", r"Clipped", 1, 1, 0, r"memory:nodes")

I get this error

However, I use this "kind of code" multiple times in my plugin, and only get the error for this line...

For example,

I execute

processing.runandload("qgis:clip", "parallel", "polygon", r"memory:Clipped")

and

processing.runandload("qgis:orientedminimumboundingbox", "polygon", 1, r"memory:Oriented_MBBox")

without any problem, in the same program.

Anbybody has any idea why ?

Also: if I define the output layer as "memory:LayerName", the layer is created but the real layer name is not equal to LayerName... QGIS always gives the layer a standard name, depending on the function..

How can I give the layer the name I want?


texts:
I have a tif file of 7378 X 6923(widthXheight) resolution. I resized it to 1200 X 1200 pixels using Image class in python.

 img = Image.open('path to tif file')
 img = img.resize((1200,1200))

But after resizing, the pixel values show unusual values. Before resizing, a certain pixel's value was 125.28. Now, after resizing, that pixel's value has turned to 1123715645. I am reading the pixel value with:

with rasterio.open('path to tif file') as dataset:
    band = dataset.read(1)
    print(type(band))
    value_of_interest = band[row_of_interest, column_of_interest]
    print(value_of_interest)

Why does this happen?

How can I preserve the height values in same pixel position?


texts: Layer group actions (e.g. transparency control)

texts:
I know this question has been asked already... long time ago, though.
When loading a large set of raster data and grouping them in QGIS i miss the opportunity to set basic properties for that group, such as transparency or contrast. Is there a more comfortable way to achieve this functionality apart from using console statements?


texts: Adding custom widget into Mapguide Maestro 4.0.1?

texts:
I am developing a Web-Mapping application using Mapguide Maestro 4.0.1, in that I need to add a custom widget.

I'm looking for any suggestions for how and where to add the custom tools in the Mapguide Maestro. Is there any documentation for how to add the custom tools?

I am using MapGuide OpenSource 2.2, MapGide Maestro 4.0.1.


texts:
I wanted to fully uninstall QGIS (previously installed by
osgeo4w). The problem is that it does not appear on the list of programs to uninstall neither in the control panel>programs and features or in settings>apps & features. What I'm I missing?


texts: Error when kriging: "output variance prediction raster"

texts:
I am trying to run a python script for kriging, but keep getting a following error message. It doesn't make any sense as the output variance prediction raster is optional.

The code:

# Create KrigingModelOrdinary Object
lagSize = 70000
majorRange = 250000
partialSill = 180000
nugget = 34000
kModelOrdinary = arcpy.sa.KrigingModelOrdinary("CIRCULAR", lagSize, majorRange,                                               partialSill, nugget)
RadiusVariable = 12
# Local variables:
Stations2010o3Data__2_ = "C:\\outVarRaster = ""
fields = ['O3_avg_dai', 'O3_avg_d_1', 'O3_avg_d_2', 'O3_avg_d_3', 'O3_avg_d_4', 'O3_avg_d_5', 'O3_avg_d_6', 'O3_avg_d_7', 'O3_avg_d_8', 'O3_avg_d_9', 'O3_avg__10', 'O3_avg__11', 'O3_avg__12', 'O3_avg__13',   'O3_avg__14', 'O3_avg__15', 'O3_avg__16', 'O3_avg__17', 'O3_avg__18', 'O3_avg__19', 'O3_avg__20', 'O3_avg__21', 'O3_avg__22', 'O3_avg__23', 'O3_avg__24', 'O3_avg__25', 'O3_avg__26', 'O3_avg__27', 'O3_avg__28', 'O3_avg__29', 'O3_avg__30', 'O3_avg__31', 'O3_avg__32', 'O3_avg__33', 'O3_avg__34', 'O3_avg__35', 'O3_avg__36', 'O3_avg__37', 'O3_avg__38', 'O3_avg__39', 'O3_avg__40', 'O3_avg__41', 'O3_avg__42', 'O3_avg__43', 'O3_avg__44', 'O3_avg__45', 'O3_avg__46', 'O3_avg__47']
for field in fields:
    counter += 1
    count_file = open("counter.txt", "w")
    count_file.write(str(counter))
    count_file.close()
# Process: Kriging
    tempEnvironment0 = arcpy.env.extent
    arcpy.env.extent = "56662.1788730418 308952.917606764 1698384.17317964 1912328.52737378"
    rasterKrig = "run{}_Krig_{}".format(counter, field)
    outKrig = arcpy.sa.Kriging(Stations2010o3Data__2_, field, kModelOrdinary, {2000}, RadiusVariable, outVarRaster)
    outKrig.save(rasterKrig)
    arcpy.env.extent = tempEnvironment0

The error message:

  Error:
  Runtime error
  Traceback (most recent call last):
    File "", line 54, in
    File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 2405, in Kriging
      out_variance_prediction_raster)
    File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py", line 47, in swapper
      result = wrapper(*args, **kwargs)
    File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 2397, in wrapper
      out_variance_prediction_raster)
    File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing_base.py", line 484, in
      return lambda *args: val(*gp_fixargs(args, True))
  RuntimeError: Object: Error in executing tool


texts: Is it possible to import ArcPy module to PyPy interpreter?

texts:
I'm trying to use PyPy for my code because of speed, but my code needs ArcPy module to run and I can't find out how to install or import this module to ArcPy. I'm using PyPy 2.3, Python 2.7 and Win 7.

Any suggestions?


texts: Trying to make xml data readable for Arcgis

texts: I am trying to convert xml data into tabular form for the process of geocoding.  We have the xml file with info that we are going to run through an address locator.  Do I need a script for this? Can we import directly into arc without coverting to xls, maybe using Python?


texts:
We do subsea surveys - seabed mapping etc. What software do you use to produce charts? I am especially interested in producing chart series.

It would be interesting to hear how you guys produce charts...


texts: Generate near table for x number of neighbours using arcpy

texts:
There are 3 likely scenerios that I am trying to capture near distances for:

An interchange subway station, which has 2 or more neighboring stations. That is, the station in question connects 2 or more major routes and has 2 or more neighboring stations.
A terminal subway station, which has only 1 neighboring station. This is the station at the end of the line.
An inline subway station, which has exactly 2 neighboring stations, one of either approach.

I am attempting to calculate a value one might call "average distance between neighboring stations"
The arcpy.GenerateNearTable_analysis() can handle two options: Distance to closest feature, and Distance between all features.
Does anyone have a clever method for solving for these scenarios? Note that each station is designated as "Interchange", "Terminal" or "Inline" in the attribute table under the field "StationType".
Added:
Here is some psuedo code based on @whuber's suggestion in the comments. I don't have time to figure this out just yet, so if anyone wants to take a stab at it you'll be rewarded with a checkmark! ;)
I have taken a look at the NetworkX library and it seems to work as I want it to.
Given the graph:
A —― B ―― C ―― D
     |
     E

as well as the nodes and links:
Nodes = ["A", "B", "C", "D", "E"]
Links = [("A", "B"), ("B", "C"), ("C", "D"), ("B", "E")]

def myFunction(node):
    identify the links that node belongs to
    count the number of links
    calculate the total link lengths
    divide the total link lengths by the number of links
    return someValue


texts:
I want to create a thematic map of my data. I downloaded a shapefile from Statistics Canada website displaying the various agricultural regions. I have data in Excel

How do I add this data and match it each polygon?

I've been reading the manual and seem to be getting more confused.

I am using Mapinfo Pro v8.5


texts:
I have georeferenced an old OS map and I have now imported the OS data shapefile for boundaries for England. How can I get rid of the extra data that extends beyond my map and just leave my scanned map with the boundary layer?

This is a screenshot of what it looks like.


texts: Shapefile to GeoJson conversion using C#

texts:
I have some shapefile which I want to convert to GeoJSON and store in MongoDB for querying.

Is it possible to take a shapefile and convert it to GeoJSON format in C#? I have been able to convert shapefiles to GeoJSON from the drivers available online and then import the documents in MongoDB, but I wish to get this all done in C# only. Reason being the ability to let me do some changes to the document on the fly.

If there is no solution available as such, then pointing me in the right direction would be helpful.


texts: How I can merge different objects of a theme by equal meaning of one attribute of a field (ArcMap)?

texts:
Could somebody advise, how I can combine (merge)  different objects of a shp theme with equal meaning of one attribute of a field (in ArcMap)?
I have shp theme with many objects, in the field "A" some of them have value “1”, some others “2”, others – “3”, etc
I need to combine all object with A=1 in a single object, all objects with A=2 – in a single as well, etc for all objects of the shp theme.

It seems very simple operation, but I can’t find this possibility in ArcMap. ((


texts: How to get data for multiple dates with OWSLib?

texts:
I am querying sentinel-2 data for multiple dates and want to get all results for these dates.

For example, I queried Arizona crater

from owslib.wms import WebMapService
import matplotlib.pyplot as plt
from pyproj import Proj, transform
import json

instance_id = ...

# in meters
region_size = (2000, 2000)

outProj = Proj(init='epsg:3857') # https://epsg.io/3857

inProj = Proj(init='epsg:4326') # https://epsg.io/4326

arizona_crater=(-111.02278,35.02722)

loc = arizona_crater

x, y = transform(inProj,outProj,loc[0],loc[1])

xupper = int(round(x - region_size[0] / 2))
xlower = int(round(x + region_size[0] / 2))
yupper = int(round(y - region_size[1] / 2))
ylower = int(round(y + region_size[1] / 2))

bbox = (xupper, yupper, xlower, ylower)
#bbox=(-11375363,3757392,-11372226,3759102)
size=region_size
#size=(1333, 716)

#wms_version = '1.1.1'
wms_version = '1.3.0'

wms = WebMapService('http://services.sentinel-hub.com/v1/wms/' + instance_id, version=wms_version)

print 'Type: ', wms.identification.type
print 'Title: ', wms.identification.title
print 'Contents: ', list(wms.contents)
print 'Operations: ',[op.name for op in wms.operations]
#print 'Capabilities: ', wms.getcapabilities()

info = wms.getfeatureinfo(
    layers=['1_NATURAL_COL0R'],
    srs='EPSG:3857',
    bbox = bbox,
    size = size,
    info_format='application/json',
    time='2015-01-01/2017-06-01/P1D',
    xy = (0, 0)
)
info = json.loads(info.read())

but getting only one result in feature list, although there are many

http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/#tiles/12/S/VD/

How to get all results for all dates or query possible dates with OWSLib?


texts:
I have a road network and a geometry for each road direction. One Geometry from A->B and another Geometry from B->A. This geometries are covering each other. Is there a tool available in QGIS to select a geometry based on the start and visualise the start of the geometrie? At the end I like to have a list of all geometries in one direction which I selected.

Did I have to write a own selection tool for this usecase?


texts: How to correlate raster map and elevation data?

texts: I have been provided with a set of raster maps. Thanks to another helpful user here I have been able to acquire some elevation data for the same region from gis.utah.gov in asc format. Now I need to correlate the two data sets. Can this be done in Matlab? If not, I am a proficient programmer and would appreciate instructions on how to do this in any language.


texts: Converting layout view graphics to features in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.4.

I have access to an .mxd file with a bunch of point locations placed on the map project as graphics. I'd like to convert these to features. I am aware of "Convert Graphics to Features" command under drawing (http://desktop.arcgis.com/en/arcmap/latest/map/page-layouts/converting-graphics-into-features.htm). However, no matter what I try (selecting, not selecting the elements; layout or data view), I am unable to get the command to be accessible (ie. not grayed out). I believe the problem is because the graphics were originally drawn on the layout view and the conversion seems to require they be on the data view.

Is there a way to transfer these graphics to the data view?


texts: Regression Model and Spatial Autocorrelation

texts:
I have used OLS and GWR for validating the dependencies between two distinct database. The Residual Square for GWR is 0.82 and thereby making it the right regression model to be used for determining the relationship between the two datasets.

What I wanted to know that, GWR being Local regression and OLS being global regression, which should be used where and when?

Also what does it really mean if Moran's I for the GWR model is Random?


texts:
I've been trying to clean up some vector data sets so that I can reproject them to WGS84 but I keep finding topological errors that fail the process. The one error I can't resolve involves a node on a polygon that shared by more than two (2) of that polygon's edges. For example:

Ideally, I would want to remove these types of topology errors so that I can reproject the shapefile but the QGIS 3 "Check Geometries" tool doesn't offer any suggestions on how to do it. I'm wondering what the best course of action is for these situations and if there is a tool/method that works best for these situations.


texts:
I'm currently trying to take a .tif file and polygonize it, I can do this if the output file format is a .shp but when I try to make a .gkpg file I get this error:

Traceback (most recent call last):
  File "C:\Program Files (x86)\GDAL\gdal_polygonize.py", line 169, in <module>
    dst_ds = drv.CreateDataSource( dst_filename )
AttributeError: 'NoneType' object has no attribute 'CreateDataSource'

I am unsure as to why this is occurring as http://www.gdal.org/ogr_formats.html says this is a possible file format.

As for the code I'm using:

call([sys.executable, 'C:\Program Files (x86)\GDAL\gdal_polygonize.py',infile.tif, '-f' ,'GKPG', outfile.gkpg])


texts: Creating Proj4 string for local coordinates

texts:
I am fairly new to coordinate system conversion, so am rather lost on the following problem.

We have been using Proj4 to convert between various coordinate systems which is fine for known coordinate systems such as MGA94 and WGS84 as the Proj4 strings are already defined.

The problem we're having is that we have some maps in local coordinates that we need to convert to a known coordinate system.

What we have

4 control points (X, Y, and Z) in Local coordinate system
4 control points (X, Y, and Z) in MGA94 (zone 51)
4 control points (X, Y, and Z) in WGS84

Each of the 4 control points for each coordinate system are the same 4 points on a map, just in the respective coordinate systems.

Is it possible using the known points of a defined coordinates system (such as WGS84 or MGA94 (Zone 51) to determine what the Proj4 string would be that represents the Local coordinate system to enable us to convert between them?


texts: Sub-watershed boundary is partial overlap the admistrative boundary

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

I would like to extract the sub-watershed for a particular area. The data I used was ASTER 30m DEM and a shapefile of my study site. The first step was to clip the DEM using the shapefile of my study site.Then, the general steps included Fil, Flow Direction, Flow Accumulation, Basin. Figure 1 shows the sub-watershed, and there are 15 sub-watersheds I extracted.

Here what my concern is that the sub-watershed boundary is usually partial overlap the admistrative boundary. And the purpose of extracting the sub-watershed is to calculate the water yield in each sub-watershed. The mode I used to calculate the water yield is the InVEST water yield model. The accuracy of sub-watershed boundary can affect the results of water yield simulation.

So, I created a large-scale DEM mosaic dataset (Figure 1b) and then extract the watershed boundary from this big DEM mosaic dataset. I found, for example, that sub-watered No.1 in Figure 1a just account for a very small part in the watershed No.2 in Figure 1b. Here I am not sure that if I can use the sub-watershed boundary shown in Figure 1a to calculate the water yield in each sub-watershed in my study site.

How can I get better results?


texts: ArcGIS Pro unable to export models as Python code

texts:
ArcGIS Pro 2.1 does not allow saving / exporting of a model constructed with ModelBuilder as a Python script.
It seems that the only way we can circumvent this limitation is to click on each individual geoprocessing history item, but, it is not very efficient to do so.
Is there another workaround?


texts:
PROJCS["3900_SC_Nad_83",GEOGCS
        ["GCS_North_American_1983",DATUM
            ["D_North_American_1983",SPHEROID ["GRS_1980",6378137,298.257222101]],
                PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]
            ],
                    PROJECTION["Lambert_Conformal_Conic"],
                        PARAMETER["False_Easting",1999996.000],
                        PARAMETER["False_Northing",0.000],
                        PARAMETER["Central_Meridian",-81.00000000],
                        PARAMETER["Standard_Parallel_1",32.50000000],
                        PARAMETER["Standard_Parallel_2",34.83333333],
                        PARAMETER["Latitude_Of_Origin",31.83333333],
                        UNIT["FOOT_US",0.3048006096]
        ]

I got this .prj file from a local county government in SC to update their parcel information and I am trying to import the shapefile into PostgreSQL and be sure that I have the correct SRID. This file indicates that it uses US feet, and I want to be sure that I am using the correct SRID does anyone know what the correct SRID is? 102733 which is in SC State Plane is supposed to be using the International Ft and not US Ft.


texts: How do I create a polygon layer from multiple touching line segments in QGIS?

texts:
I'm trying to create a base map for a Maine Coastal area, and wanted to create a polygon layer of the land cover from the coastline in QGIS.

Whether I created the coastline vector layer from the USGS national map or downloaded it from NOAA, both options leave me with a multipart line.  Ignoring islands the USGS data has ~100 almost connected line segments and the NOAA layer (below) is about 800 connected line segments.

I had thought I could disolve these and run linetopoly, but that seems to have some more processing needed that I don't understand, does anyone see what I should be doing to create or get my polygon layer?

The raw data looks like this, maybe I could make the coast into one segment?

856
 -69.340984,44.013133
...
 -69.345531,44.014013
END
857
 -69.356459,44.008586
...
 -69.356459,44.008586
END
858
 -69.340544,44.012986
 -69.341057,44.013719
etc...


texts: Excluding Specific Feature classes from ModelBuilder iteration?

texts:
I am creating a model the will iterate through the feature classes in the dataset.

It will select attributes from each feature class, the only issue is, two of the feature classes will cause a break in the process due to being void of the attribute.

Is there a way to exclude the two said feature classes from the iteration?


texts:
I would like to download some of the ESRI World Imagery using this however when I try to add the layer to an online map I get "The layer, World Imagery (for Export), cannot be added to the map". I've tried manually adding it directly to the map however this doesn't work either. I know questions about downloading ESRI World Imagery have been asked before but I have not seen any questions regarding using the ESRI World Imagery (for export) layer. I can't seem to find a way to use this service, is anyone familiar with this process?

In simpler terms, how can I export ESRI World Imagery using its World Imagery (for Export)?


texts: Adding values to landuse class raster using ArcMap

texts:
I am currently working on a map that holds my animals' projected radio-locations from my work this summer.
I have a current land class map downloaded from the USGS of the site; however, some of the pixels classifications wash over certain classes that are present within the field.
For instance, I am working with turtles that frequent a vernal pool and a watering hole. Within the map, there is no value for those landclass features in those locations present within the map. They are instead just a part of the larger surrounding landclass "i.e., deciduous forest, or scrub/shrub".
How do I delineate these areas and add them to my landuse class feature raster file?


texts: Calling %Variable% in Python Code Block in ModelBuilder?

texts: I have a parameter in a ModelBuilder model I am running that defines the geodatabase that all of the temporary files will be written to. I need to be able to call that variable in a field calculator that I run at one point. I used arcpy.da.TableToNumpyArray and need that geodatabase to be called in as the path to find the file / field. I have what I believe should be functioning code below but I keep getting a syntax error at line 3.


texts:
I'm trying to create a label that will show the following:
MP 23.6
The "23.6" is a double type column in my table, the "MP" is the text I want to add in front of it of the number (milepost).

I am new to QGIS.


texts: CORS issue with OpenLayers3 getGetFeatureInfoUrl

texts:
In an Angular directive I'm trying to create an Information tool with OpenLayers 3 using the getGetFeatureInfoUrl method. I've set up a function which is using $http to return the info as JSON:

 //GetFeatureInfo
        scope.getinfo = function() {
            var wmsSource = new ol.source.TileWMS({
              url: 'http://xxxxxx.com:8080/geoserver/yyyy/wms',
              params: {'LAYERS': 'zzzzz'},
              serverType: 'geoserver'
            });

            var wmsLayer = new ol.layer.Tile({
              source: wmsSource
            });

            map.on('singleclick', function(evt) {
              //document.getElementById('info').innerHTML = '';
              var viewResolution = /** @type {number} */ (view.getResolution());
              var url = wmsSource.getGetFeatureInfoUrl(
                  evt.coordinate, viewResolution, 'EPSG:3857',
                  {'INFO_FORMAT': 'application/json',
                    'propertyName': 'formal_en'});
              if (url) {
                  console.log('here');

                  var parser = new ol.format.GeoJSON();
                    $http({
                    url: url,
                    dataType: 'jsonp',
                    jsonpCallback: 'parseResponse'
                    }).then(function(response) {

                    var result = parser.readFeatures(response);
                    if (result.length) {
                        var info =[];
                        for (var i = 0, ii = result.length; i < ii; ++i) {
                            info.push(result[i].get('formal_en'));
                            }

                            console.log(info);

                         } else {
                             console.log('info else');
                    }
                  });
              } //
            });
        };

This is the error code from Console:
Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://xxxxxx.com:8080/geoserver/yyyy/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=zzzzz&LAYERS=zzzzz&INFO_FORMAT=application%2Fjsonp&propertyName=formal_en&I=34&J=135&WIDTH=256&HEIGHT=256&CRS=EPSG%3A3857&STYLES=&BBOX=-95393.41129989922%2C6826743.870205661%2C-94170.4188473364%2C6827966.862658224. (Reason: CORS header 'Access-Control-Allow-Origin' missing).

from what I understand I have to add a filter. I have tried adding this to end of C:\Program Files (x86)\GeoServer 2.9.0\webapps\geoserver\WEB-INF\web.xml

<filter>
    <filter-name>cross-origin</filter-name>
    <filter-class>org.eclipse.jetty.servlets.CrossOriginFilter</filter-class>
</filter>
<filter-mapping>
    <filter-name>cross-origin</filter-name>
    <url-pattern>/*</url-pattern>
</filter-mapping>

... and then restarting GeoServer but then GeoServer stopped working. Is this the correct way to fix the issue and if so am I editing web.xml correctly? Can anyone help please?


texts: ArcCatalog Registering a Feature Class into Geodatabase: Getting an Error

texts:
This is a SQL Server GeoDatabase. I have created a couple of tables using SSMS; then imported data to them. Then, inside ArcMAP, I created Layers and pointed their data to these tables via SDE. So now I have an MXD with Feature Classes.Then I defined Relationships between the Feature Classes.

So far so good. What I now need is to enable Attachments on a couple of these Feature Classes. I suppose I can't do that until I Register these Feature Classes into the GeoDatabase using ArcCatalog because the Enable Attachment option is not available?

Assuming that was the problem, I tried to do the 'Register with GeoDatabase' from ArcCatalog but am getting an error: 'Invalid Entity Type'? I have Googled but can't find anything which would help me. This is just a local test environment but I need to learn enough from it to implement an enterprise wide setup.

Any idea?


texts: Problem compiling QGIS against trunk - georeferencer

texts:
I got problems compiling QGIS (trunk) this morning. The georeferencer seems to be the source of the problem:

[ 81%] Building CXX object src/plugins/georeferencer/CMakeFiles/georefplugin.dir/qgsgeorefplugingui.cpp.o
/usr/local/src/Quantum-GIS/src/plugins/georeferencer/qgsgeorefplugingui.cpp: In member function ‘void QgsGeorefPluginGui::createActions()’:
/usr/local/src/Quantum-GIS/src/plugins/georeferencer/qgsgeorefplugingui.cpp:864:3: error: ‘mActionLocalHistogramStretch’ was not declared in this scope
/usr/local/src/Quantum-GIS/src/plugins/georeferencer/qgsgeorefplugingui.cpp:868:3: error: ‘mActionFullHistogramStretch’ was not declared in this scope
/usr/local/src/Quantum-GIS/src/plugins/georeferencer/qgsgeorefplugingui.cpp: In member function ‘void QgsGeorefPluginGui::addRaster(QString)’:
/usr/local/src/Quantum-GIS/src/plugins/georeferencer/qgsgeorefplugingui.cpp:1124:3: error: ‘mActionLocalHistogramStretch’ was not declared in this scope
/usr/local/src/Quantum-GIS/src/plugins/georeferencer/qgsgeorefplugingui.cpp:1125:3: error: ‘mActionFullHistogramStretch’ was not declared in this scope
make[2]: *** [src/plugins/georeferencer/CMakeFiles/georefplugin.dir/qgsgeorefplugingui.cpp.o] Error 1
make[1]: *** [src/plugins/georeferencer/CMakeFiles/georefplugin.dir/all] Error 2
make: *** [all] Error 2

Bit weird as I've been successfully compiled QGIS trunk on this machine running Ubuntu 12.04.

Any help would be most appreciated.


texts: Is it possible to perform a join on the fly within OpenLayers, to join a filtered tabular service to a polygon WFS layer?

texts:
We have around 2 million species observation records we want to display to users in an OpenLayers-based front-end.  I created a VIEW in our PostgreSQL database that summarizes these records by species, by township (an arbitrary set of 3,000 polygons covering our state, for those not familiar), and indicates the number of records and most recent record, per species, per township.  I would like to display these to end users, color-coded in a manner like the one shown in the image. The tabular data needed is simply a township ID, species ID, count of occurrences, and year of most recent occurrence.  The spatial data comprises the polygons representing townships.

Since the geometry for the township boundaries are the same in each case, I would prefer to have two materialized views that I serve via GeoServer and join together using OpenLayers:

1) A township view containing only the geometry and township ID, published as WFS

2) A tabular observation record view containing the township ID, species ID, and the summary info (count of records, most recent records).

The reason for going with a materialized view is display speed -- querying to summarize data in this way across 2,000,000 observations for 6,000 species takes a significant amount of time.  So, my plan is to materialize a view of this query, and refresh on a weekly basis.  But, since it is a materialized view, we'd essentially be storing a (large) duplicate set of township geometry for all 6,000 species.

Instead, I'd hoped that I could have a materialized view of species observation summary data that I could join to the materialized view of township geometry, so that we're not storing or conveying duplicate geometry.

I saw this question, but it doesn't quite answer my question, and I thought it was possible that something has changed in OpenLayers since then that might make this possible.

Is this something that is possible to do, without a significant amount of hacking?


texts:
I created a custom toolbox in ArcGIS to add new Python scripts to be used by multiple persons.
When I save the tool box and the add it to their ArcToolbox they can edit the source code and the parameters.
How can I hide it or make it inaccessible?


texts: Merge Management Not Working with more than two input

texts:
I have code to merge two featureclasses with Merge_Management geoprocessing tool using ArcObject.

Below is the code :

 for (int i = 0; i < map.LayerCount; i++)
            {
                if (map.Layer[i].Name == "HT Overhead Line")
                {   ILayer HTOverheadLine = map.Layer[i];
                    IFeatureLayer fl1 = new FeatureLayerClass();
                    fl1 = HTOverheadLine as IFeatureLayer;
                    IDataset dataset = (ESRI.ArcGIS.Geodatabase.IDataset)(HTOverheadLine);
                    STR = dataset.Workspace.PathName + "\\" + ((IDataset)fl1.FeatureClass).BrowseName;
                }
                else
                {
                    if (map.Layer[i].Name == "HT Underground Line")
                    {
                        ILayer HTUndergroundLine = map.Layer[i];
                        IFeatureLayer fl2 = new FeatureLayerClass();
                        fl2 = HTUndergroundLine as IFeatureLayer;
                        IDataset dataset = (ESRI.ArcGIS.Geodatabase.IDataset)(HTUndergroundLine);
                        STR = STR + ";" + dataset.Workspace.PathName + "\\" + ((IDataset)fl2.FeatureClass).BrowseName;
                    }

                }

            }
            IGeoProcessor gp = new GeoProcessorClass();
            gp.AddOutputsToMap = true;
            gp.OverwriteOutput = true;
            IVariantArray parameters = new VarArrayClass();

            parameters.Add(STR);
            //MessageBox.Show(STR);
            parameters.Add(@"outpath\Merge.shp");
            gp.Execute("Merge_management", parameters, null);

This is working fine and merges the two featureclasses, however when I am trying to merge three featureclasses, it gives a strange error with below code

 for (int i = 0; i < map.LayerCount; i++)
        {
            if (map.Layer[i].Name == "HT Overhead Line")
            {   ILayer HTOverheadLine = map.Layer[i];
                IFeatureLayer fl1 = new FeatureLayerClass();
                fl1 = HTOverheadLine as IFeatureLayer;
                IDataset dataset = (ESRI.ArcGIS.Geodatabase.IDataset)(HTOverheadLine);
                STR = dataset.Workspace.PathName + "\\" + ((IDataset)fl1.FeatureClass).BrowseName;
            }
            else
            {
                if (map.Layer[i].Name == "HT Underground Line")
                {
                    ILayer HTUndergroundLine = map.Layer[i];
                    IFeatureLayer fl2 = new FeatureLayerClass();
                    fl2 = HTUndergroundLine as IFeatureLayer;
                    IDataset dataset = (ESRI.ArcGIS.Geodatabase.IDataset)(HTUndergroundLine);
                    STR = STR + ";" + dataset.Workspace.PathName + "\\" + ((IDataset)fl2.FeatureClass).BrowseName;
                }
                else
                      {
                          if (map.Layer[i].Name == "LT Underground Line")
                          {
                              ILayer LTUndergroundLine = map.Layer[i];
                              IFeatureLayer fl3 = new FeatureLayerClass();
                              fl3 = LTUndergroundLine as IFeatureLayer;

                              IDataset dataset = (ESRI.ArcGIS.Geodatabase.IDataset)(LTUndergroundLine);
                              STR = STR + ";" + dataset.Workspace.PathName + "\\" + ((IDataset)fl3.FeatureClass).BrowseName;

                          }

            }

        }
        IGeoProcessor gp = new GeoProcessorClass();
        gp.AddOutputsToMap = true;
        gp.OverwriteOutput = true;
        IVariantArray parameters = new VarArrayClass();

        parameters.Add(STR);
        //MessageBox.Show(STR);
        parameters.Add(@"outpath\Merge.shp");
        gp.Execute("Merge_management", parameters, null);

I do not understand whats wrong between these two codes?


texts: Using a file .shp with ESRI Geo Tools for Hadoop

texts:
I've a file .shp that I would like to load into a TABLE using HIVE with ESRI Geo Tool for Hadoop.

These tools read JSON files and I've tried to convert my file into a .json format using GDAL Library with command ogr2ogr -f geoJSON com.json com.shp

This format isn't suitable for ESRI GeoTools for Hadoop that need of unenclosed JSON. How can I correctly convert my shapefile?