texts: leaflet choropleth map using dynamic value

texts:
I am doing a choropleth map using leaflet and JSON data. I got a clean tutorial in leaflet site. In that tutorial the used us-states.js JSON file where the choropleth was based on the feature.properties.density feature of the JSON if I am not wrong. If we see the JS file we can found there is a field called "ID" which is state id. I have a database where I have 5 different categories of data based on state ID. Categories are population, density, male_population, female_population, literacy_rate. I am calling those data through ajax and get a GeoJSON data state-wise as follows (Its a dummy data).

[{"state_id":"01","population":"123456","density":"1234","male_pop":"65432","female_pop":"57421","literacy_rate":"98"}]
...
[{"state_id":"50","population":"123456","density":"1234","male_pop":"65432","female_pop":"57421","literacy_rate":"98"}]

I want to integrate this data as choropleth value. when I call on literacy function the variation of color will be based on literacy_rate data. I can make the changes in getcolor function category wise.

Thanks


texts: QGIS Zonal Statistics: Sum of Partial Area Larger Than Sum of Total Area

texts:
I've been trying to use zonal statistics to get the proportion of a population in an area relative to the larger administrative district that the area falls under. My workflow is:

Run zonal statistics(sum) for red polygon.
Run zonal statistics(sum) for blue polygon.
Divide output of #2 by output of #1.

Strangely, the proportion comes out as greater than 100%, i.e. output of #2 is larger than #1. Why would this be happening and how can I solve this? I thought the issue was being caused by invalid geometries, but that doesn't seem to be the issue.


texts: Which buildings have access to nature where more than 50% of the nature area is closer than 3km?

texts:
I have 1) a file with all the residential buildings in a small municipality and 2) a file with all the surrounding green areas.
I'm having a hard time solving this because the residential buildings are so many that doing a buffer for individual buildings produces way too many overlapping circles. In spite of this, I tried following this tutorial Percentage of polygon in one shapefile within polygon of another . I used the overlap analysis and set the green areas as input layer, and the 3km buffer around residential buildings as the overlay layer. But in the attribute table the analysis returns a column with the area percentage where all the values are 99 or 100. Which makes me think that the analysis isn't performed on each individual buffered circle, but on the whole shape instead. I tried the other way around as well (buffer layer as input and green areas as overlay) but the percentage values it returns are very small and none larger than 49%, which can't be right either.
I have also tried performing a 3km buffer on the green areas and see which buildings it intersects with but that doesn't give the full picture since there are other buildings that have access to less than 100% of the green area within 3km.
Here is a screenshot of my layers.
Anyone have any tips on how I should actually go about this?


texts:
I made a classification mask of the layer and get the following figure:

I got several vertex pixels. A vertex pixel is a pixel that attaches to a bigger group by just one vertex (shown in red circles).
Is there a way to remove these pixels? I have tried the connectedPixelCount method but it looks like it will group the entire thing into one instead of masking the vertex pixel.


texts: Using Row Iterator and IDW in ModelBuilder?

texts:
I am using ModelBuilder to perform inverse distance weighting. I would like to interpolate values using a point file with values for several unique dates - resulting in a raster dataset for each unique day. My data is in long format where each row has a value for each coordinate/date.

I have tried to create a model that performs the IDW on each unique day, using row iterator to iterate through each "date". It runs, but I don't think it's correct. Here is the model.

This results in a raster dataset for each unique day. The problem is, each raster is identical. They should not be - the interpolated values should vary each day. Maybe I am using the wrong iterator.

Can you please suggest ways to get this to work?


texts: How to reuse memory pointer of GDAL Memory driver

texts:
I am doing the same as in [GDAL API]: can't save image in some formats

That is, saving a JPEG image with GDAL. I know I'd have to create a MEM buffer and then CreateCopy to JPEG. However, I'd like to avoid the allocation of two buffers since my image might be big and I use my own memory buffer.

Is there a way to pass the data memory address that I'd like to use to the MEM driver?

I've seen this mysterious article http://www.gdal.org/frmt_mem.html which references the DATAPOINTER option. I've seen another reference here
which does

sprintf(filename,"MEM:::DATAPOINTER=%d,PIXELS=%d,LINES=%d,BANDS=1,DATATYPE=%d,PIXELOFFSET=0,LINEOFFSET=0,BANDOFFSET=0",datpt,pixels,lines,datatype);
dataset = GDALOpen(filename,GA_Update);

How would this translate to c++? Maybe

    papszOptions = CSLSetNameValue( papszOptions, "DATAPOINTER", datapointer );

and so on?


texts: Combine leaflet-tilelayer-mask and Leaflet.TileLayer.ColorFilter plugins

texts:
I am trying to color my tileLayer mask differently and I found the colorFilter plugin. I can't seem to figure how to use the color filter tile created in my mask tileLayer.
let map = L.map("map", {
  zoom: 16,
  center: [35.9734, 117.7926]
});

let sampleFilter = [
  'bright:85%',
  'hue:83deg',
  'saturate:97%'
];

// background layer -- color filter plugin
let myTileLayer = L.tileLayer.colorFilter('https://{s}.tile.thunderforest.com/pioneer/{z}/{x}/{y}.png?apikey=xxxxxxxxxxxxxxxx', {
  attribution: '',
  filter: sampleFilter,
}).addTo(map);

// mask layer -- tileLayer mask plugin
let mask = L.tileLayer.mask('https://{s}.tile.thunderforest.com/pioneer/{z}/{x}/{y}.png?apikey=xxxxxxxxxxxxxxxx', {
   maskSize: 512,
   maxZoom: 18,
   maxNativeZoom: 17
 }).addTo(map);

mask.setCenter(e.containerPoint);


texts: Copy multiple Feature Classes using ONE CopyFeatures command

texts:
My script below works fine, however, how can I combine each copy.features task into one?
Can I use a tool where I can simply list all the 14 feature classes I want to export in ONE LINE? Then run arcpy.CopyFeatures_management once underneath.

print "Script started..."
print "***"
print "This monthly script will back up all grounds maintenance feature classes as shapefiles."
print "***"
print "The backup folder is:"
print "K:\GIS2016\GIS_Data\Curo_data\Estates_data\Grounds_Maintenance\Grounds_Maintenance_Backups"
import arcpy
import os
import datetime
from arcpy import env

Cur_Date = datetime.datetime.now().strftime("%d%B%Y")
print Cur_Date

print "***"
print "Setting script to overwrite previous files..."
# Overwrite pre-existing files
arcpy.env.overwriteOutput = True

newpath = r'\\somernt\curo\GIS\GIS2016\GIS_Data\Curo_data\Estates_data\Grounds_Maintenance\Grounds_Maintenance_Backups\GM_' + Cur_Date
if not os.path.exists(newpath): os.makedirs(newpath)

print "***"
print "Starting backup process..."
env.workspace = r"\\somernt\curo\GIS\GIS2016\GIS_Data\Curo_data\Estates_data\Grounds_Maintenance\Grounds_Maintenance.gdb"
arcpy.CopyFeatures_management("LITTER_BINS", newpath + "\\LITTER_BINS.shp")
print "Bins"
arcpy.CopyFeatures_management("AMENITY_GRASS", newpath + "\\AMENITY_GRASS.shp")
print "Amenity grass"
arcpy.CopyFeatures_management("AMENITY_GRASS_SHELTERED", newpath + "\\AMENITY_GRASS_SHELTERED.shp")
print "Sheltered grass"
arcpy.CopyFeatures_management("HEDGES_1SIDE_ONLY", newpath + "\\HEDGES_1SIDE_ONLY.shp")
print "Hedges- 1 side only"
arcpy.CopyFeatures_management("HEDGES_1SIDE_TOP", newpath + "\\HEDGES_1SIDE_TOP.shp")
print "Hedges- 1 side and top"
arcpy.CopyFeatures_management("HEDGES_2SIDES_TOP", newpath + "\\HEDGES_2SIDES_TOP.shp")
print "Hedges- 2 side and top"
arcpy.CopyFeatures_management("LEAF_CLEARANCE", newpath + "\\LEAF_CLEARANCE.shp")
print "Leaf clearance"
arcpy.CopyFeatures_management("MEADOW_GRASS", newpath + "\\MEADOW_GRASS.shp")
print "Meadow grass"
arcpy.CopyFeatures_management("ROUGH_GROUND", newpath + "\\ROUGH_GROUND.shp")
print "Rough ground"
arcpy.CopyFeatures_management("SHELTERED_GRITTING_AREAS", newpath + "\\SHELTERED_GRITTING_AREAS.shp")
print "Sheltered gritting areas"
arcpy.CopyFeatures_management("SHRUBS_ORNAMENTAL", newpath + "\\SHRUBS_ORNAMENTAL.shp")
print "Shrubs"
arcpy.CopyFeatures_management("SHRUBS_WALL", newpath + "\\SHRUBS_WALL.shp")
print "Shrub wall"
arcpy.CopyFeatures_management("SOMER_GRIT_BINS", newpath + "\\SOMER_GRIT_BINS.shp")
print "Grit bins"
arcpy.CopyFeatures_management("SWEEPING", newpath + "\\SWEEPING.shp")
print "Sweeping (hard surfaces)"

print "Script finished"


texts: Ogr2Ogr: Failed to create GPX file output.gpx

texts:
When i try to export an shapefile to gpx using ogr2ogr.

ogr2ogr -f GPX output.gpx D:\mypaste\input.shp -sql "SELECT Id AS name FROM layername"

It gives me the following error "Failed to create GPX file output.gpx." "GPX driver failed to create output.gpx." Does anyone knows why this happens? How can i rectify this?


texts: Is there a way to automate PDF reports in Grass GIS through the web?

texts:
I am an experienced programmer, but new to GIS. My question: I would like to automate the creation of PDF reports from a website, using Grass GIS.

I've come across Web Processing Service (WPS) (http://www.opengeospatial.org/standards/wps/) but am unsure how this fits into any possible report automation.

Is WPS the way to automate report generation, or is some kind of direct scripting (such as with Python) a better solution?


texts: Calculating length of just portion of linestring that intersect with polygon using PostGIS?

texts:
I want to calculate the portion of the linestring that lies within the limits of a polygon. Better with an image

What is the length of the portion of linestring within Province A? and within Province B? I'm looking for a query that returns something like:

province_name | length   |
--------------------------
Province A    | 30.5342  |
Province B    | 19.4321  |

My approach is a PL/pgSQL script that loops over all the points in the linestring and set a mark in the first point that is in province B (point N). Then, calculate the length from point 0 to point N (distance within province A) and the length from point N+1 to last point (distance within province B). But sounds rusty, and difficult to extend to an indefinite number of provinces.

Is there any better approach? Suggestions welcome. I'm using PostGIS 2.0, by the way.

@HeyOverThere suggests a query that almost works, with a little modification:

SELECT a.name, sum(ST_LENGTH(ST_Intersection(a.geom, b.geom))/1000)
FROM polygons a, lines b
WHERE ST_Intersects(a.geom, b.geom) group by a.name order by a.name

That query returns 2 columns, with the province name and a length. But the length is exactly twice the value it should be. Does it make sense?

More info: the provinces are really multipolygons (but all share a the province name as common field), and the line has been built from gps tracks (points) and st_makeline.

As @GetSpatial suggested, I've tried using subqueries to sum and group the results of the first query. Same result. The length of each intersection is correct. So, maybe the intersection is two times as long as I expect (?).

Here, the 3 things I've tried. All of them get a distance 2x the expected one (link to real data at the end):

First try, with one query:

SELECT a.province, sum(ST_LENGTH(ST_Intersection(a.geom, b.the_geom))/1000) as length_per_province
FROM provinces a, line_tracks b
WHERE ST_Intersects(a.geom, b.the_geom) group by a.province order by a.province;

Second try, with CTE:

with line_segments as (
SELECT a.province, ST_LENGTH(ST_Intersection(a.geom, b.the_geom))/1000 as length
FROM provinces a, line_tracks b
WHERE ST_Intersects(a.geom, b.the_geom) order by a.province
) select province, sum(length) as length_per_province from line_segments group by province;

Third one, with subquery, as suggested:

select province, sum(length) as length_per_province from(
SELECT a.province, ST_LENGTH(ST_Intersection(a.geom, b.the_geom))/1000 as length
FROM provinces a, line_tracks b
WHERE ST_Intersects(a.geom, b.the_geom)) as foo group by province;

I get this, in all cases:

"province";"total_length"
"p1";62.8687677452213
"p2";40.254207910339

The total length of the line:

select sum(st_length(the_geom))/1000 as total_length from line_tracks

The  result

"total_length"
51.5614878277802

And the 2 tables here (SET client_encoding = 'LATIN1'): https://dl.dropboxusercontent.com/u/6599273/gis_data/data.sql.tar.bz2


texts:
I have a buffers and I want to select all points within each of the buffer and output it as a table. That means one buffer may have 0 to about 10 points in it.
I have tried spatial Join and intersect in ArcMap 10.1. Its returning where there is zero points correctly. In instances with more than one, it lists only one randomly instead of all of them
What other tool can I use?


texts: How to use a .kml file to create a 3D mesh from contours?

texts:
I have downloaded a .kml file from this source: http://www.carpworld.co.za/index.php/dam-contours/boskop-dam.
This is how it looks in Google Earth:

This file contains contours of Boskopdam (a local lake). The contours don't have height data/attributes, or non that I am aware of, but the task I'm trying to complete is converting or displaying these contours in a 3D fashion like this:

I've seen a few videos on how to use a .dem file in Blender to create a 3D mesh, but thats about all I know.

I have QGIS and have little to no knowledge on using it.


texts: QGIS tool for aggregating population data in fishnet

texts:
I'm looking for a tool or method to aggregate a vector fishnet based on a minimum value. Specifically, I have a fishnet containing a population value and I would like to join together cells that contain fewer people than X in an intelligent way.

Is there a straight forward way to do this in QGIS? Preferably without too much command line, but I'll get my hands dirty if that's what it takes.


texts: QGIS 2.14 WMS. Can I add WMS services from this site to QGIS?

texts:
I'm interesting by this http://mapy.geoportal.gov.pl/imap/?gpmap=gp3&locale=en site content and I need to add it to QGIS.
For example, Ortofoto wms

When I'm trying to add this WMS to QGIS I catch error

  Failed to download capabilities:

  Download of capabilities failed: Operation canceled My settings

My QGIS is 2.14.7


texts: Is there any alternative for reshapeGeometry in Qgis 2.4 version

texts:
I am splitting the polygon feature if it intersects with the line feature.
I am using reshapeGeometry() but it is nt wrking in qgis2.4. Is  there any alternative to use for this one.it is wrking perfectly fine in qgis 2.2.

for feature in polygonpr.getFeatures():
    geometry = QgsGeometry.fromPolygon(feature.geometry().asPolygon())
    for line in linepr.getFeatures():
        t = feature.geometry().reshapeGeometry(line.geometry().asPolyline())
        print t
Traceback (most recent call last):
  File "<input>", line 4, in <module>
Exception: unknown

I am getting the following error.

Is there any way to replace reshapeGeometry? Any solution


texts:
What is wrong with my code?
Template = {
    render : function(){

        if (!Map.map){
            Map.create_map();
            Map.set_tile_layer();
        }

        API.get_businesses();
    }
}
API = {
    get_businesses : function(){
        $.ajax({
            url : '/api/mapviz_get_businesses',
            type: 'GET',
            datatype: 'json',
            data: {},
            success: function(d){
                $(d).each(function(i, item){
                    var marker = L.marker([item.coords.coordinates[0], item.coords.coordinates[1]]).addTo(Map.map)
                })
            }
        })
    }
}
Map = {
    map : null,
    create_map : function(){
        Map.map = L.map('map').setView(
            [32.7766642,-96.79698789999998],5);
    },
    set_tile_layer : function(){
        L.tileLayer('https://api.mapbox.com/styles/v1/{id}/tiles/{z}/{x}/{y}?access_token={accessToken}', {
            attribution: 'Map data &copy; OpenStreetMap contributors, Imagery © Mapbox',
            maxZoom: 18,
            id: 'mikejohnsonjr/ckmbpaffj6dmw17s3h4ll5wje',
            tileSize: 512,
            zoomOffset: -1,
            accessToken: 'pk.eyJ1IjoibWlrZWpvaG5zb25qciIsImEiOiJjaWY5dXJ4bnAxaXVnc3VrbnpocXd3MXF5In0.iVbz_2py7DdmvpoQTzVe0A'
        }).addTo(Map.map);
    }
}

Template.render()

The map loads but the markers do not, not even when I open up console and do the L.marker() function manually.
Can anyone help?
I'm following this very simple guide and cannot see where I'm going wrong: https://leafletjs.com/examples/quick-start/
Everything returned by the ajax function is fine.


texts: Calculating kernel density for polygons in ArcGIS Desktop

texts:
I have a population data set at the census block level. My goal is to use the data set to identify some population centers in the region based on population density. The simplest way to do it is to calculate the density by census blocks and use an arbitrary density cut-off value to determine those centers. One shortcoming of this method, however, is that it relies on census blocks, which have irregular size and shape.

The other approach I am thinking of is to covert census blocks into points based on their centroids, and calculate kernel density based on these points. But I am not sure if this makes sense.

In other words, does this approach help me overcome the problems that are associated with the first approach?

Also, if it makes sense to do that, how the output raster size should be determined?


texts: MOLUSCE plugin issue- getting same LULC raster for all simulation

texts:
I'm using MOLUSCE plugin in QGIS to predict Land Use Land Cover (LULC) for future.
LULC data for 2005 and 2010 is given as initial and final raster along with spatial variables DEM, distance from road, distance from stream and slope raster files. LULC for 2015 is simulated with number of simulation iteration 1 and validated. When I'm trying to simulate for future years by changing the number of simulation iterations to 2, 3, 4, ... I'm getting the same LULC raster file simulated for 2015.


texts: PostGIS comes with a table of spatial reference systems with each one being identitfied by a SRID. However, from a user's perspective, SRID is just a number and doesn't tell them much. Is there a resource that associates SRID with  more descriptive names?


texts: Unable to open broken data set in MapInfo 17.0

texts:
When I try to open an workspace/map in MapInfo 17.0 that has a broken dataset or it can't find the direction to an dataset, my MapInfo freezes.
Is there a way I can skip past the broken dataset and still open the workspace/map without MapInfo freezing?


texts: Adjusting extents so that same-CRS shapefiles overlap in QGIS?

texts:
I have two files with identical projections/CRS. Using both OTF on & off I cannot get them to overlap: one shapefile contains coastal polygons of Massachusetts, and the other shapefile contains coastal lines for MA.

 +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs

 +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs

however the extents are

 xMin,yMin 2321544.64,258328.06 : xMax,yMax 2333464.74,277553.12

vs.

 xMin,yMin -176.685,17.9108 : xMax,yMax -65.2246,71.3413

I've followed the steps in http://qgis.spatialthoughts.com/2012/04/tutorial-working-with-projections-in.html , tried with OTF turned on & off: the layers still do not overlap.

EDIT:

I've saved the MA coast line SHP ( Layer > save as ) using each of the Layer CRS, Project CRS and Selected CRS options: the extent remains the same. As far as I can tell, the chosen projection (EPSG:2163) is a PCS, not a GCS.


texts: Optimizing Google Earth Engine cloud cover for Sentinel 5p

texts:
I am trying to generate cloud free Sentinel 5p images, I have a basic cloud detection code working but it runs very slowly (most likely due to nested map and reduce features).

Is there a way to do this better?

var coll,coll1,median,median1;
button.onClick(function(){
  chart_panel.clear();
  Map.clear();
  first();
  second();
});

var first=function(){
  coll=p_var3.filter(ee.Filter.calendarRange(start_month_var,start_month_var,'month') ).filter(ee.Filter.calendarRange(start_year_var,start_year_var,'year'))
  .filterBounds(city_var1).select(p_var2);
  chartfunction1();
  median = coll.median().clip(city_var1);
  var Extrema=median.reduceRegion(ee.Reducer.minMax(),city_var1,250);
  var vis_params=ee.Dictionary({
    min: Extrema.getNumber(p_var4),
    max: Extrema.getNumber(p_var5),
    opacity: 1.0,
    palette: ["blue","red","yellow","green","cyan","black","purple"],
    bands:[p_var2]
  });
  Map.addLayer(median,vis_params.getInfo(), 'true-color composite');
};

var second=function(){
  coll1=p_var3.filter(ee.Filter.calendarRange(end_month_var,end_month_var,'month') ).filter(ee.Filter.calendarRange(end_year_var,end_year_var,'year'))
  .filterBounds(city_var1).select(p_var2);
  chartfunction2();
   median1 = coll1.median().clip(city_var1);
  var Extrema1=median1.reduceRegion(ee.Reducer.minMax(),city_var1,250);

  var vis_params1=ee.Dictionary({
    min:Extrema1.getNumber(p_var4),
    max:Extrema1.getNumber(p_var5),
    opacity: 1.0,
    palette:["blue","red","yellow","green","cyan","black","purple"],
    bands:[p_var2]
   });
   Map.addLayer(median1,vis_params1.getInfo(), 'true-color composite2');
};

var chartfunction1=function(){
  var chart = ui.Chart.image.seriesByRegion({
    imageCollection: coll.select(p_var2),
    regions: city_var1,
    reducer: ee.Reducer.median(),
    band:[p_var2],
    scale: 30,
    xProperty:'system:time_start',
    seriesProperty: 'label'
  }).setChartType('ScatterChart')
  .setOptions({
    title: 'Pollution',
    vAxis: {title: p_var1},
    hAxis: {title: 'Month'},
    pointSize:3,
    series: {
              0: {color: 'green'}
            }
  });
  chart_panel.add(chart);
};
var chartfunction2=function(){
  var chart1 = ui.Chart.image.seriesByRegion({
    imageCollection: coll1.select(p_var2),
    regions: city_var1,
    reducer: ee.Reducer.median(),
    band:[p_var2],
    scale: 30,
    xProperty:'system:time_start',
    seriesProperty: 'label'
  }).setChartType('ScatterChart')
  .setOptions({
    title: 'Pollution',
    vAxis: {title: p_var1},
    hAxis: {title: 'Month'},
    pointSize:3,
    series: {
              0: {color: 'blue'}
            }
  });
  chart_panel.add(chart1);
};


texts: Understanding NOAA/GFS0P25 dataset of Google Earth Engine

texts:
I'm trying to get time-series data from NOAA/GFS0P25 data set on GEE.
My code:

var m2 = "NOAA/GFS0P25"

var band1 = 'total_precipitation_surface'
var data_collection = ee.ImageCollection(m2);
var dataset = data_collection.filterBounds(plots)
                    .filterDate('2019-03-04', '2019-03-05')
                    .filterMetadata('forecast_hours', 'greater_than', 0)
                    .select(band1);

print (dataset)

where plots is just a point with x,y (8.228528757789945,49.390447480404966).
I was expecting to get only 4 elements (1 days * 4 times)  accrotding to the docs: The 384-hour forecasts, with 3-hour forecast interval, are made at 6-hour temporal resolution .
However, I received 688 elements, which makes it impossible to obtain a larger time series than one day due to the 5000 elements limit.
A subset of the elements looks like this:

I don't understand why are there 688 elements and not 4. I also plotted the data and I got only 4 points of data as expected:

Why is this discrepancy? how can I eventually get a time series similar to the plot but not for just 1 day but for 6 months?


texts: Openlayers 4: Add OSM labels without geometries

texts:
Is there a way to add a layer in Openlayers 4 containing only OpenStreetMap labels (for cities and roads mainly).?
I know how to display OSM layers with and without labels using these links:

https://tiles.wmflabs.org/osm/${z}/${x}/${y}.png
https://tiles.wmflabs.org/osm-no-labels/${z}/${x}/${y}.png

I would like to display a layer containing only labels and not geometries over a satellite image.

Any idea on how I can achieve that?


texts: Adding the image file to KML file using simplekml

texts:
I have made simple tool to convert CSV file containing latitude and logitude to KML　file.
Adding the image using the relative path is can be done easily for each location using following code:
pnt.style.balloonstyle.text = "<![CDATA[

]]>"]
There is any way to save image permanently to the KML file not in the relative path,so that next time KML file will not loaded from file relative path but from own file system.


texts: Add WMTS to Leaflet using its GetCapabilities URL

texts:
I am using the Leaflet.TileLayer.WMTS plugin (https://github.com/mylen/leaflet.TileLayer.WMTS). However, I believe that this only works with WMTSes which add parameters to GetTile requests as querystring variables in key=val pairs (eg http://www.provider.net/path?var1=val1&var2=val2 etc). I need to add a WMTS which uses a URL path for its parameters (eg http://www.provider.net/path/val1/val2 etc). The WMTS in question is https://www.basemap.at/wmts/1.0.0/WMTSCapabilities.xml.

The way other clients handle this is by parsing the GetCapabilities URL, which provides the parameter mapping for GetTile requests. I see that OpenLayers has a method for parsing GetCapabilities to add the WMTS (https://openlayers.org/en/latest/examples/wmts-layer-from-capabilities.html).

I'm aware that at least one other question has this as the underlying issue (https://stackoverflow.com/questions/44300140/wmts-layer-syntax-for-leaflet), but it doesn't make it clear what the underlying problem is, so I thought it was worth asking a new question.

I know that I can either read the GetCapabilities XML myself (or look at network requests from a working client) and determine the structure to add a plain XYZ tileserver. However, I am trying to code this into qgis2web, which needs to accomplish this programmatically.

So, is there any way of adding a WMTS to a Leaflet webmap if you only have the URL for its GetCapabilities method?


texts: Cannot read X Y dimension from OCO-2 Satellite Dataset NetCDF4

texts:
I have downloaded some data from the OCO-2 satellite and what to show it spatially in ArcGIS (10.5). I have done this with the AIRS and AQUA satellites of the same file types and there is not issue. The file extension is ".nc" and the file type that I downloaded was "NetCDF 4".

Below is screenshot of the "NetCDF to Feature Class" tool. When I select the variable that I want, there is nothing to select in the lat and long fields.
The lat and long variables are available in the "Variables" drop down, but not in the X and Y.

This guy is having the same problem as me, but his question is unanswered:
https://geonet.esri.com/thread/168798

Here is where I downloaded the file from:
https://mirador.gsfc.nasa.gov/cgi-bin/mirador/granlist.pl?page=1&dataSet=OCO2_L2_Standard&version=7&allversion=7&version=7r&allversion=7r&keyword=OCO&pointLocation=(-90,-180),(90,180)&location=(-90,-180),(90,180)&searchType=Location&event=&startTime=2009-03-31&endTime=2017-01-27%2023:59:59&search=&CGISESSID=7869eb21869b19ce4393bb1280781991&nr=21361&temporalres=16%20days&prodpg=http://disc.gsfc.nasa.gov/datacollection/OCO2_L2_Standard_V7r.html&longname=OCO-2%20Level%202%20geolocated%20XCO2%20retrievals%20results,%20physical%20model,%20Retrospective%20Processing%20V7r&granulePresentation=ungrouped


texts: Splitting self-intersecting line at only one point using QGIS?

texts:
I have a self-intersecting line that I tried to split using the Split Features tool in the Advanced Digitizing Toolbar. The problem is that when I try to split the line at one point, the line splits at every point where it crosses over itself.

This line was originally a tracklog generated in WGS 84 using a Garmin GPS. I imported it into QGIS in GPX format, and saved as a shapefile.

The line is styled with a red marker at the first node.  Before splitting at the blue circle:

After splitting:

What's going wrong here?


texts: Is there a way to detect that the extent has changed in an app using the ESRI iOS SDK?

texts: Is there a delegate, callback, or other method that can be used to detect that the extent of the map has been changed in the ESRI iOS SDK?


texts: ArcGIS Linear Referencing successive overlay operations

texts:
I have 8 tables with RouteID, BeginMeasure and EndMeasure.  I want to successively overlay them starting with the first two and then overlaying the result and the next table then the result of that overlay and the next table until all tables have been run.  I am doing UNION as overlay type.  For some reason when I get to the 4th overlay the tools stops overlaying and just seems to start appending.

If two records have the same measure range and all other values are ok the operation doesn't seem to combine the records.  Are Null vs blank strings considered conflicts for unioning records?

EDIT: The steps are are to use the Result of the first overlay and as the input table and the next table to be overlayed as the overlay table (second parameter).  Is there an ArcMap memory issue?  What is the difference between the input table and the overlay table other than just attribute ordering?

EDIT 2: Looks like there may be memory issues in Arcmap.  After testing/executing 6 or 7 overlays with different datasets it seemed to break down and not want to work.  Closed ArcMap and re opened and started working again.


texts:
controlPoint = function (feature) {
                        alert(feature);
                        pointNumber = 1;
                        draw.deactivate();
                        noktasil(feature);

                    };
                     draw = new OpenLayers.Control.DrawFeature(
                        wfs_layer, OpenLayers.Handler.Point,
                        {
                            title: "Draw Feature",
                            displayClass: "olControlDrawFeaturePolygon",

                            multi: true,

                        }
                        );
                     draw.events.register("featureadded", ' ', controlPoint);

That's my codes  for creating new feature, it works. But i want to delete feature without click. And as you can see,  i store feature which i create and send it for deleting.

function noktasil(feature) {
                    // if feature doesn't have a fid, destroy it
                    var sil = new DeleteFeature(wfs_layer, { title: "Delete Feature" });
                    sil.noktasil(feature);

                }

                //set up the modification tools
                var DeleteFeature = OpenLayers.Class(OpenLayers.Control, {
                    initialize: function (layer, options) {
                        OpenLayers.Control.prototype.initialize.apply(this, [options]);
                        this.layer = layer;
                        this.handler = new OpenLayers.Handler.Feature(
                            this, layer, { click: this.clickFeature }
                        );
                    },
                    noktasil: function (feature) {
                        // if feature doesn't have a fid, destroy it

                        if (feature.fid == undefined) {
                            this.layer.destroyFeatures([feature]);
                        } else {
                            feature.state = OpenLayers.State.DELETE;
                            this.layer.events.triggerEvent("afterfeaturemodified",
                                                       { feature: feature });
                            feature.renderIntent = "select";
                            this.layer.drawFeature(feature);
                        }
                    },
                    clickFeature: function (feature) {
                        // if feature doesn't have a fid, destroy it

                        if (feature.fid == undefined) {
                            this.layer.destroyFeatures([feature]);
                        } else {
                            feature.state = OpenLayers.State.DELETE;
                            this.layer.events.triggerEvent("afterfeaturemodified",
                                                       { feature: feature });
                            feature.renderIntent = "select";
                            this.layer.drawFeature(feature);
                        }
                    },
                    setMap: function (map) {
                        this.handler.setMap(map);
                        OpenLayers.Control.prototype.setMap.apply(this, arguments);
                    },
                    CLASS_NAME: "OpenLayers.Control.DeleteFeature"
                }

);

Everythings is same with clickFeature. i wrote new function and i call it and i give this function feature which i stored. but i gave this error :Uncaught TypeError: Object # has no method 'destroy' How can i do that ?


texts: Tilestache WGS84 UTFGrid output isn't recognized by Leaflet UTFGrid

texts:
I'm attempting to use TileStache as a server for a UTFGrid created from a shapefile in WGS 84, but the Leaflet UTFgrid plugin doesn't seem to be able to use the data. There are no javascript errors or server errors, and when I directly request the json files, they seem to be in order, though there are strange characters () which might be causing problems. When using a similar setup with the shapefile converted to spherical mercator, everything works fine, except for the output being shifted too much to be usable.

I'm using the newest version of TileStache, Mapnik 3.0.12, the version of Leaflet included with Mapbox JS 2.4.0, and the latest version of the UTFGrid plugin.

My TileStache CFG:

{
 "cache":
 {
  "name": "Test",
  "path": "/tmp/stache",
  "umask": "0000"
 },
 "layers":
 {
  "test":
  {
    "provider": {"name": "mapnik grid", "mapfile": "test.xml", "layers":[[0,["ZCTA5CE10"]]]},
    "projection": "WGS84",
    "allowed origin": "*"
  }
 }
}

My Mapnik XML:

<?xml version="1.0"?>
<!DOCTYPE Map>
<Map srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">

<Layer name="layer" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">

    <Datasource>
        <Parameter name="type">shape</Parameter>
        <Parameter name="file">/home/zip/test2.shp</Parameter>
    </Datasource>
</Layer>
</Map>

My HTML:

<!DOCTYPE html>

<link href='https://api.mapbox.com/mapbox.js/v2.4.0/mapbox.css' rel='stylesheet' />


texts: Programmatically download (seamless) NED tiles?

texts:
I used to download 1/3 arcsecond elevation geotiffs from extract.cr.usgs.gov using web services.  But a few days ago, that stopped working, probably because USGS reorganized their site.  I tried emailing USGS and haven't heard back yet.

I could, of course, manually download states or degree tiles using the seamless server, but I want to do it programmatically so that users of my trip planner can get elevation data for their region without human intervention.  And SRTM isn't high enough resolution for bicycle/wheelchair trip planning.

Any ideas?


texts: How can I generate a multiband tif image using getDonwloadUrl?

texts:
I use the GEE API in a Python environment and I want to upload numerous small image to my local machine. I try to implement the solution suggested here. Before upgrading my code I wanted to do some test in EE code editor but I cannot make it work.
So I tried to use the ee method getDownloadUrl to retrieve multiband images on very small surface using the following code:

function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 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 image.updateMask(mask).divide(10000);
}

var aoi = ee.FeatureCollection('users/bornToBeAlive/FRA_Grid_2')
var dataset = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterDate('2020-01-01', '2020-01-30')
  // Pre-filter to get less cloudy granules.
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
  .map(maskS2clouds)

var visualization = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

//clip the dataset on the aoi
var clip = dataset.mean().clip(aoi)

// get the first tile
Map.addLayer(clip, visualization, 'S2')
Map.addLayer(aoi.first().geometry(), {}, "france")
Map.centerObject(aoi.first())

The displayd image is very snowy but I manage to see things.
now to upload this image I want to get a downloadable link:
var link = clip.select(['B4', 'B3', 'B2']).getDownloadURL({
    'name': 'toto',
    'region': aoi.first().geometry(),
    'filePerBand': false
});

print(link)

Unfortunately the link I get is giving me an image that is not scaled properly (the square is way too big) and I only get the same value for each pixel on each band.
Did I do something wrong? Does anyone knows how to use  this function?


texts:
I would like to apply a style to the KML that gets generated from the write() operation of OpenLayers.Format.KML. The current style that is on the vector layer does not get exported with this code:

function GetKMLFromFeatures = function (features) {
    var format = new OpenLayers.Format.KML({
        'maxDepth':10,
        'extractStyles':true,
        'internalProjection': map.baseLayer.projection,
        'externalProjection': new OpenLayers.Projection("EPSG:4326")
    });

    return format.write(features);
}

I would like to add something like this to the KML:

<Style id="PolyStyle">
    <PolyStyle>
        <color>ff0000cc</color>
        <fill>0</fill>
        <!-- boolean -->
        <outline>1</outline>
        <!-- boolean -->
    </PolyStyle>
</Style>

So that i can specify the fill property. How annoying is this:

I want this:


texts: Calculating daily travel distance from large dataset using QGIS?

texts:
My dataset is around 200 days worth of data on the daily distance travelled by a troop of monkeys. I've never properly used QGIS on my own so don't really know where to start. All of the days are together in one dataset.
Basically, what I need to do is calculate how far the monkeys travel each day.

The picture is a sample of my dataset. The monkeys were followed throughout the day and a waypoint of their daily travel was taking every 20 minutes. So, I have an update of how far they have moved every 20 minutes and basically want to calculate how far they have traveled throughout the whole day using these waypoints.


texts: GeoServer VendorOption for SLD to place labels overlapping and out of bounds

texts:
I use GeoServer 2.1.1. and have an SLD to serve a layer with labels in OpenLayers. I noticed that some labels are not displayed because other labels are to near. GeoServer wants to make your labels look nice. Weird enough: I don't.
Using the VendorOptions in the SLD I figured out some of the problems. I read the SLD cookbook.
I am now using:

<TextSymbolizer>
 ...
 <VendorOption name="spaceAround">-1</VendorOption>
 <VendorOption name="conflictResolution">false</VendorOption>
 ...
</TextSymbolizer>

These two options makes GeoServer draw multiple labels on top of each other even though it will render the labels unreadable (which is what I require).

The problem that still remains is the positioning of the labels. I have the labels set to:

<LabelPlacement>
 <PointPlacement>
  <AnchorPoint>
   <AnchorPointX>0.5</AnchorPointX>
   <AnchorPointY>0.5</AnchorPointY>
  </AnchorPoint>
 </PointPlacement>
</LabelPlacement>

But as soon as the labels are to near of the edge of the image/bounds the it no longer centers the labels but will anchor it left or right depending on the position.

My question if someone knows how I can make GeoServer ignore trying to make it look so damn nice (I know it sounds weird!) and just let my labels run out of the image/bounds.


texts: How to Make Lines Automatically Adjust same with Another Line (Route or Road)

texts:
I have 1000 lines and 1 route or road line. I want to automatically my lines adjust to route or road without erasing all data in the attribute table. My expected results, 900 lines will same with the route. I used QGIS 2.18

This is my lines

and this is the route/road line


texts: Dealing with ArcSDE connection limit exceeded?

texts:
We routinely exceed ArcSDE's connection limit, about once a week. Many of the "connections" are nothing but orphaned records in the SDE.PROCESS_INFORMATION table. Is there any way to clean up these connections periodically? I just service packed ArcSDE to 9.3.1 SP2, hoping this would resolve the problem. (It did not.) Currently I restart the arcsde service to clean up the connections, a pretty inelegant solution.

This is the error, by the way:

Failed to connect to database. Maximum number of connections to instance exceeded

My workaround:

My solution was simple: restart the SDE service weekly, and double the number of allowed connections. This is of course not a solution; shame on ESRI for shipping enterprise software that requires such a workaround. After all, it's still possible to exceed the maximum connections; I can only hope that our current level of use doesn't lead to this.

I altered the init.d script to prevent sdemon from prompting the user to confirm they want to stop the service (this is what the -N argument does). I added a restart option as well.

#!/bin/bash
#
# arcsde Init file for starting and stopping ArcSDE 9.3
#
# chkconfig: 35 90 40
# description: ArcSDE startup script

# Source function library.

. /etc/rc.d/init.d/functions

SDE_OWNER="sde"
SDEHOME="/home/sde/sdeexe93"

case "$1" in
start)
echo -n $"Starting ArcSDE:"
su - $SDE_OWNER -c "$SDEHOME/bin/sdemon -o start -p pwd"
echo "OK"
;;
stop)
echo -n $"Stopping ArcSDE:"
su - $SDE_OWNER -c "$SDEHOME/bin/sdemon -o shutdown -p pwd -N"
echo "OK"
;;
restart)
cd "$CWD"
$0 stop
sleep 5
$0 start
;;
*)
echo $"Usage: $0 {start|stop|restart}"
esac

Then I created a cron job to restart the service once weekly:

# restart arcsde once weekly at 5am Saturdays
0 5 * * 6 service arcsde restart

I also increased the maximum number of allowed connections from 64 to 128 by editing the $SDEHOME/etc/giomgr.defs file:

CONNECTIONS      128     # maximum number of connections
                         # NOTE:  On windows machines, you may need to
                         # increase server non-interactive desktop memory.
                         # Consult the ESRI support site for more information.

Then I imported the new settings:

$ sdeconfig -o import -f $SDEHOME/etc/giomgr.defs -i esri_sde -u sde

That was it. We'll see how it goes.


texts: Why is my Natural Earth map not shown correctly in QGIS? (strange circle; only Caucasian Sea shown)

texts:
Using the world map data from Natural Earth I have some problems with the correct displaying in QGIS.

Shapefiles added:

1:10 Admin Countries
1:10 Oceans

As you can see in my screenshot, there are some mistakes:

A strange circle appears (small red frame)
Russia ist not shown
The ocean layer is only shown for the Caucasian Sea (blue polygon)

On some scales or projections the ocean is illustrated but up to now, I have not detected when and why this happens. The added shapefiles are not edited or rendered dependend on scale.


texts: Cost distance of one raster using routes of separate cost distance in ArcGIS Desktop

texts:
Is there a way, using ArcGIS Desktop, I can calculate the cost distance for one raster of cost values but instead of using the least cost distances from that raster, using the routes/paths determined from a separate cost raster?
For instance, say I have a cost raster for energy expenditure based on slope but I want energy costs taking the least cost distances for another raster such the vegetation density. The product would essentially be the total sum of said energy values but while walking through the thinner vegetation if possible.
I need this for an entire surface rather than a path.


texts: Publishing oracle spatial table in geoserver

texts:
I published oracle view in geoserver. using this link. But when I see this in layer preview I get following xml response:

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://192.168.70.65:80/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" >   <ServiceException code="internalError">
      Rendering process failed
java.io.IOException
null
ORA-00904: &quot;NEPALCLUSTER&quot;: invalid identifier

</ServiceException></ServiceExceptionReport>

 CREATE OR REPLACE FORCE VIEW "myview"."d_DAILY" ("nepalCluster", "nepalclutter", "ZONE", "MI_PRINX", "MI_STYLE", "GEOLOC", "VAL") AS
  select
  c.nepalCluster,c.nepalclutter,c.zone,c.mi_prinx,c.mi_style,c.geoloc,NVL(coalesce(optz.VAL,opth.VAL),-1) VAL
from
CLUSTER_REGION c
left join
(select
*
from gDaily a) optz
on
optz.CELL_CLUSTER = c.nepalCluster
left join
(select
*
from zDaily b) opth
on
opth.CELL_CLUSTER = c.nepalCluster;


texts: Converting DataFrame to GeoDataFrame with Polygon

texts:
I would like to convert my DataFrame to GeoDataFrame, but every post is about points. I need to convert polygons. Coordinates are in df['geometry']. How can I do that?
My DataFrame is below:
I was trying:
df = gpd.GeoDataFrame(df, geometry='geometry')

But I'am getting error:
TypeError: Input must be valid geometry objects: 495495.84 516704.74 495487.95 516706.54 495490.65 516718.33 495498.54 516716.52 495495.84 516704.74


texts: Leaflet GeoJSON style function feature only colors one GeoJSON feature

texts:
I need to style every feature differently that I pull from a GeoJSON file.  However this anonymous style function only changes the style on the first feature it comes across and stops. I am new to Leaflet.  Have seen several demos but can't find a reason why this anonymous style function only changes the style of the first feature and not the rest. The getTableData function returns a color from table generated on page load.

Feature example

var countyData = {
  "type": "FeatureCollection",

  "features": [
{ "type": "Feature", "id": 0, "properties": { "ENTITYYR": 2010.0, "NAME": "COUNTY1", "FIPS": 25.0, "FIPS_STR": "49025", "COLOR4": 4 }, "geometry": { "type": "Polygon", "coordinates": [ [ Long List Of Coordinates here ] ] } }
, etc

    geojson = L.geoJson(countyData, {
        style: function(feature) {
    switch (feature.properties.NAME) {
        case 'COUNTY1': return {color: '#' + getTableData(feature.properties.NAME)};
        break;
        case 'COUNTY2': return {color: '#' + getTableData(feature.properties.NAME)};
        break;
        case 'COUNTY3': return {color: '#' + getTableData(feature.properties.NAME)};
        break;
        case 'COUNTY4': return {color: '#' + getTableData(feature.properties.NAME)};
        break;
        case 'COUNTY5': return {color: '#' + getTableData(feature.properties.NAME)};
        break;
    }
        },
        onEachFeature: onEachFeature
}).addTo(map);


texts: Exporting district level precipitation data for series of months using Google Earth Engine Python API

texts:
I am trying to export precipitation data for a range of months at the district level of Somalia using GEE's phython api. I'm using the following code
#import Somalia district shape file and convert to ee object

districts = geemap.shp_to_ee(r"C:\Users\Leslie\OneDrive\Raw GIS Data\Admin_District_PESS_2014\Original dbf File\SOM_Population_PESS_2014")

# import the RS products
chirps = ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD')

# Define time range
startyear = 2000
endyear = 2001

# create list for years
years = range(startyear,endyear);

# make a list with months
months = range(1,12);

# Set date in ee date format
startdate = ee.Date.fromYMD(startyear,1,1)
enddate = ee.Date.fromYMD(endyear+1,12,31)

# Filter chirps
Pchirps = chirps.filterDate(startdate, enddate).sort('system:time_start', False).select("precipitation")

# calculate the monthly mean
def calcMonthlyMean(imageCollection):
    mylist = ee.List([])
    for y in years:
        for m in months:
            w = imageCollection.filter(ee.Filter.calendarRange(y, y, 'year')).filter(ee.Filter.calendarRange(m, m, 'month')).sum();
            mylist = mylist.add(w.set('year', y).set('month', m).set('date', ee.Date.fromYMD(y,m,1)).set('system:time_start',ee.Date.fromYMD(y,m,1)))
    return ee.ImageCollection.fromImages(mylist)

# run the calcMonthlyMean function
monthlyChirps = ee.ImageCollection(calcMonthlyMean(Pchirps))

#function to apply reduceRegions over the districts feature collection
def eeAlgo(image):
    reduced = image.reduceRegions(collection=districts, reducer=ee.Reducer.median(), scale= 5000)
    return image

vals = monthlyChirps.map(algorithm=eeAlgo)

I am not sure if this is right especially since 'vals.getInfo()' does not indicate a property called district.


texts: Multiple Layers+Retrieval of data from DB?

texts:
I have some big & complicated task this time. I have 2 html pages(page1 & page2). Page1 is for viewing tha map. And page2 is for inputing the values(lat-lon) in a form. Any user can input values and any user can view the map(for now there are no security issues). page1 is provided with a text box and a button. text box(mostly a drop down menu) is used to select the layer and button is clicked to view the map. page2 is a simple form with boxes to give layer name,lat-lon values.

Many people suggested me to use Geoserver,postgis,postgreSQL...  but I dont understand how they can be used for my project

Task to be done:

whenever any user(lets say user2) inputs layer name and lan-lot values in page2, the values should be stored in a table(with table name as layer name).This storage can be done using PHP into PostgreSQL. Is there a better way to do this?
Now any user(say user1) can select the layer name and click on button to view the map, the markers of that layer(of values present in corresponding table) should be projected on the map. I have to give option to the user to view more than one layer at a time(say with different colored markers for layers).  Do I need Geoserver for this?

PostGIS is used to get the values in OL accepted format(JSON) from postgreSQL, but why do I need GeoServer?

In other words,how does geoserver & postgis help me to fulfill this task?

Your explanation in words helps me a lot instead of any links. Any examples(probably with some code) will also help me.


texts: How to get lat/long from Google Analytics when country but no city given?

texts:
About a year ago, I created a web-based map of the visitors to a web site. I've been asked to switch it from using Urchin and web server log data, to retrieving data from Google analytics.

When I was using the older data sources, I had the IP address of each visit, and I would use the IP address to look up  city name / country code / latitude / longitude.
Sometimes the city was unknown in the database I was using, but I would almost always get the country code, with a centroid latitude/longitude.

This is important data, because I would put markers on a map showing where visitors were from, and sometimes a large proportion of visits from a given country were from unknown cities. For example, I have results like

Country: Anguilla  Visits: 2  City: (not set) Latitude: 0.0000 Longitude: 0.0000

This is the only data point from Anguilla. So if I don't plot it, Anguilla is completely omitted.

Now with the Google Analytics Data Export API, IP address is unavailable. Thankfully, city name / country name / latitude / longitude are available most of the time, so I don't have to look those up by IP. But in the cases when the city is unknown, but country name is given, lat/long is given as 0.0000,0.0000.

How can I plot these country points on a map? It's difficult to look up country locations in a database because the country name is given, rather than a country code. Country names of course are somewhat variable. It really surprises me that Google would not make country code available at all, but I don't see it listed under the allowed "dimensions" for an API query.

I guess I could build up my own lookup table, fine-tuned to match Google's country names, but it would be subject to breakage if Google ever changed spelling or phrasing (e.g. "Democratic Republic of the Congo" vs. "Congo-Kinshasa").

Thanks for any suggestions.


texts:
How can i add a bridge as a normal street in the config of osm2po 5.1.0?

When i try to calculate the route from "333842" to "291763" using austria-latest.osm.pbf an error occurs about failing to calculate the route. After analyzing the route i found that there is a bridge within this route.

I tried to add this in the config file using

wtr.tag.highway.bridge       =   1,  26, 50,  car|foot

unfortunately without success. Any ideas how to solve this?


texts: Problem with shapefile created using Arcgis 10.3

texts:
I have create  a shapefile using Arcgis 10.3 . the shapefile contains attribute data in Arabic language . the problem came when I open it in Mapwell (an app  specializes in creating navigable maps) Mapwell displays it as series of question marks and other symbols , even if I set the encoding to Arabic .
Here another I have noticed , When I export the same data using Arcgis 10.0 it just works flawlessly.

I did some workarounds like changing the encoding scheme of some files (.cpg and .dbf) and other things, and nothing has worked


texts:
I have some kind of noob question: how do I plot spatialpoints dataframe in R using spplot (or ggmap) based on column values?

let's say we have this:

library(sp)
data(meuse)

v <- SpatialPointsDataFrame(meuse, coords = meuse[,1:2])

How would I plot 'copper' column here (with triangles pointed up or down based on a random threshold - let's say 31.0)?

I couldn't find a specific anwser so far.

Or how do I do that using proportional symbols?


texts: Cleaning shapes of polygon layer with overlaps and gaps using ArcGIS Desktop?

texts:
I have a layer of polygons where the shapes are quite messy with gaps and overlaps, see example:

It's clear that the polygons are supposed to be contiguous and are supposed to roughly follow the streets.

Does anyone know of a pre-made tool, package, script, etc, or programmatic strategy, that would allow me to clean a dataset like this?

I've heard of the Integrate tool

However this does not completely fix this problem since it modifies nodes of the shape all along it's boundary (not just along the sides where there are gaps/overlaps).


texts: For loop for access to a directory files in ArcPy

texts:
I have a Python Toolbox and I want to use a directory to read different shapefiles. But when I run the script It only works for the first shape.
fd is the directory and the code is:
fd= parameters[0].valueAsText
raster= parameters[1].valueAsText

mar="Z:\\Datos_centralizados\\OSM_coastlines\\water_polygons.shp"
out= "Z:\\EMSV-057_AOIs_Reshape\\5_Process\\1_TaskA\\ModelBuilder\\MB_pruebas.gdb"
#-------------------------

#leemos los shapes de la carpeta y los metemos a una lista

flist = []
for root,folder,files in os.walk(fd):
        for file in files:
            if os.path.isfile(os.path.join(root,file)) and file.endswith('.shp'):
                flist.append(os.path.join(root,file))

for z in flist:
    shape=z

    # Cogemos datos alfanumericos del AOI
    desc = arcpy.Describe(shape)

    # aplicamos la herramienta Erase para quedarnos con la tierra
    erase=arcpy.Erase_analysis(shape,mar, out+"\\Tierra", cluster_tolerance="")

    # Estadisticos
    tZONA=arcpy.gp.ZonalStatisticsAsTable_sa(erase, "FID", raster, out+"\\tabla_Zona", "NODATA", "SUM")


texts: How can I join attributes by nearest, filtering the results by an another attribute present in both tables?

texts: I have the centroid layer of climate polygons and I want to find the closest weather station of the same type of climate. The two layers (the one of centroids and the one of weather stations) contain a field with the type of climate. I am working with QGIS, the "Zanzibar 3.8" version.


texts:
I have created a route by hand using the Tracing Tool available in QGIS, placing the waypoint on top of an Open Street Map layer. I wanted to export all the coordinates to another program, so I saved them in a CSV file. The problem is that I can only see the latitudes and longitudes, but not the altitudes associated to each pair of latitude-longitude point. I think that my problem is that OSM doesn't work with the z-axis.

I want to create a track specifying only the latitude and longitude of each waypoint in a way that the program automatically tells the altitude of the road/street at each waypoint. Someone knows how could I do it?


texts: Using ogr2ogr to convert Natural Earth data to Polar Stereographic projection

texts:
I’m trying to build some maps of the North Pole area for my website. I’ve flown over it a few times and want to make a map of the flight paths I took. As you can see, a normal Mapbox map doesn't look so great for these polar routes:

I want a map like this, which I’ve recently learned is called a Polar Stereographic projection.

And I want to put it in a Mapbox map. Mapbox only supports Mercator projections, so if you want to work with different projections, they need to be cast back to Mercator somehow (I know this is tacky/hacky).
How can I do the following with GDAL/ogr2ogr:

Get Natural Earth Shapefiles
Transform to a Polar Steerograph projection like ESPG:3995
Make those coordinates work on a Mercator projected map like Mapbox uses
Get GeoJSON of the shapefile to go into Mapbox Studio with and build a map?

Previous research I’ve done:
I found an old post on the Mapbox blog about doing something like this, but the gdal commands are outlined generally rather than written out, and the associated Github repo runs all the shapefile changes with Postgres/PostGIS which is overkill, I’m just trying to run a few ogr2ogr commands to transform the shape file into gsojson.
I tried a few ogr2ogr commands with t_srs, a_srs, s_srs, but my geojson always comes out like this:
"geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ -6607704.23424893990159, 2768583.037965904921293 ], [ -6606072.226559722796082, 2767507.039690425619483 ],

Where the coordinates appear to be in https://epsg.io/3995 format, but since they’re not between -180 and 180, Mapbox can’t process them (I’m failing to get this file converted back to Mercator format I think)
I also found the Dirty Reprojectors app, but it doesn’t have data in Polar Stereographic projections, presumably because D3 doesn’t support that projection transformation.


texts:
I am trying to develop a GIS web application. I am using ArcGIS for desktop as the database. I  need to use existing ArcGIS Restful services but also need to publish some new services. I think that ArcGIS for server would be the best choice for me. But I cannot afford it.
And I want to know:

If there is a free lightweight version of ArcGIS for server?
Does a open source alternative such as GeoServer have APIs to use ArcGIS Restful services?


texts:
I have a shapefile which defines for me the area of each oil field in the North Sea (~300 fields). I also have a shapefile which defines the location of every well that has been drilled (~13000) wells. I am looking for a way to associate the wells to the fields.

I have managed to do this for a single field by splitting the oil fields shapefile into a shapefile for each oil field, then taking the intersection with the wells shapefile, and copy-pasting the attributes table out. It will take me a very long time to do this for 300 fields!

Is there a better way?


texts: Cannot connect to WMTS layer from NASA GIBS API in QGIS

texts:
QGIS Version: 2.18.3

I am trying to add some basemaps layers from NASA's GIBS API, but I cannot seem to get the URL right to connect to it.

According to the site, I think I should just be able to use the REST endpoint at https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/ Under "Add Layers from a WM(T)S Server" in QGIS. When I do this, however it is not recognized as a
WMS server, because I get the error "The server you are trying to connect to does not seem to be a WMS server. Please check the URL.".

If I try the GetCapabilities link at https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/1.0.0/WMTSCapabilities.xml, then I do get a list of services to choose from, but the map canvas is blank when I try to add anything.

The Layer Properties then shows this under the General > Layer Source tab:

contextualWMSLegend=0&crs=EPSG:4326&dpiMode=7&featureCount=10&format=image/png&layers=VIIRS_Night_Lights&styles=default&tileDimensions=time%3D2016-01-01/2016-01-01/P1Y&tileMatrixSet=500m&url=https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/1.0.0/WMTSCapabilities.xml

Specifically I am trying to add the "Night Lights" layer.


texts:
Can anybody guide me to a website or a source to download the Geographical UK road network data set. The data set should include the road distances, since I want to apply shortest path search algorithm on it.

I searched the web for such data but I could not find it in readable format.


texts: Imported shapefiles will not load correctly in Android app MapIT

texts:
I am able to import shapefiles created with ArcMap, but the features do not appear on my map.  When I zoom to them I am taken to the North Pole.  I assume this is a projection issue.

MapIT has a much more limited selection of projections, mostly European, than ArcMap and the programs do not share any.  I have tried reprojecting the .shp files created with ArcMap to similar projections available in MapIT to no avail.

For instance, I tried NAD1983 UTM to NAD1927 UTM which I assume is similar to EPSG UTM.  Running Android 7.1.2 (nougat) by the way.


texts:
I am just trying to follow Anita Graser's examples in https://anitagraser.com/2017/04/08/a-guide-to-geometry-generator-symbol-layers/ - where the "Creating a line between features in different layers" part caught my attention.

Is there a way to select more than one feature based on an adequate function using wildcards or expressions for numerical fields ... like ReturnAnyFeatrue('from_layer', 'with_field', >= 200) ?

get_feature  just returns the first element that matches the "attribute"/"value" settings - thus the following expression returns a set of lines from any object of the current layer to a specific STATION - i.e. 'Millennium Tower':

make_line(
centroid( $geometry ),
geometry( get_feature( 'Citybike', 'STATION', 'Millennium Tower' ) )
)

I am interestet in a way to draw geometry-generator-lines from current layers' points to any station starting with - for instance - "Mi" based on an expression like 'Mi*'
is there a way that just makes use of geometry generator and simple expression ?

(unfortunately, there is still a lack of documentation and examples concerning geometry generator ... )


texts: Splitting multigeometries in a Virtual Layer (SQL query)

texts:
I have a virtual layer already set up and working, the only thing that could make it better is getting rid of multipart geometries. The query is:
SELECT ejido, circ, radio, mzna, cc, ten, registrado,
    RANK() OVER (ORDER BY CalcRegistrado(ten,registrado)) AS color,
    st_buffer(st_union(geometry),0.1) AS geometry
FROM "VIRTUAL_ORIGEN_PARCELAS"
WHERE registrado IS NOT NULL AND registrado <> ''
GROUP BY circ, radio, cc, ten <> 'S', mzna, CalcRegistrado(ten,registrado)

The buffer is meant to join polygons because I have to work with spaghetti layers, but sometimes parcels are indeed separated by other parcels and I need them to be separate polygons in the resulting layer. I have tried to generate a numbers table and selecting from a join with my geometry result (as seen here) but with no success.

To do so, first create a table containing the numbers from 1 to
100,000 (which should be plenty to accommodate the number of
individual elements in a multielement geometry), as follows:
SELECT TOP 100000 IDENTITY(int,1,1) AS n
INTO Numbers
FROM MASTER..spt_values a, MASTER..spt_values b;

You can then join from
your spatial table to the numbers table to retrieve each individual
geometry using a query as follows:
SELECT
  id,
  geom.STGeometryN(Numbers.n)
FROM
  #MixedData JOIN Numbers
  ON Numbers.n <= geom.STNumGeometries()

How can I apply that to my query?


texts:
I have some mxds from a few years back that I have been tasked with re-purposing.  I need to recreate the original symbology, now missing, which was apparently done by loading a lyr file from parts unknown.  Is there any way I can tell which lyr file was used to create the symbology?

If I create a legend it pulls the symbology I need, with the proper names & colors, but the source data doesn't contain the basis for this symbology.  What the source data has is a massive field of numbers, which in the lyr file (I assume) were encoded into what I need (ie 174367 = Room X).  I hope I won't have to manually recreate this encoding.


texts: How to import handdrawn polygon using json file in earthengine

texts:
Is there any way to import a data, say handdrawn polygons (saved as coordinates) as json or text files? I have hard time reading a file into the code editor, even ones created locally, as the code editor doesn't allow use of packages such as FileReader and XMLHttpRequest (that are generally used in javascript programming). So is there any way to save to file, and read data from text/json files in earthengine code editor? other than importing them using asset manager which restricts us to import only raster and shape data.

edit: I am new to the earthengine or in fact any gis programming. But I am not sure why the question is unclear, it is simple. In every programming language/library there is a way to save/read data locally or to drive (in this case) using standard text, csv, json etc., formats. So is there a way here or are we restricted just to use asset manager to import the two types of data? (if you are aware even the python interface supports that (as it runs locally), which I came to know recently). Because it is very obvious need when I program unless my whole code going to be written in the code editor, which is unlikely, especially if I want to do any serious computing involving other 3rd party packages, plus it feels the programming is very restricted when using code editor.


texts: Extracting single band from multiple band raster with PyQGIS?

texts:
I'm approaching naw in the script and I'm trying to select many multi-band files from a folder and extract a single band from the file, putting everything in an output folder with the single-band transformed files.

Currently I concluded the following script:

 import sys, os

 output = "C:/../folder"

 layers =[layer for layer in QgsProject.instance().mapLayers().values()]

 for lyr in layers:
    if (lyr.type()==QgsMapLayer.RasterLayer):

    cmd1 = "gdal_translate -b 1 "+ lyr.source() + "C:/..raster_input" + outputDir + "C:\..output" + lyr.name() + "test.ecw"

    os.system(cmd1)

    cmd2 = outputDir + "C:\..output" + lyr.name() + "test.ecw"

    iface.addRasterLayer(cmd2)

But it does not work!


texts: ArcObjects C# memory management - Out of memory exception

texts:
I'm a noob in C# and ArcObjects (started two weeks ago).

I'm trying to write a stand-alone console program that will change label placement quality in all data frames for all MXD files in a given folder. What I've got now kind of works, but after 40 or 50 MXDs that are changed, my program throws "out of memory exception". I suppose that standard garbage collection fails in this case. When I watch memory usage in task manager for this process, it fluctuates (so some memory is freed from time to time) but with a steady tendency upwards.

I have searched Google and SDK help as well for this topic. In general they say that all the COM objects created by oneself should be explicitly destroyed, as well as cursors and enumerators.

See the code below. The only time I use new explicitly is when creating MapDocument object and I immediately add it to ComReleaser - so, the way I understand it, it should be taken care of at the end of using block. I don't have any enumerators nor cursors, as I deal with layout only. What else should I release?

foreach (System.IO.FileInfo file in files)
    {
        counter++;
        if (!log.Exists(file.FullName))
        {
            using (ComReleaser cr = new ComReleaser())
            {
                bool changed = false;
                Console.WriteLine("({0}/{1}) Opening: " + file.FullName, counter, total);
                //Open mxd
                MapDocument mapDoc = new MapDocument();
                cr.ManageLifetime(mapDoc);
                mapDoc.Open(file.FullName);
                //Iterate data frames
                for (int i = 0; i < mapDoc.MapCount; i++)
                {
                    Map map = (Map)mapDoc.Map[i];
                    Console.WriteLine("Dataframe {0}", map.Name);
                    IAnnotateMap2 annoMap = (IAnnotateMap2)map.AnnotationEngine;
                    IMapOverposter overposter = (IMapOverposter)map;
                    //check if Maplex labeling engine
                    if (overposter.OverposterProperties is IMaplexOverposterProperties)
                    {
                        IMaplexOverposterProperties overposterProps = (IMaplexOverposterProperties)overposter.OverposterProperties;
                        Console.WriteLine("Current placement quality: {0}", overposterProps.PlacementQuality.ToString());
                        //change placement quality if necessary
                        if (overposterProps.PlacementQuality != quality)
                        {
                            Console.WriteLine("Changing...");
                            overposterProps.PlacementQuality = quality;
                            changed = true;
                        }
                        else
                        {
                            Console.WriteLine("No need to change.");
                        }

                        Console.WriteLine("Current placement quality: {0}", overposterProps.PlacementQuality.ToString());
                    }
                    else
                    {
                        Console.WriteLine("Dataframe does not use Maplex.");
                    }
                }

                //Save document if needed
                if (changed)
                {
                    mapDoc.Save();
                    Console.WriteLine("Saving...");
                }
                //close mxd
                Console.WriteLine("Closing...");
                mapDoc.Close();
                log.Add(file.FullName);
            }//ComReleaser
        } else
        {
            Console.WriteLine("({0}/{1}) Skipping: " + file.FullName, counter, total);
        }
    }


texts: Handling an irregularly shaped elevation raster file

texts:
We are trying to process a set of elevation rasters of increasing resolution. We have a 30m resolution file that covers a large area, and a 5m resolution file that looks something like this image.

There is an area of valid elevation, but the rest of the file contains 'null' elevation data. We need to blend this elevation with the lower resolution and avoid visible artifacts. If we had a regular rectangle elevation file, then I'm pretty sure that any number of tools would do the trick - QGIS, gdalwarp etc..., but I haven't found a recipe that will work completely in our case.

I would be perfectly happy if we could upsample the lower resolution elevation to fill in the null area of the 5m file, and blend between those areas. That would probably be the best outcome.

I have been trying gdalwarp which has given me fairly decent results, doing something like this:

gdalwarp -srcnodata -32767 -r cubic  30m.tif 5m.tif outputcubic.tif

That's not too bad, in the merged file it upsamples the null elevation areas of the 5m with data from the 30m. I might be able to live with this, but it would be even nicer if the line between the higher res elevation and the upsampled lower res elevation were blended.  I see blending options for gdalwarp, but aren't these options for cutlines?

Any ideas about how to accomplish this with either gdalwarp - or with something else?


texts: Why will arcpy.Project_management not create output?

texts:
I have 2 feature datasets in a .gdb (CNMS_Inventory and CNMS_Inventory_Albers). I set the dataset for CNMS_Inventory to 4269(NAD83) when I created it and it contains a feature class called S_Studies_County_Intersect_Merged. My goal is to basically change the projection of S_Studies_County_Intersect_Merged to North America Albers Equal Area Conic and write the output to the dataset (in the same gdb) called CNMS_Inventory_Albers (which was set to North America Albers Equal Area Conic). I first tried just using CopyFeatures but it through an error:

cnms_inventory = os.path.join(all_regions_gdb, "CNMS_Inventory")
new_studs = os.path.join(cnms_inventory,"S_Studies_County_Intersect_Merged")
albers_ds = os.path.join(all_regions_gdb, "CNMS_Inventory_Albers")
albers_studs = os.path.join(albers_ds, "S_Studies_Ln_Merge_Albers")
albers_unmapped = os.path.join(albers_ds, "S_Unmapped_Ln_Merge_Albers")
arcpy.CopyFeatures_management(new_studs, albers_studs)

  ExecuteErrorTraceback (most recent call last)
       in ()
        2 albers_studs = os.path.join(albers_ds, "S_Studies_Ln_Merge_Albers")
        3 albers_unmapped = os.path.join(albers_ds, "S_Unmapped_Ln_Merge_Albers")
      ----> 4 arcpy.CopyFeatures_management(new_studs, albers_studs)
     C:\Program Files (x86)\ArcGIS\Desktop10.5\arcpy\arcpy\management.py in
      CopyFeatures(in_features, out_feature_class, config_keyword, spatial_grid_1, spatial_grid_2, spatial_grid_3)
      2584         return retval
      2585     except Exception as e:
      -> 2586         raise e
      2587
      2588 @gptooldoc('DeleteFeatures_management', None)
      ExecuteError: ERROR 000224: Cannot insert features
      Failed to execute (CopyFeatures).

I figured that this maybe due to the fact that the output feature, new_studs (S_Studies_Ln_Merge_Albers) needs to be projected into Albers first (I was thinking I maybe cannot just write a NAD83 projected feature to an Albers projected dataset). So I tried the following which gives be another error:

sr = arcpy.SpatialReference(4269)
projection = arcpy.SpatialReference('North America Albers Equal Area Conic')
arcpy.Project_management(new_studs, albers_studs, out_coor_system=projection, in_coor_system=coord_sys)

  ExecuteError: ERROR 000210: Cannot create output
  T:\CCSI\TECH\FEMA\Datasets\CNMS\FY18Q3\CNMS_AllRegions_FY18Q3.gdb\S_Studies_Ln_Merge_Albers
  Failed to execute (Project).

Looking at this error, it seems like Project is trying to write my output to "T:\CCSI\TECH\FEMA\Datasets\CNMS\FY18Q3\CNMS_AllRegions_FY18Q3.gdb\S_Studies_Ln_Merge_Albers" instead of "T:\\CCSI\\TECH\\FEMA\\Datasets\\CNMS\\FY18Q3\\CNMS_AllRegions_FY18Q3.gdb\\CNMS_Inventory_Albers\\S_Studies_Ln_Merge_Albers", which basically leaves out the feature dataset CNMS_Inventory_Albers. Can anyone see where I am going wrong here or how I can fix this?


texts:
I started using QGIS yesterday since ArcGIS refuses to work, and I have my difficulties with it.

In short, I want to extract the information on GCP from this netCDF file, so that I can perform geo-spacial analysis: http://dataservices.gfz-potsdam.de/pik/showshort.php?id=escidoc:2740907 (download button is in the top left)

I have only experience working with shapefiles, which is why I would prefer to transform it to shapefile format. However, I would also be happy with a text-file featuring xy values and 26 additional columns for the GCP values (the data features GCP per cell from 1850 to 2100 for each decade -> 26 decades).

Could you kindly help me do either with QGIS?

N.B.: The dataset is huge and QGIS needs ages just to load it. I would not mind working with a slightly lower resolution, if that would make it easier to handle.


texts:
Is there a way to figure out if every tile of an ol.source.TileWMS is loaded?

My attempt so far is the following:

var tileColl = new ol.Collection()
tileColl.on('change:length, function(e){
     if(e.target.getLength === 0){
        console.log('ready!');
    }
});
tileWms.on('tileloadstart, function(e){
    tileColl.push(e.tile);
}
tileWms.on('tileloadend, function(e){
 tileColl.remove(e.tile);
}

Unfortunately, this fails because sometimes all tiles in the list are loaded before all tileloadstart-Events are fired -eventually the collection becomes empty. (Probably because the number of request send by the browser is limited)

I removed all interactions, so the user cannot move the view. I know I can solve the problem by using (for example) ol.source.ImageWMS but I prefer using ol.source.TileWMS because it has some visible progress when the tiles are loading.

Some ideas I had but could not make them work:
Calculate the number of tiles that are to load,
Figuring out the last tile which is to load.

Is there way to achieve this?


texts: Convert Multipolygon Geojson to multiple geojson polygons?

texts: I see a lot of questions going from Polygons to MultiPolygon, but is there some easy way to go the other way around? It would also be useful to preserve the attributes from the MultiPolygon to apply to the new Polygons.


texts: QGIS object selecting if any segment length is longer than

texts:
I have a polygon layer with generated triangles covering for example district or city area. I need to select all the triangles if any triangle (object) segment is longer than 100 meters. Picture: selected triangles has a that segment (or many segments)

I prefer not to explode triangles to lines, calculating length and selecting lines, and run spatial query with triangles because every exploded line may touch a triangles with shorter segments.

Tnx. Tomek


texts: Creating MXD From Standalone ArcPy Snippet?

texts:
Is it possible to create some MXD from a list using Arcpy stand alone snippet? For example if I have a list like

counties = ['Benton', 'Caroll','Polk','Scot']

can I run a code like

import arcpy

for mxd in counties
    MXD = arcpy.mapping.MapDocument(r"C:\\"+ mxd+ '.mxd')
    MXD.save()


texts: Adding a button from an add-in to a TOC context menu

texts:
I would like to add a button to a group layer context menu in the table of contents. The button is from an add-in that also contains an extension. Should be straight-forward enough. The extension monitors the IActiveViewEvents_Event events and checks for a group layer being added by a specific name. When a layer is added, the IActiveViewEvents_ItemAddedEventHandler() is fired and calls a method that adds the button to the group layer context menus if the name is right. That mechanism is working. The snippet of code that is failing is below:

// Set a uid for the group layer context menu.
UID gUID = new UID();
gUID.Value = "esriArcMapUI.GroupLayerContextMenu";   // "{863A0D98-73DC-4331-8658-ED0E22247E36}"

// Get the command bar denoting that context menu.
ICommandBar groupLayerContextMenu = ArcMap.Application.Document.CommandBars.Find(gUID, false, false) as ICommandBar;

// If there is a group layer context menu.
if (groupLayerContextMenu != null)
{
    // Change the UID to the add-in command.
    gUID.Value = ThisAddIn.IDs.TestContextMenuButton;

    ICommandItem cmd = groupLayerContextMenu.Find(gUID);

    if (cmd == null)
    {
        cmd = groupLayerContextMenu.Add(uid, ref missing);  'Command not available' error thrown here...
        cmd.Refresh();
    }
}

To me, it seems that the command is not found because ICommandBar.Add(uid) is trying to dig for it in the registry but because it resides in an add-in it has not been added to the registry. Is there another way to add a button, or any item for that matter, when it lives in an add-in? Or am I missing something simple?


texts: Tools for importing historical maps into modern GIS systems?

texts:
What tools exist to derive geometries for use with modern GIS systems from older scanned maps?

I am referring more to libraries that generate actual data rather than present an interface for manually editing or tracing features.

If no such tools exist, are there any academic papers or research (or shameless hacks!) from related fields that might shed light on the problem?


texts: How to get a GML overlay on Google Maps?

texts:
Anyone know how to get a GML file to show as an overlay on Google Maps? To be honest I don't know my GML from my KML, or my ETRS89 from my WGS84, as this topic is completely new to me.

A bit of reading that I understood suggests that I have to convert GML to KML first, though it seems that process may corrupt the data and move the polygons, allegedly by up to 15m, which isn't acceptable.

I've downloaded some free software: GML Viewer from Snowflake Software, which seems to show the GML file, but it's a bit klunky, and what I really want is to overlay it over Google Maps' hybrid view, so I can see roads and the satellite image too.  Maybe QGIS which I've also read about would be a better tool to use?

Any ideas or pointers to tutorials very welcome.


texts: Identifying geographic coordinate system?

texts:
I'm looking at the following map and I can't figure out what kind of coordinate system is being used. This isn't lat/long. I assume once I know which coordinate system it is, I can easily search for a converter online.

As one can see, the coordinates are in the range of 784600N to 786000 N and 295800E to 297400E.

Any ideas?

The area in question is Gabon. Around the lat/long coordinates of  1.69°S, 13.41°E


texts: I have a problem with the TIN. I just started using QGIS because I no longer have access to ArcGIS. I want to create a reservoir bathymetry map from depth survey points. I created a polyline and converted it to a polygon so I can contain the interpolation within the shoreline boundary but when I run the TIN, it does not recognize the shoreline polygon as the extent and it interpolates out of that. What is wrong? Also is the TIN the most suitable for bathymetric interpolation or is there any other method? It is for Lake Diefenbaker 180 km long, more similar to a river than a lake. Thanks


texts: Adding points when joining two datasets together in ArcGIS Desktop

texts:
I am trying to join an excel file with a shapefile based on a zipcode. My problem is that I have multiple persons with the same zip code, but it joins only 1 of those persons.

Is it possible to automatically add multiple points when trying to join several persons from the excel file with my shapefile (which only has one zip code point)?


texts: Moving attributes from one column to another on joined layer using ArcPy

texts:
In ArcMap you can calculate a field on a joined table to just move the attributes from one column to another with just an = sign. Since my layer has feature linked anno on 7 different layers the calculation has to be done in an edit session. I tried this and it ran over three hours and not sure if it was going to finish. This won't work since I have to do this on about 30 fields.

I am trying to do this in ArcPy. I was using arcpy.CalculateField_management. I don't think that works based on the error I received. The error is invalid property. I think this means that it wants to make a calculation not just move attributes.

Here is my sample code:

import arcpy
from arcpy import env

anno_conn = "\\\\any database"
env.workspace = anno_conn

#Variables
Table_A = "data_1"
Table_B = "data_r3"
Layer_A = Table_A + "_lyr"
Layer_B = Table_B + "_lyr"

try:
    #join tables
    arcpy.MakeFeatureLayer_management(Table_A, Layer_A)
    arcpy.MakeFeatureLayer_management(Table_B, Layer_B)

    arcpy.AddJoin_management(Layer_A, "visi_link", Layer_B, "prior_visi")

    # Using an edit session - feature linked layer can't be edited outside an edit session
    # inside try edits will be saved if successful, will roll back if an exception is thrown
    with arcpy.da.Editor(env.workspace) as edit:
        arcpy.CalculateField_management(Layer_A, Table_A + ".link_id", "!" + Table_B + ".link_id!", "PYTHON")

except Exception, err:
    traceback.print_exc()

It seems to work up until the Calc field line.


texts: Adding a shapefile with a double-click to an opened project in QGIS in Nautilus on Ubuntu

texts:
On my Ubuntu computer, I have GIS file types like .shp associated with QGIS in Nautilus. When I double click on a .shp while I have a QGIS project open, it opens a new project in a new instance of QGIS, even if I have a project open.

On my Mac, if I have QGIS running and I double-click a .shp in Finder, the file will be added to the existing project.

Other files associated with other programs will open in a running instance rather than a new instance: text files will open in a new tab in my text editor, for example.

I'd like QGIS to add the shapefile in my opened project instead of launching a new instance. How can I do that ?


texts: Transferring ArcGIS Server directories to another location (D drive) on AMI EC 2

texts:
Updated question:

I am using ArcGIS Server 10.1 + ArcGIS desktop 10.1 on AMI EC2 instance & I am facing space problem so trying to move existing server directories to different drive.

Currently our dept is under maintenance so I am copying “arcgisserver” folder (Currently it is on C:\arcgisserver) to my “D drive”.

So my question is that which things should I take care (like data ,rest service etc) ?

I am familiar with ArcGIS Server 10 and it was easy to perform above settings in 10.0 as compare to 10.1

Currently I am going through ESRI web help and researching on it

Any help will be great !!!  Thanks in advance : )


texts: ArcGIS Online & Mobile Apps Projection Support?

texts:
What is the level of projection/co-ordinate system and transformation support accross the ArcGIS Online family of apps ;

ArcGIS Online Web Map Viewer (Desktop)
ArcGIS Online Web Map Viewer (Mobile)
ArcGIS for Windows Mobile
ArcGIS for Android / iOS
ArcGIS Collector

Is there any documentation or a support matrix such as ESRI provide for data sources ?


texts: Using r, readShapeSpatial fails with simplified mapshaper.org shape

texts:
So I'm mapping out the postal regions of my country and it's an enormous, 5.4million row dataframe.  I don't need every single little twist and turn of the boundries, so I was super excited to discover mapshaper.org.  Seemed easy enough, I just threw my .shp file into it, and it spit back a file 10% of the size!  Excellent!

Nope..  now it doesn't play with readShapeSpatial...  Hopefully I won't need to provide the shapefile (i can track it down if that parts really necessary)

Original code, works:

canada.shp <- readShapeSpatial("canada", verbose=TRUE, delete_null_obj=TRUE)

New map, doesn't:

canada.shp <- readShapeSpatial("simplecanada", verbose=TRUE, delete_null_obj=TRUE)
Shapefile type: Polygon, (5), # of Shapes: 1621
Error in `rownames<-`(`*tmp*`, value = c("0", "1", "2", "3", "4", "5",  :
  attempt to set 'rownames' on an object with no dimensions
In addition: Warning message:
In .Map2PolyDF(Map, IDs = IDvar, proj4string = proj4string, force_ring = force_ring,  :
  Null objects with the following indices deleted: 895

traceback() reveals:

4: stop("attempt to set 'rownames' on an object with no dimensions")
    3: `rownames<-`(`*tmp*`, value = c("0", "1", "2", "3", "4", "5",
   "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
   "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27",
   "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38",
   "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49",
   "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", (ETC!!!!)
2: .Map2PolyDF(Map, IDs = IDvar, proj4string = proj4string, force_ring = force_ring,
       delete_null_obj = delete_null_obj, retrieve_ABS_null = retrieve_ABS_null)
1: readShapePoly("simplecanada",
       verbose = TRUE, delete_null_obj = TRUE)


texts: Comparison of geocoding solutions (such as Gisgraphy & PostGIS 2.0)

texts:
As I am sure many of you are aware, PostGIS has finally released Version 2.0 and included within this version they have also built a geocoder (see: http://postgis.net/docs/Extras.html).

Although I realize V2 was just released I am curious to know if anyone has had the opportunity to compare other opensource solutions such as gisgraphy.com.

I personally am interested in understanding the differences in terms of speed, features/options, accuracy, and the ease of being able to use ones own dataset(s).


texts: Undo button undoes last couple of edits instead of last one in ArcMap?

texts:
I've made an add-in in ArcPy, using add-in assistant, to update attributes of selected features. It has two buttons to populate either 'U' or 'D' as attribute.

The Python code is:

import arcpy
import pythonaddins

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

        if sel != "":
            if arcpy.Describe (sel).FIDSet:
                with arcpy.da.UpdateCursor (sel, "Remark") as curs:
                    for row in curs:
                        row = ("D",)
                        curs.updateRow (row)
                    del curs
                    del row
            else:
                print("No selection")
            arcpy.RefreshActiveView()
        else:
            print("No selection")

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

        if sel != "":
            if arcpy.Describe (sel).FIDSet:
                with arcpy.da.UpdateCursor (sel, "Remark") as curs:
                    for row in curs:
                        row = ("U",)
                        curs.updateRow (row)
                    del curs
                    del row
            else:
                print("No selection")
            arcpy.RefreshActiveView()
        else:
            print("No selection")

This add-in works perfectly. The only problem is that whenever I press "Undo/Ctrl+Z" it undoes last 2-3 edits done by this button instead of only last edit. However, if I manually update something before pressing Undo, then this problem doesn't occur.

How do I rectify this problem?

I am using ArcMap version 10.3.


texts: How to use DateDiff() Function across day transition in ArcGIS Field Calculator?

texts:
I am trying to calculate the time difference between two fields using the DateDiff() function in minutes. I have successfully calculated the times that fall within the same day but I am having trouble with the times that are occurring in day transition (after 11:59 PM). For example: My date fields are as such,

MM\dd\yyyy hh:mm:ss tt

So between these two times

01\17\2007 11:57:08 PM and 01\18\2007 12:15:01 AM

The field calculator would apply a huge negative number. I am thinking a conditional statement is needed for basic arithmetic but am at a loss.

Please be as detailed as possible. I am not a developer.


texts: Is there a way to embed the filter UI into a map (for the user to use)?

texts:
Is there a way to embed the dynamic filter UI into a map for the user to use? I want them to be able to filter data in the same way that I can do it from my UI.

Is the first map on the following page just a "photoshopped" image? Or is there a way to let the user filter the data on their end?
http://cartodb.com/analyze


texts: How to set a minimum point snap distance?

texts:
There are two points: A and B.  I need to snap any points within a 65 m search radius (blue buffer) of point A.  However, the snap distance of B needs to be exactly 5m apart from A (i.e. located anywhere on the red buffer).  I know the Snap or Integrate tools can perform the snapping, however, I need to find a programmatic solution to set the minimum snap distance.  How can I set a minimum snap distance for this type of analysis?  I will need to eventually scale this up to accommodate thousands of points.

Desired Outcome (below)


texts: "Extract Error writing to file oraociei1.dll" while installing QGIS version 2.0.1

texts:
I have been desperately trying to install QGIS version 2.0.1 Dufour, but in the middle of the installation I always get an error ("Extract Error writing to file oraociei1.dll"). And automatically installation is aborted ("Setup was not completed succesfully").

I have tried so many times and also tried to download from different websites, not only from qgis.org, but nothing helps. Sooner or later I get the same type of error and that was it...

Any ideas what I could do?


texts: Is there a way to use the GDAL C++ API to find areas of overlap between two SAR images?

texts:
I have a bunch of SAR images (about 1 GB each) on which i need to do some image processing. Each image is georeferenced and there is significant overlap between each pair of images in the set ie. one piece of terrain appears in 2 or 3 of the images. To reduce the amount of work in the image processing step, is there some way to extract out the approximate areas of overlap between the images ?  Specifically, does the GDAL C++ API provide any such functionality ?

I have tried reading the documentation of the API but i couldn't find anything.


texts: Getting link or file from uploads folder of ArcGIS Server feature service?

texts:
I have succeeded in adding a file to arcgis server with the upload method of a feature service on arcgis server. But i can't see how to download what i have uploaded to.

How can I have a link (like the printTask button) or get the file ?


texts: Auto increment by group in ArcGIS ModelBuilder

texts:
I want to calculate auto increment a field name by group order Column "ID" into column "Name" in ModelBuilder and my input is shapefile.
Example:


texts: tilemill - styling map similar to google at max out zoom level

texts:
I have osm bright in tilemill
When i zoom all the way out to the world, the style on the maps looks very bland and basic, I was hoping to achieve something like how google maps look (I'm not after a direct copy), also mapquest has something similar

I.e. when you zoom all the way out, google maps looks a bit more interesting (nicer colours and shading) instead of one bland colour

My questions are
1) is there something open source that already has a less bland style (not sat photo view) - use google maps max zoom out to see what i mean

2) is it possible to combine my existing osmbright skin with something else but which only kicks in at certain zoom levels ?

Edit
Here is a good example of something i'd like
https://cloud.maptiler.com/customize/#streets/0.76/0/120.5


texts:
I have ArcGIS Pro (2.3.x) installed on my PC. I have python scripts which call ArcGIS toolbox models and runs them. Everything works fine.

Today I installed GDAL and Python bindings for GDAL. I also followed installation instructions and created a new system variables for GDAL_DATA (C:\Program Files (x86)\GDAL\gdal-data) and GDAL_DRIVER_PATH (C:\Program Files (x86)\GDAL\gdalplugins). I also edited the PATH variable and added 'C:\Program Files (x86)\GDAL' as an environment variable. Gdal runs fine on cmd prompt.

However, when I went back to run my arcpy scripts, I had a load of errors like 'can't load requested DLL GDAL C:\Program Files (x86)\GDAL\gdalplugins\xxxx.dll is not a valid win32 application'. It seems like installing a separate GDAL install has pointed ArcGIS in the wrong direction for any GDAL related queries? I didn't even know ArcGIS used GDAL...

When I then uninstall GDAL, the ArcGIS scripts work fine.

Any advice to what has gone wrong here and what I should do to install an independent GDAL on my ArcGIS machine?


texts:
I have a database with Easting and Northing points of New Zealand. The data structure is as follows

Easting  |   Northing

2595425  |   6622565

2594585  |   6622595

2595000  |   6622085

2594500  |   6622005

2595300  |   6622915

I need to create shapefile from this data. My Final goal is to convert this data to WGS84 system.


texts: What are data visualization and analyses pros/cons of QGIS over ArcGIS for Desktop?

texts:
Many of my colleagues use GIS on a very sporadic basis for simple data visualization and basic analysis (clipping, buffering, etc.) and I am hoping to transition these non-core GIS people to QGIS so that we aren't relying solely on ArcGIS and it's associated licensing costs.

I've used QGIS a bit in the past, but was hoping that some of you could explain the pros/cons of using it over ArcGIS. What if someone wanted to do more advanced raster analysis? What is the support like for data formats? How about interfacing with databases?

I'm just looking for a high-level overview and any hints/tips/pitfalls you may be able to provide information on.


texts: Unable to install plugins in QGIS 3.10 for Mac. No python support detected

texts:
I'm aware this question has already been asked numerous times in various forms, but I'm not finding any solutions that work for me. "No python support detected" in QGIS.

Issue:

Tried:
1. Updating my .bash_profile with the following, per the solution in the comments to this question: Python plugins not available, OS X 10.10.5, QGIS 2.8.4, which did not work:

Edit: Adding those lines ^ to my bash profile (and not undoing it when it didn't work for me) caused me issues later on. When I tried to use Python 2 I would get an error that my PYTHONPATH was set to QGIS Python 3. I forgot that I had made this change, so I was baffled for quite a while until I remembered. facepalm

2. A commenter on this question: Enabling (python) plugins on QGIS3 (Linux)? suggested installing qgis-python but I have been unable to find this package to try this solution.

Versions:
macOS Sierra 10.12.6, QGIS 3.10, Python 2.7 (original OS install)

Update 1:
Per this thread on a QGIS Github forum, I think my problem may be my Mac OS. I didn't realize that the version of QGIS I installed (3.10) is for Mac High Sierra 10.13 and newer. I only have Sierra 10.12. I'm going to install a newer MacOS and see if that works. Will post an update...

Update 2:
Yes! Installing a newer MacOS worked! I'm going to add it as an answer. Hopefully this will help someone else with an older MacOS who didn't catch the install requirement.