texts: Label properties not switching to selected layer for labelling

texts: In ArcGIS Pro my labelling properties is "stuck" on one layer i.e. not switching to allow other layers, fields and classes to be manipulated (the field do not change based on the selected layers attributes either). The attached screenshot even shows that this happens on projects where the layer "stuck" in layer properties isn't even within the project I have open. Is this a bug in the labelling?


texts: Getting Borders as SVG using Peters Projection?

texts: Is there a place where I can easily download the borders of all countries and continents as SVG using Peters Projection? Or what is the easiest way to get those? (I haven't installed anything on my machine related to GIS right now...)


texts: How to plot points on maps using ggplot2 and R?

texts:
Thanks for help me in GIS 101 Problem #1, now I have geocoded a few hospital in Connecticut using google map, I have able to visualize them in QGIS, now I am trying to do it in R.
ct <- readShapeSpatial("housect_37800_0000_2010_s100_census_1_shp/wgs84/housect_37800_0000_2010_s100_census_1_shp_wgs84.shp")

ct_mod <- fortify(ct,region="SLDLST10")
chart <- ggplot(data=ct_mod,aes(long,lat,group=group))
chart <- chart + scale_x_continuous(limits=c(-73.8,-71.7),breaks=seq(-74,-71,0.1))
chart <- chart + scale_y_continuous(limits=c(40.9,42.1),breaks=seq(40,43,0.1))
chart <- chart + geom_polygon(fill="grey80")
chart <- chart + geom_path(color="white")
chart1 <- chart + geom_point(data=hosp.list,aes(x=coord_x,y=coord_y))
chart2 <- chart + geom_point(aes(x=-72,y=41))

My hosp.list is a dataframe with hospital name, address, lat, long.
Now I can show chart and chart2, but not chart1, because of the error Error in eval(expr, envir, enclos) : object 'group' not found, can anyone help?
update 01
Thanks for celenius's comment, I tried the following instead:
chart1 <- chart + geom_point(aes(x=hosp.list[,"coord_x"],y=hosp.list[,"coord_y"]))

And this gives me an error below:
Error in data.frame(x = c(-72.4509975, -72.4509975, -92.059305, -73.16584,  :
  arguments imply differing number of rows: 58, 133341

What I am trying to do is as below:

I plot the Connecticut 2010 Census State Legislative District as a layer
I created a dataframe containing several hospitals in Connecticut, geocoded their addresses into long and lat, and put the long/lat into the data frame also (i.e. the dataframe has 3 column)
I want to create a new layer on top of the "district layer", showing the hospitals by dots
the next step I will do somethings like: display name of each dot, highlight the region with hospitals, etc... which should by me future question poster here.

Thanks again.


texts: How to install shp2pgsql and raster2pgsql from source?

texts:
I just installed PostGIS 2.4.0alpha from source by following osgeo instructions on Ubuntu 16.04.

After installation, PostGIS seems to work. However, command line tools that are supposed to come with PostGIS are not installed.

For example, shp2pgsql is missing after sudo make install.

$ shp2pgsql
The program 'shp2pgsql' is currently not installed. You can install it by typing:
sudo apt install postgis

(similarly raster2pgsql is also missing).

Looking into the postgis-2.4.0alpha/loader folder, the executables shp2pgsql and pgsql2shp are present.

What steps do I need to do beyond sudo make install to install the command line tools such as shp2pgsql and raster2pgsql?


texts: I would like to have guidance on how to generate a cost surface for the data I have. In essence, I got a 1x1 DEM for the territory under study (namely, an island) and I have also digitized the spring-line using poly-lines. My aim would be to calculate for any give cell across the whole island the distance from the poly-lines but also taking the terrain elevation into account. I am familiar with Euclidean Distance, but the latter would not be useful since it does not take into account the elevation variable. I have tried to consult ArcGIS's help but I had no luck in finding something that could fit my situation. May be that the solution could possibly boil down to weighting the Euclidean Distance by some value (elevation?), but I do not manage to get my head around the issue.


texts: JS errors on cartoDB map editor - tracker.js

texts:
When editing my map in the CartoDB interface the page breaks with errors - see screenshot.

Any ideas how to resolve this?

tracker.js:14cartodb/table/views/tooltip/custom_templates/infowindow_light not found

infowindow_html_pane.js:202Uncaught TypeError: cdb.templates.getTemplate(...) is not a function(…)


texts: No 'Access-Control-Allow-Origin' header is present error when requesting json response from Geoserver

texts:
I am still working to get my GeoServer WFS displayed in my Leaflet map and cannot get the Ajax response to work properly. I have already run the "getFeature" request in my browser and it returned json as the output. I have also enable jsonp in the web.xml file to avoid cross-origin issues (I'm not using CORS). There appeared to be a problem initially with using dataType: 'jsonp' in my ajax request, and it was recommend that I switch it to json (which I did). Now I am getting the message XMLHttpRequest cannot load http://localhost:8080/geoserver/ows... No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'null' is therefore not allowed access. Seeing as I cannot set the response.addHeader() in GeoServer(at least I don't know how to do that), what is my recourse to solve this issue? Note: it was recommended that I use outputFormat: 'text/javascript' instead of application/json and that I use jsonpCallback: 'parseResponse' instead of getJson. Is that correct? Below is my current code:

var url = 'http://localhost:8080/geoserver/ows';

var defaultParameters = {
  service: 'WFS',
  version: '1.1.0',
  request: 'GetFeature',
  typeName: 'FGCMS:1_Deg_Grid_polygon',
  //maxFeatures: 200,
  outputFormat: 'text/javascript',
  format_options: 'callback:getJson',
  SrsName : 'EPSG:4326'
};

var parameters = L.Util.extend(defaultParameters);
var URL = url + L.Util.getParamString(parameters);

var WFSLayer = null;
var ajax = $.ajax({
  url: URL,
  dataType: 'json',
  jsonpCallback: 'parseResponse',
  //async: false,
  success: WFSLayer,
  error: function() { console.log('Oh snap!'); }
});

function WFSLayer(data) {
  L.geoJson(data, {
    style: function (feature) {
    return {color: 'black'}
    }
  }).addTo(map);
}


texts:
I've been struggling with this problem for awhile.
Here is a curvature raster to illustrate the problem http://i.imgur.com/IIjYAov.jpg

I'm using 10m DEM rasters from the NRCS Data Gateway for a large portion of California. The problem is some strange tiling that appears after merging, but isn't present in the original data (before merging). The grids create low points, so flow accumulation etc is drawn into them.

I've tried mosaic to new raster, create raster catalog and filter in ArcGIS 10.2, merge and build vrt in Qgis2.4 but the problem persists throughout.

When the license frees up, I'm going to try Merge and Filter (production mapping) and see if that works.

Am I just merging files wrong? Are these edge effects from merging tiles?


texts: QGIS GPS Information Menu Automatically Add Points

texts:
We are using a Surface Pro 3 with Windows 8.1 64bit, Bluetooth GPS, and QGIS (2.4 standalone download).  We have opened a new QGIS project and are able to connect to the GPS unit within the GPS Information menu.  The map centers on the current GPS location and track points are automatically added to the screen.  Our problem is that the points are not being saved to the point file.  In fact, the track points on the screen just disappear after closing out of the QGIS project.  Below are the steps that we took.  Can anyone please help?

Added a point vector file
Connect to GPS
Opened an editing session for the point vector file
Within GPS Information menu checked “Automatically save added features” and “Automatically add points”
5.We can see the track develop on the screen as we drive.  When we close the editing session, and close out of QGIS, none of the points have been stored in the shapefile.

In a nutshell, we want to be able to automatically add and save points to a point file during an editing session without having to physically press “add point” each time we need to add a point.


texts:
I have installed Python and Grass6.4 using Macports. When I tried to use "import grass.script" in python environment, it gives "no module found" error. I'd spend much time to fixed it but could not succeed.
Please help me to fix this issue in MacOS Maverick Lion 10.9
When U issued the following commands I get,

which python
/usr/bin/python

which grass64
/opt/local/bin/grass64


texts: MapInfo to QGIS style generator ; error message

texts:
I'm working on Windows 10 with MapInfo Pro 16.0, QGIS 2.16.3.

I tried to use the usefull following plugin called "MapInfo to QGIS style generator":

https://github.com/NathanW2/MapInfo-to-QGIS-style-generator/wiki/Using-MapInfo-to-QGIS-style-generator.

But when i launch this following command via osgeo4w shell :

python mapinfoToQgis.py F:\temp\commerces.TAB F:\temp\out.qml -c ID_SYMBOLOGIE --UseMapInfo

An error message appears. Please reading the following error message:

F:\temp\plugin\plugin> python mapinfoToQgis.py F:\temp\commerces.TAB F:\temp\out.qml -c ID_SYMBOLOGIE --UseMapInfo
Opening MapInfo...
Traceback (most recent call last):
  File "mapinfoToQgis.py", line 308, in <module>
    gen.createQMLFromMapInfoTable(args[0],args[1],columnName)
  File "mapinfoToQgis.py", line 242, in createQMLFromMapInfoTable
    mapinfo = Dispatch("MapInfo.Application")
  File "C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\win32com\client\__init__.py", line 95, in Dispatch
    dispatch, userName = dynamic._GetGoodDispatchAndUserName(dispatch,userName,clsctx)
  File "C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\win32com\client\dynamic.py", line 114, in _GetGoodDispatchAndUserName
    return (_GetGoodDispatch(IDispatch, clsctx), userName)
  File "C:\PROGRA~1\QGIS2~1.16\apps\Python27\lib\site-packages\win32com\client\dynamic.py", line 91, in _GetGoodDispatch
    IDispatch = pythoncom.CoCreateInstance(IDispatch, None, clsctx, pythoncom.IID_IDispatch)
pywintypes.com_error: (-2147221005, 'Chaine de classe incorrecte', None, None)

Could you throw light for me regarding the meaning of the error message and could youhow to help me in order to find a solution?

Sorry but I rejoiced too quickly yesterday.
The Peter 's procedure allows indeed to no longer have an error message during the execution of the command.
However, the plugin does not seem to work. No task is done. I let the process run all the night. And the status stays always on "Opening Mapinfo ... Exporting style table"

...
But nothing else. No qml file has been generated.
It's strange because i managed to use this  plugin in march 2015 last year (with a previous mapinfo version) ... According to you, the reason of the problems could be  from updates realised beteween previous version and last new mapinfo v.16?


texts:
I am trying to utilize the C++ code examples from blog.qgis.org, specifically Tutorial 4, Adding Rasters to QGIS Map Canvas.

I am successfully able to run the example loading a GeoTIFF with QGIS 1.4, however with 1.8 I get "Layer is Not Valid".  As a note, the same GeoTIFF file was able to be loaded in the QGIS 1.8 application just fine, but not in the example code. I've tried with several different files, with the same result. Is there additions to the 1.8 API which require additional calls to get this functionality to work? Problem code below. The full example code can be found at http://blog.qgis.org/node/94

Operating System: Fedora 17 64-bit, QGIS 1.8.0

void MainWindow::addLayer()
{
  QString myFileName = QFileDialog::getOpenFileName(this, tr("Open File"),
                         QCoreApplication::applicationDirPath () + "/data",
                         tr("GeoTiff (*.tif)"));
  QFileInfo myRasterFileInfo(myFileName);
  QgsRasterLayer * mypLayer = new QgsRasterLayer(myRasterFileInfo.filePath(),
  myRasterFileInfo.completeBaseName());
  if (mypLayer->isValid())
  {
     qDebug("Layer is valid");
  }
  else
  {
     qDebug("Layer is NOT valid");
     return;
  }


texts: How to create a custom polygon fill pattern with only one icon?

texts: I have created a custom SVG icon in QGIS. When I zoom in it repeats the shape over and over within the polygon. I want to only see one icon even when I zoom in. For example, the icon is a single diamond, but when I zoom in it fills the poly with twenty diamonds....I just want to see one large diamond.


texts: Pre-loading features in QGIS from geodatabase via XML file

texts:
The background: I've inherited a pretty handy plugin for QGIS and an associated an XML file. The plugin (among other things) when installed, will ensure that each instance of QGIS that you open has a selection of pre-loaded shapefiles - these are saved in various folders on the network drive. (For clarity, these are not loaded in the layer panel. Rather, they're nesting in a newly created 'contents panel' (for lack of a better term). To load them, all one simply has to do is click on the shapefile name in the 'contents panel' and presto, the shapefile will load in the layer panel and the map window.) Below is a screenshot of the 'contents panel':

This all saves time for a lot of the team who are less familiar with QGIS from having to search various network files and folders.

If we receive a new shapefile that would be useful to have pre-loaded, all I have to do is edit the accompanying XML file. (i.e. I don't need to edit any of the plugin code). Screenshot of the XML code header and code that I was intending to use to pre-load GDB data:

As far as I can tell, the syntax and structure are what is specified in the header. (Despite what the screengrab shows, I'd tidied the bracketing and code to match that of the header code so bracketing isn't the issue and I have ensured that the syntax that I've used is correct.)

The issue: Which brings me to the problem that I'm having. I now wish to migrate some of these shapefiles to a geodatabase (I'm trying ESRI's file geodatabase as a proof of concept), however, when I edit the XML file to include for this (or any) geodatabase, the 'contents panel' doesn't update with the new feature class as I'd expect.

What I've tried: I've noticed in the XML file that the FGDB that is referenced is of a raster catalogue type rather than a vector type. I've edited the XML to account for this with no luck.

I've also tried adding the suffix .table' to the desired feature class, but this also doesn't help out.

Lastly, I tried referencing a geodatabase that was on the same network drive as the xml file, still no luck.

I am by no means a code or script person but can tenaciously plug away until something works, but was hoping that someone might have come across this problem before and have a ready solution/ explanation.


texts: How to call the GeoJSON dynamically according to area selection in OpenLayers 3

texts:
I'm trying to get GeoJSON according to the selection of a particular area. For example, if I have selected a built up area it will only show that area . I'm calling GeoJSON data through a PHP page, and I have created the plunker for this (this is the example).

The main problem I'm having is how to assign the vector value from another page to show the GeoJSON.

See this is my JavaScript page where I'm trying to assign the data:

$(document).ready(function(){
       $.getJSON("php/drpdwn.php", success = function(data){
            var items="";
            for(var i = 0; i< data.length; i++){
                items +="<option value='"+data[i].id+"'>" + data[i].name+"</option>";
            }
            $("#a1_title").append(items);
            $("#a1_title").change();
        });
        $("#a1_title").change(function(){
      var selectedValue = $(this).find(":selected").text();
   if( selectedValue == "Landuse/landcover") {

              $.getJSON("php/lulc_db.php",success = function(data){
                  var items="";
                  for(var i = 0; i< data.length; i++){
                      items+="<option value='"+data[i]+"'>" + data[i] +"</option>";
                  }
                  var myNode = document.getElementById("a2_title");
                  while (myNode.firstChild) {
                      myNode.removeChild(myNode.firstChild);
                  }
                  // ADD NEW ITEMS
                  $("#a2_title").append(items);
              });
              $("#a2_title").change(function(){
    $("#save").on('click',function(){
        $.ajax('php/geojson.php?layer=Built Up').then(function(response) {

  var geojsonFormat = new ol.format.GeoJSON();
  var features = geojsonFormat.readFeatures(response,
      {featureProjection: 'EPSG:4326'});
  vector.setSource(features);
});
    });
});
          }

And this is the way I'm adding vector layer:

         var vectorSourcec = new ol.source.Vector({
              format: new ol.format.GeoJSON(),
              });
    var vector = new ol.layer.Vector({
          source:vectorSourcec

          });
var map = new ol.Map({
  layers: [new ol.layer.Tile({
            source: new ol.source.OSM()
          }), vector],
  renderer: 'canvas',
target: document.getElementById('map'),

            view: new ol.View({
            Projection: 'EPSG:4326',
                center: ol.proj.fromLonLat([79.17436, 21.01668]),
  zoom: 5
                                }),
        });

Please improve me here where I'm doing.


texts: Is there a way to convert an RDL filetype to a SHP filetype

texts: I Have an open records file from a county that provided an RDL file extension which we can view in MS 8 Viewer but is there any way of getting that into SHP format so we can normalize the data fields to align with the rest of our dataset?


texts: How to mask and scale raster using rasterio?

texts:
I have been trying to read a raster file in netcdf format which I will later on sample. I need to read using the DatasetReader as float.
When I read:
ds = rasterio.open(f'netcdf:{file}:AOT_L2_Mean', dtype=rasterio.float64)
aot = ds.sample(120.57577514648, 16.003829956055)

The sampled data are still in Int16. These data are aerosol optical thickness thus only have value 0 - 5 and not 0 - 20,000. How to correctly read it as float and have values that are 0 - 5 like as I open it on QGIS?
Here is a sample data: https://drive.google.com/file/d/11ZRiLCIId1G1Dfjlc3TbvdgLOb4B8mit/view?usp=sharing


texts: Slope calculated with QGIS vs slope in GoogleMaps

texts:
I want to calculate the average slope in percent of a vector line in QGIS. I have a DEM layer, converted it to a slope layer and then I use the plugin Profile Tool to calculate the slope of the line.
Is this the correct way? I created a few lines that follow big streets and compared my results with GoogleMaps. For GoogleMaps it's slope in % = elevation-delta / distance. But the results don't match.
Am I just stupid or is google or the plugin wrong?


texts: Populate Integer Field based on Text Field

texts:
Using ArcGIS 10.2.2

In the picture below, I want to populate the Spd_Class field (short integer). The value of the Spd_Class attribute will be based on the TYPE field -

residential = 30
secondary = 35
primary = 50

is there a quick way to use the text field in TYPE to populate Spd_Class? I've been unable to think of a way to do this.


texts: Chirps Monthly mean Data for multiple points or region using GEE

texts:
Using the this Code (https://code.earthengine.google.com/3491eb14052114811c048687b7885bf7) I am trying to obtain monthly mean time series of precipitation data (chirps) of multiple location. I have 2 queries

how can I get time series in same graph for these six location? I tried using ui.Chart.image.seriesByRegion but I am getting error stating - Error generating chart: No features contain non-null values of "system:time_start"

How can I get Single Mean value of the month of a year for each of these.

'''
var r = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(77.15067,9.34339), {'label': 'Anathode Dam'}),
  ee.Feature(ee.Geometry.Point(77.06994606,8.954271389), {'label': 'Kallada dam'}),
  ee.Feature(ee.Geometry.Point(76.74478393,9.852621715), {'label': 'Mallankar dam'}),
  ee.Feature(ee.Geometry.Point(75.82326883,11.59669815), {'label': 'Kuttiyadi Dam'}),
  ee.Feature(ee.Geometry.Point(76.53852053,10.98574468), {'label': 'kanhira Dam'}),
  ee.Feature(ee.Geometry.Point(76.79667658,10.62448299), {'label': 'Meenkara Dam'}),
  ee.Feature(ee.Geometry.Point(77.14592432,8.535029363), {'label': 'Neyyar dam'})
]);
Map.addLayer(r);

var precip1year=chirps.filterDate('2018-01-01', '2018-12-31');
var TS2 = ui.Chart.image.series(
  precip1year, r,  ee.Reducer.mean(), 1000, 'system:time_start')
  .setOptions({
        title: 'Monthly precipitation',
        hAxis: {title: 'Time'},
        vAxis: {title: 'Precipitation (mm)'}
  });

print(TS2)

var precip1year=chirps.filterDate('2017-01-01', '2019-01-01');
var TS2 = ui.Chart.image.seriesByRegion(
  precip1year, r,  ee.Reducer.mean(),'precipitation', 1000, 'system:time_start', 'label')
  .setChartType('ScatterChart')
  .setOptions({
        title: 'Monthly precipitation',
        hAxis: {title: 'Time'},
        vAxis: {title: 'Precipitation (mm)'},
        lineWidth: 1,
        pointSize: 4,
          series: {
           0: {color: 'FF0000'},
           1: {color: '00FF00'},
           2: {color: '0000FF'},
            3: {color: '#27c8ff'},
            4: {color: 'Orange'},
            5: {color: 'Purple'},
            6: {color: '#ff58fe'},
}});
print(TS2);


texts: SyntaxError: invalid syntax () in Pl/Python (PostgreSQL)

texts:
I have python 2.7 and I have weather script which returns temperature infos, I would like to implement this script into PostgreSQL. I always got this error: DETAIL:  SyntaxError: invalid syntax (<string>, line 10)

code:

    CREATE OR REPLACE FUNCTION GetWeather(lon float, lat float)
    RETURNS float
    AS $$
    import urllib2
    import simplejson as json

    data = urllib2.urlopen(
    "http://api.openweathermap.org/data/2.1/find/station?lat=%s&lon=%s&cnt=1"% (lat, lon))
    js_data = json.load(data)
    if js_data['cod'] == '200':
    if int(js_data['cnt'])>0:
    station = js_data['list'][0]
    print 'Data from weather station %s' %station['name']
    if 'main' in station:
    if 'temp' in station['main']:
    temperature = station['main']['temp'] - 273.15
    else:temperature = None
    else:temperature = None

    return temperature

    $$ LANGUAGE plpythonu;

I tried this version too, it is working in pure python but not in postgresql

CREATE OR REPLACE FUNCTION GetWeather(lon float, lat float)
    RETURNS float
    AS $$
    import urllib2
    import simplejson as json

    def get_temp(lat, lon):
    data = urllib2.urlopen(
    "http://api.openweathermap.org/data/2.1/find/station?lat=%s&lon=%s&cnt=1"% (lat, lon))
    js_data = json.load(data)
    try:
    return js_data['list'][0]['main']['temp']
    except (KeyError, IndexError):
        return None

    $$ LANGUAGE plpythonu;


texts:
I have several thousand markers stored in a kmz file that I load into a Google Maps (v3 API).

I've looked around to find a way to cluster the overlay data, but to no avail.
MarkerClusterer looks nice, but it only seems to work with a marker array?

function initialize() {
    geocoder = new google.maps.Geocoder();
    var latlng = new google.maps.LatLng(25,-20);
    var myOptions = {
        zoom: 2,
        center: latlng,
        mapTypeId: google.maps.MapTypeId.ROADMAP
    };
    map = new google.maps.Map(document.getElementById("map_canvas"), myOptions);
    var ctaLayer = new google.maps.KmlLayer('thedata.kmz');
    ctaLayer.setMap(map);

    //cluster
    var mcOptions = {gridSize: 50, maxZoom: 15};
    var mc = new MarkerClusterer(map, ctaLayer, mcOptions);
}

How can I cluster the markers in the overlay, just like markerclusterer does it?


texts: Creating grid to divide large map into smaller sections using ArcGIS for Desktop?

texts:
I have a map of small remnants of forest along creeks for an entire county in California. I need to divide the whole map into equal sections so I can print it so you can actually see the sections of forest. I'd like a map with a grid which I can look at, and then look at a zoomed in map of (for example) section A1.

I've tried the Grids and Graticules wizard, which gives me a grid for the map, but when I zoom, the grid remains the same. I want the grid sections to be georeferenced with the map, so when I zoom in, the grid zooms also.

I'm using ArcMap 10.1.


texts: Problems with .ecw files using gdalbuildvrt

texts:
I used

gdalbuildvrt C:\Users\Vector_district.vrt D:\177_SwissImage_10cm_2017_190410\2548\*.ecw

trying to create a vrt.
The console gives back following error for each file in the folder:

  ERROR 4:
  `D:\177_SwissImage_10cm_2017_190410\2548\file.ecw'
  not recognized as a supported file format. Warning 1: Can't open
  D:\177_SwissImage_10cm_2017_190410\2548\file.ecw.
  Skipping it

When I use an other folder with .tif files it works perfectly.

How do I have to adapt my code to create a vrt with .ecw-files?


texts: Export KML from GeoJSON in QGIS and coordinates change?

texts:
I'm very new in this area. I have a geojson vector layer file, with coordinates like

{ "type": "Feature",  "geometry": { "type": "Point", "coordinates": [ 487677.2993499999866, 112018.26101 ] } },

I open it in QGIS, I change the default SRC of the couche and from the project to match with my data (EPSG:4149).

When I export then this in KML, I have those coordinates for the same point:

  <Placemark>
          <Point><coordinates>-122.700650000013411,112018.26101</coordinates></Point>
  </Placemark>

I don't understand why the first coordinate changes. Can you help me?


texts: ST_Within between 2 Multilinestrings not matching in PostGIS/QGIS

texts:
I'm working with QGIS 3.4 LTR.
I have 2 geometric tables with Multilinestring type geometry in the same SRID : Table1 (thin black lines with thin black labels) and Table2 (thick colored lines with buffer zone labels).

I want to get the id of Table2 features in Table1 column named id_table1
 when : geometries of Table1 are within geometries of Table2 OR geometries of Table1 contain geometries of Table2 OR when the geometries of the 2 tables features overlaps.

I have this query :

SELECT table1.id AS "id_table1",
string_agg(table2.id::text,',') AS "id_table2"
FROM schema.table1
LEFT JOIN schema.table2 ON ST_WITHIN (table1.geom,table2.geom) OR ST_CONTAINS (table1.geom,table2.geom) OR ST_OVERLAPS (table1.geom,table2.geom)
GROUP BY table1.id

In the table extract, we can see that id_table1 80, 83 and 87 did not get id_table2 :


texts: I have a Web Application (ASP.NET) made in ArcGIS Server Manager. How to create a button for the user to export his vision for PDF / Image?

texts:
I am using ArcGIS Server Manager to create Web Applications (ASP.NET).

I would like the customer to access has an option to export your vision to PDF or Image, or both.

I would like to generate the file, keep up the zoom, pan, measurements and other iterations that the user made. It would be like a "Print Screen" on your computer, and generates the PDF/Picture file.

I thought about using some javascript API, but not found.

I thought about using Model Builder, but not tried it yet because I think it will not keep the iterations of the user did on the screen: zoom, pan, etc..

I have ArcGIS 10, with SP5.

Does anyone have a clue how to do this export functionality PDF / Image file on a Web Application?

Update:

I'll try to do this:

www.codeforge.com/read/132008/ExportMap.cs__html

In my Web Application standard:

    \wwwroot\ArcGIS\Manager\Modules\Applications\Templates\mv_cs

Now I'm fixing my computer to to do this. When I do then I'll share with you. ;)


texts:
I tried using the get_googlemap function but I have an error that I cannot fix.
I was wondering if anyone could help me in this regard.
ggmap::register_google(key = "xxxxxxxxxxxxx (my API Key)")
d <- function(x=-95.36, y=29.76, n,r,a){
  round(data.frame(
    lon = jitter(rep(x,n), amount = a),
    lat = jitter(rep(y,n), amount = a)
  ), digits = r)
}
df <- d(n=50,r=3,a=.3)
map <- get_googlemap(markers = df, path = df, scale = 2)
ggmap(map)

It is an example provided at https://www.rdocumentation.org/packages/ggmap/versions/3.0.0/topics/get_googlemap
but I encounter this error every time I run this code

Error in aperm.default(map, c(2, 1, 3)) :
invalid first argument, must be an array
In addition: Warning message:
In get_googlemap(markers = df, path = df, scale = 2) : HTTP 400 Bad Reques

even using this code results in the same way.
map <- get_googlemap(c(-97.14667, 31.5493))


texts: Rasterizing vector shapefiles using multiple burn in values

texts:
I have a vector shapefile of a geographic area that contains several different object codes that I need to rasterize all as one image.

I am working with QGIS, but the Raster>conversion tool seems to be unable to do this (see below).  In the tool, it seems as though you cannot select more than one value, not to mention the tool - even when it works - appears to ignore whatever value you have set and rasterizes all values at once, in the same image.  That would be ok, but I have a large volume of output to complete, and I don't want to separate every needed raster mask into a separate layer, as this would be extremely inefficient and time consuming

So, I tried running the GDAL command it generates, and even the GDAL command, when entered in the python console, errors out.  I have never worked with GDAL before and I am way out of my depth, so I am looking for pointers.

Below is the command I started with (and that QGIS ran successfully when I attempted to rasterize using the form), and no variation of this seems to work in the console, and any edit I perform also returns a syntax error.

  gdal_rasterize -l bornholm_final_OSM -a code -ts 100492.0 100492.0 -te
  14.510849428030344 54.922004004836744 15.299793055831714 55.37191520877061 -ot Float32 -of GTiff "E:/My Documents/OSM Stuffbornholm_final_OSM.shp"
  "E:/My Documents/OSM Stuff/point-5m test_qgis2.tif"

I have looked at the GDAL documentation about "gdal_rasterize", but I am truly stumped at this point.  I just need the raster to do multiple classes in one image.

Using @BenW'S answer I am now able to create a raster using GDAL, but I now need to ensure that it retains the colors.


texts: How to add a progress bar to an ogr2ogr script

texts:
We have started using ogr2ogr to upload MapInfo tables into Microsoft SQL server in my workplace. Some of the tables are massive (>1,20,000 records) and there is no way of knowing how far the script has progressed without a status bar.

I was wondering if there is a command that can be used within the ogr2ogr script which actually tells us how much more of the processing is still left/pending (e.g. 50% complete etc.).


texts: Executing Styles based on rule Geoserver

texts:
I would like to execute a style based on some parameter passed through WMS Rest request. What I am trying to achieve is, while requesting the tile, I want a particular country to be filled in separate color. So I would pass "USA" in ENV parameter and the tiles for USA should be colored in different color than the rest of the world. So in the final output I would have the whole world map but only USA highlighted.

Regards
Sudheesh


texts:
What i've got is a map served by the tilecache and shown by openlayers and geoext.

Now i need to print this map, so what i've done is:

made a button which creates (window.open()) a new window with my map in it and window.print() on body onload.

This works well exept that i'm getting a second blank page. So the question is: how (and where) do i set map size and proportions to have it printed correctly?

map is created this way:

map = new OpenLayers.Map("map", {
    maxExtent : new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34),
    maxResolution : 156543.03125 / 64,
    numZoomLevels : 11,
    units : 'm',
    projection : "EPSG:900913",
    displayProjection : new OpenLayers.Projection("EPSG:4326"),
});


texts:
I want to use the remicres/otb-mosaic plugin (https://github.com/remicres/otb-mosaic) in QGIS.

Fig. 1 is the usage introduction given on the webpage.

However I cannot find the function in QGIS after I copied the file "Mosaic.xml" into the QGIS installation path as instructed on the webpage.

How to find the function in QGIS?


texts: Using a virtual layer, display polygons based on location of another polygon layer

texts:
QGIS 3.4

I have a .gml file of every freehold title in a city outlined red below. I want to create a virtual layer of the freeholds which are contained within the blue-green isochrones. Even if a small portion of the freehold title intersects the isochrone, I need the entire freehold to be displayed.

my code so far:

SELECT *
FROM Freeholds, Isochrones
WHERE st_within(Freeholds.geometry, Isochrones.geometry) = 1

I am a novice at SQL and recieve the error 'Referenced table Freeholds in query not found!' How can I correct the query?


texts: Numpy Array to ArcGIS GRID gives dimensional error with masked array

texts:
I have an 84x73 masked numpy array to load into ArcGIS (version 10.3.1, Windows 8.1).

To overcome the fact that Masked Numpy arrays cannot be saved to .npy file as of yet, I am saving my array from Jupyter as

import numpy

myarray.dump("C:\MyName\MyFolder\Array_to_raster.npy")

Then, I load this array in ArcGIS, as suggested at Numpy Array to ArcGIS GRID gives ValueError?:

myarray = numpy.load(r"C:\MyName\MyFolder\Array_to_raster.npy")

And then I attempt to create a raster out of it as

arcpy.NumPyArrayToRaster (myarray, arcpy.Point(447328.4200,357465.6840), 1500, 1500, -2147483648)

where -2147483648 is the fill_value of my masked cells.

However, this throws an error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 2292, in NumPyArrayToRaster
    return _NumPyArrayToRaster(*args, **kwargs)
ValueError: Argument in_array: A two or three dimensional NumPy array is required.

But as I said, I already have a 2D array:

>>> myarray.shape
(84, 73)

How can I create the raster out of myarray?


texts: How to use variable in expression in arcpy CalculateField tool?

texts:
I want to use the following expression in the arcpy.CalculateField_management tool.

I have a variable called "value" and a field called "Field01"

I want to use (value - Field01) / 20

I can´t seem to find the correct syntax with arcpy.CalculateField_management.

arcpy.CalculateField_management (fc, "test",expression, "PYTHON")

Given error is: TypeError: cannot concatenate 'str' and 'int' objects

I have tried the following and many other combinations:

expression = value - "!XKoordInt!" / 20

expression = "("+value+"-!XKoordInt!)/20")"


texts: Any way to create minimal example shapefile with fake data in R?

texts:
I am wondering how to create minimal example shapefile of certain countries such as US, Germany, or Brazil with fake data. I am new with geospatial data and it is nice to know how to create minimal geospatial data in R. I need reproducible shapefile because I want to create my minimal solution on it and test my approach.

How can I make this happen minimal shapefile of certain countries with fake data in R? Any quick approach to do it? Any idea?

Update:

I am going to create shapefile of certain countries such as US, Germany or vice-versa by state or county level but fill it with fake random data for reproducible use. Any better idea?


texts: Error in GRASS using i.segment.uspo function

texts:
So I was using this function called i.segment.uspo and I got this error

  ERROR: Insufficient number of non-NULL cells in current region
  Traceback (most recent call last):
    File "/usr/lib/grass72/scripts/i.segment.uspo", line 935,
  in
      main()
    File "/usr/lib/grass72/scripts/i.segment.uspo", line 838,
  in main
      for mapname, lv, autocor, threshold, minsize in
  iter(result_queue.get, 'STOP'):
  ValueError: need more than 1 value to unpack

now both of my raster and vector are in the same coordinate reference system which is 2100. My vector data are polygons and where created in QGIS so they could represent different land-covers and my raster was clipped, again, in QGIS so that NULL values would be removed. Is there a way to fix that error?

EDIT

the code:

CATS="$(v.db.select -c map=train_vec columns=cat)"
for i in $CATS
do
where="cat=$i"
regionname="subset_uspo_gera_$i"
v.extract --overwrite input=train_vec type=area where=$where output=$regionname
g.region --overwrite vector=$regionname save=region_uspo_$i
step_gera=$step_gera"region_uspo_gera_corine_$i,"

i.segment.uspo --overwrite group=OBIAgera output=/home/user/Desktop/res/obiatest_$i.csv segment_map=best_gera_$i region=region_uspo_$i threshold_start=0.601 threshold_stop=0.901 threshold_step=0.001 minsizes=8 memory=4000 processes=4

r.to.vect --overwrite input=best_vec_"$i"_region_uspo_vec_corine_"$i"_rank1 output=temp_best_$i type=area
v.out.ogr --overwrite input=temp_best_$i type=area output=/home/user/Desktop/resaults/res_$i.shp
done

The error most likely comes from i.segment.uspo but I don't know if it's because of my training areas (basically the regions) or because of something else

P.S the code was based on a python algorithm which I decided to do it in shell script for my thesis.

EDIT 2
Here is the source code of the command and it stops somewhere around here:

837         if rg:
838             for mapname, lv, autocor, threshold, minsize in iter(result_queue.get, 'STOP'):
839                 regional_maplist.append(mapname)
840                 variancelist.append(lv)
841                 autocorlist.append(autocor)
842                 threshlist.append(threshold)
843                 minsizelist.append(minsize)


texts: Grass cmd file blocked by group strategy

texts:
I encounter difficulties whenever I try to run any GRASS plugin from within QGIS (mainly) : it looks like my company has set up a group strategy preventing the call to .cmd files from programms.
This is only empiric guess since I am not really skilled with such strategies and I could be indeed wrong about the cause matching the symptom!

Anyway : here's the message grom the GRASS plugin (translated from french to english for a better understanding):

  Execution of
  C:\Users\ME\AppData\Local\Temp\processing_3eb0b5419fcc4fcfb2518bc4d51d5a42\grassdata\grass_batch_job.cmd
  ...  This programme is blocked by a group stategy. For further
  information, contact your system administrator.

And I am quite upset about it since GRASS holds very usefull tools..

Does anyone have a workaround to have GRASS work in this penitentiary environnement?


texts:
Possible Duplicate:
Converting between KML and shapefile (SHP) format

I have SHP files, I want KML files, what tools can convert?


texts: Python IndentationError when trying to edit graphical modeler in QGis. Can it be fixed?

texts:
I am trying to edit a model I built in the QGis graphical modeler and I get the following error.

IndentationError: unindent does not match any outer indentation level

2017-11-27T14:53:42 1   Traceback (most recent call last):
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\modeler\EditModelAction.py", line 42, in execute
                dlg = ModelerDialog(self.alg.getCopy())
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\modeler\ModelerDialog.py", line 190, in __init__
                self.repaintModel()
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\modeler\ModelerDialog.py", line 374, in repaintModel
                self.scene.paintModel(self.alg)
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\modeler\ModelerScene.py", line 115, in paintModel
                arrow = ModelerArrowItem(self.algItems[depend], -1,
            KeyError: u'SCRIPTADDFEATUREWITHOUTGEOMETRY_1'

I added this script to work around an extract by location error and when i tried to remove it Qgis closed the modeler and threw a warning.  Now I only get this error when trying to edit.  It's a pretty complicated model with 19 processes and I'd hate to have to go back and rebuild it.  Is there a way to fix this?


texts: How to use feature collection after esri leaflet query.run()?

texts:
35 years coding but new to GIS leaflet and Esri.leaflet.  Javascript I can do but it is not one my long term primary languages.  I am trying to make a query result available for processing after some dependencies are taken care of.  cQuery here is pre-defined using L.esri.Query using layer 14 of the service it is pointed at.
var rslt;
cQuery.run(function (error, fset) {
    // if (error) display alert, return
    rslt = fset;

    if (fset.features.length = 0 && check2ndLayer) {
        cQuery.layer(15).run(function (error, fset2) {
            // same thing here for error
            rslt = fset2;
        });
    };

    // A
});

// B

The code runs but...inserting some alerts I have determined layer 14 query completes, condition to run query for 15 engages.  Code placed at "A" can process "fset" or "rslt" and does so if layer 14 found something.  If not, code runs past there and then I get the indicator that the layer 15 query has found something to process but it never does.  What I want is to end up with one featureCollection in "rslt" and process that at "B" but no matter what I do "rslt" seems to always be undefined at that point.
What am I doing incorrectly?


texts:
i work with EVI MODIS and i need to stack all raster i had using python but i don't know how

i'm using this backag and script to read my rasters

    import os
    import glob
    from osgeo import gdal
    import numpy as np
    import matplotlib.pyplot as plt
    import scipy.ndimage
    import pandas as pd
    import datetime as dt
    import rasterio as rio

    #repertoire de travail

    inDir = r'D:\MODIS\Doy2000/'
    os.chdir(inDir)
    #repertoir sorti

    outDir = os.path.normpath(os.path.split(inDir)[0] + os.sep + 'output') + '\\'
    if not os.path.exists(outDir): os.makedirs(outDir)

    # Create list of MODIS rasters for surface reflectance
    modis_EVI = glob.glob('MOD13Q1.006__250m_16_days_EVI_**.tif')

modis_EVI
Out[32]:
['MOD13Q1.006__250m_16_days_EVI_doy2000209_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000225_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000241_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000257_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000273_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000289_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000305_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000321_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000337_aid0001.tif',
 'MOD13Q1.006__250m_16_days_EVI_doy2000353_aid0001.tif']

but idon't know how stack all my raster in python
need help i'm very new in python


texts: Find cluster of points within a 14K dataset which are laying within a circle of a certain radius

texts:
I have a dataset of geolocated points, which are dispersed over a large area
I need to find clusters of a minimum number of X points within these sets which are laying within a circle of which the radius is R
I want to find those clusters, retreive the orginal data of the point in the cluster but also the geographical midpoint of the point within that cluster


texts: How to determine mean VIIRS Nighttime Day/Night Band Composites Version 1 for each district (county)?

texts:
How to determine mean VIIRS Nighttime Lights for each district (county) shapefile attached ?
The exported csv file does not contain required results!
var start = ee.Date('2017-01-01');
var end = ee.Date('2017-12-31');
var night_img = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG')
.filterDate(start, end)
.select('avg_rad')
.first();


texts:
When I try the geocoding function in QGIS plugin mmqgis it doesn´t recognize my columns headers and when I try to fill in the fields, it only recognizes the whole

The CSV has around 200 entries.

Does someone know a workaround?


texts: Creating tool in Python toolbox to export field list to Excel?

texts:
I am creating my first tool in a Python Toolbox (*.pyt) to export a list of fieldnames to Excel. I have the following bits of code:

import arcpy

class Toolbox(object):
    def __init__(self):

        self.label = "Toolbox"
        self.alias = "Field to Table"
        self.tools = [Tool]

class Tool(object):
    def __init__(self):

        self.label = "Tool"
        self.description = "Export Fields to Table"
        self.canRunInBackground = False

    def getParameterInfo(self):
        param0 = arcpy.Parameter(
        displayName="Shapefile",
        name="shapefile",
        datatype="DEShapefile",
        parameterType="Required",
        direction="Input")

        param1 = arcpy.Parameter(
        displayName="Field Table",
        name="datafile",
        datatype="DETable",
        parameterType="Required",
        direction="Output")

        params = [param0, param1]
        return params

    def isLicensed(self):
        return True

    def updateParameters(self, parameters):
        return

    def updateMessages(self, parameters):
        return

    def execute(self, parameters, messages):
        shapefile = params[0].valueAsText
        fields = arcpy.ListFields(shapefile)
        datafile = params[1].valueAsText
        spreds = open (datafile, 'w')
        for f in fields:
            spreds.write(f.name + "\n")

        spreds.close()
        return

Firstly, the actual execution part of the code works fine as I have tested it. I'm not sure what's wrong when I convert to a .pyt file.

Once I have got this working, is there anyway to create a parameter within the tool to create a new excel/dbf file, and run the script with that information, rather than having to create one first, before running the script. I hope that makes sense.


texts: Seeking ArcGIS Pro tool to get min and max latitude and longitude coordinates?

texts:
I am very new to ArcGIS Pro and have been struggling with getting the minimum and maximum latitude/longitude of an area in ArcGIS Pro. I am looking to get four coordinates of a rectangle that encompasses each polygon I have in a shapefile.

I have tried calculating extents of the minimal bounding geometry rectangles and vertices to points, so I thought I would reach out here. This is a description of what I have done so far.

I have about 2800 polygons in a shape file that represent conservation areas across the US. I used the minimal bounding geometry method of arcpy (script below) to get rectangles that encompass each polygon:

arcpy.MinimumBoundingGeometry_management('US_IBAs_2018Feb',"G:\IBA_Update\GIS\Bounding\Output", "RECTANGLE_BY_AREA", None)

I then tried to calculate the extent geometry attributes of each rectangle. This seemed to work, but I could not figure out why the x/y min/max weren't being displayed in decimal degrees instead of meters (see photo).

arcpy.AddGeometryAttributes_management('Output', 'EXTENT')

I tried making sure the map and the map and what is being returned from the extent were both in WGS84 or Web Mercator Aux Sphere, but had the same problem each time.

So from here I decided to just try and use the vertices to point tool. This gave me four coordinate points of each rectangle which is exactly what I need! However, all of the points have the same coordinate.

How can I correct this issue?


texts: find location details within circle buffer

texts:
Wondering how would i go about getting this to work.

I have a table with location data (lat, long, homevalue) in sql server 2008. Given a lat long say (32.113, -81.3225), i want to draw a 50 mile radius circle and get the number of locations and total home value within the circle buffer.

Thanks
vic

Mapoholic,

Thanks for your reply. If i understand you correctly, you are using the "shape" as an attribute to perform stintersect and stbuffer on the location table. My limitations are that, i have the location table which has just the lat, long, homevalue in dollars. I need to independently draw a cirlce based on the given lat long, which if i am doing it right would be like this

DECLARE @Result geography
SELECT @Result = geography::Point(32.113, -81.3225,4326);
select @Result.STBuffer(10000);

and then i need to JOIN or use some method to check if the lat longs in the location table are within the @Result.STBuffer(10000). Any ideas?

Thanks again

vic


texts:
I'm new using QGis. I read somewhere that QMS was better and more stable than the OpenLayers complement to use Google Maps. My problem is that it does not appear when I search with the plugins administrator. I use Ubuntu 14.04 and Qgis 2.0.

regards


texts: What is the real distance between positions?

texts:
I'm trying to figure out what gives me the most accurate distance between positions:

Measuring distance directly (using Geod.inv)?
Measuring distance in 2D (after transforming to a known projected coordinate reference system)?

I'm trying to project latitude/longitude positions to 2D because it is easier to work with them that way, I'm trying to work with centimeter accuracy.

As a sanity check I did some tests to check if distance between positions were preserved.

Here is an interesting data point
(note: I'm actually using DotSpatial library but I re-did this example in pyproj because it is a better-known tool):

Measuring distance "directly" (using geod.inv):

>>> from pyproj import Geod
>>> wgs84_geod = Geod(ellps='WGS84')
>>> lat1, lon1 = (27.23150120120691, -91.522575364888709) # somewhere in the gulf of mexico
>>> lat2, lon2 = (27.239416154284083, -91.520817319790083) # near lat1, lon1
>>> az12,az21,dist = wgs84_geod.inv(lon1,lat1,lon2,lat2)
>>> dist
894.1531771043997

Measuring distance in 2D (after projecting):

>>> from pyproj import Proj
>>> wgs84 = Proj('+proj=longlat +datum=WGS84 +no_defs') # http://epsg.io/4326
>>> nad27Blm15n = Proj('+proj=tmerc +lat_0=0 +lon_0=-93 +k=0.9996 +x_0=500000.001016002 +y_0=0 +datum=NAD27 +no_defs') # http://epsg.io/32065 but in meters
>>> from pyproj import transform
>>> p1_projected = transform(wgs84, nad27Blm15n, lon1, lat1)
>>> p2_projected = transform(wgs84, nad27Blm15n, lon2, lat2)
>>> p1_projected
(646301.4612381, 3012741.116573205)
>>> p2_projected
(646465.2249033558, 3013620.015186602)
>>> import math
>>> math.sqrt((p1_projected[0] - p2_projected[0])**2 + (p1_projected[1] - p2_projected[1])**2)
894.0253411891856

So, I'm getting 894.1531771043997 and 894.0253411891856 respectively which have about 0.128 meters difference.

Which one is the most correct distance? Why?

(note: I also have examples where the differences are > 1m, but only when the distances are > 2km).


texts: Remove a part of line features shapefile using a polygon shapefile in QGIS

texts:
I have a road network map which I need to cut off parts off using a polygon.
Both road network and the polygon are shapefiles. What I need to do is to use the polygon like a cookie cutter and cut off the road parts which fall into the polygon. After the operation is complete, what should be left is the road network which does not belong inside those polygons
So start:

Polygon:

Expected Result:

I tried to do this in QGIS using the Difference, intersection and other features in the Geoprocessing toolbox but they don't seem to work because I have two types of features. (Lines and Polygon)
I cant convert the road layer to polygon, that would end up in a disastrous mesh of polygons.
Is there a way to do this in QGIS? I am very new to QGIS, so assume a complete dummy.


texts: Calculate area of overlap within polygons

texts:
ArcGIS Pro 10.6

I have a polygon layer with a column of 5 unique values, class names.
I have a second polygon layer with a column of unique values, Parcel number.

Each parcel will have multiple classes within its bounds.

I need a layer where each parcel appears once and is associated with the area of each class within its bounds.

I have tried pivots, intersects, merges, tabulations... the list goes on. Somebody save my sanity please.


texts:
Possible Duplicate:
  Cannot populate a Ext GridPanel through GeoExt WMSCapabilitiesStore

I'm using Mapserver + GeoExt, but i can't not recive any data from this request, any one know how to use that?

var grid = new Ext.grid.GridPanel({
                    title: "Availables WMS Layers",
                    renderTo: document.body,
                    width: 800,
                    height: 150,
                    viewConfig: {forceFit: true},
                    store: new GeoExt.data.WMSCapabilitiesStore({
                        url: "http://localhost:8085/fcgi-bin/wms.exe?SERVICE=WMS&REQUEST=GetCapabilities&VERSION=1.1.1",
                        autoLoad: true
                    }),
                    columns: [
                        {header: "Name", dataIndex: "name", sortable: true}
                    ]
                });


texts: MapInfo Calculating the number of Points that Fall within a Boundary

texts:
I have two layers in MapInfo (points and polygons). I need to select points which fall within polygons where the income > 1500

I tryed this one but not working: http://www.sgsi.com/MIUserGroup/imgs/sql1.jpg

I also tryed this SQL code but doesn't work:
points.Obj Within polygon.Obj where income > 1500

Any idea?


texts: Routing algorithm for generating isochrones using open route service

texts: I am looking for a suitable tool for generating isochrones. I need to know the routing algorithm which is used by the tool that I will generate the isochrones using it. So I would like to know what is the routing algorithm that open route service uses to generate isochrones?


texts: locating latitude-longitude point in a piste-map with matrix transforms

texts:
I am trying to plot points of a particular latitude longitude on this image of a Swiss mountainous landscape.

What I am doing is I have some 54 known location-points plotted on this image (not visible in the linked image), which is a piste map, and is therefore an artist's impression of the landscape. I found out the latitude-longitudes and also height above sea level for these 54 known points. I opened my piste map in the GIMP, and plotted the 54 points, while at the same time reading their co-ordinates in points, in the piste-map. So, I have x-y coordinates (lat-long) and z coordinate (height above sea level) for all those 54 points (so its a 3d co-ordinate system that I have), and I also have X-Y coordinates of those places in my piste-map (a 2d co-ordinate system), and I assume a Z=1 for all those points.

To locate a particular lat-long on my piste-map(I have the coordinates for this place, i.e. the lat-long but not elevation), I need to calculate their X-Y on the piste-map. My approach is, that I calculate distances between this point which is to be located and those 54 known points , and then I sort them in ascending order. The first two points in this order will be the nearest points to the one which we want to locate. Now, I find my 3rd known point, and for that, I loop over the rest of the points, to get the nearest possible point such that this point forms a triangle with the first two closest points with the point which is to be located falling inside the formed triangle. Now, I have three known points forming a triangle say triangle ABC. For these known points, I have a triangle A'B'C' in the piste map, and the point to be located lies in the latter. I figure out the transform matrix which when multiplied with the 3X3 matrix of the three known points(in this 3X3 matrix, each row corresponds to a point, with the column values x,y and z in the 3d scape) results in a 3X3 matrix of those points in the piste-map. I multiply this transform matrix calculated with a new 3X3 matrix for three points, the first one being the point to be located, and the remaining two are the last two known points in the first 3X3 matrix. This way, the first two elements of my resulting matrix are the X and Y values for the point to be located, and I proceed ahead to draw the point on the piste-map. Is this approach correct with respect to the way I am using the transform matrix? Thanks! I would really appreciate if any alternative approaches are suggested as well!


texts: how to geowebcache my own group layer in geoserver?

texts:
I am unable to geowebcache my layer group in geoserver. But in Geowebcache demo page, default layer groups are available. How to add my group layer into geowebcache?
please guide me.


texts: How are the coordinates for multipolygons arranged?

texts:
How are the coordinates for a multipolygon (i.e. polygon with a hole in it) arranged?

Do shapefiles and geoJSON files always save the coordinates for the outer ring first?

I'm asking because I have a set of multipolygons with holes and I'm looking to convert them to standard polygons by removing the coordinates of the holes.

If the outer ring of the polygon is always saved as the first set of coordinates, then I know I can easily select those.


texts: I want to make an interactive map where the end user can filter the data

texts:
I've been looking at services like CartoDB and Mapbox for making and hosting some interactive maps. However, it doesn't seem like they provide the option for the end user to filter through the data themselves. Does you know of any out-of-the-box solutions for creating this kind of functionality? Or perhaps CartoDB ad Mapbox do this and I just don't know how.

When I talk about filtering, I mean this. Say a map shows soccer fields around the USA. The user can then toggle between different parts of the data and see only what they're interested in (outdoor fields, indoor fields, grassy fields, turf fields, etc).


texts: Grouping connected linestrings in PostGIS?

texts:
I have a table of streets I've selected based on a set of attributes (let's say it's speed_limit < 25).  There are groups of streets that are locally contiguous;  I'd like to group these sets of connected linestrings into GeometryCollections.  In the image below, there would be two GeometryCollections: one with the red lines and one with the blue lines.

I tried running a couple of "dissolve, deaggregate" queries along the lines of:

SELECT (ST_Dump(st_union)).geom
FROM
    (SELECT ST_Union(geom) FROM roads) sq

With everything I've tried, I either end up with a single feature (ST_Union) or my original geometry (ST_Dump of ST_Union).

Maybe it's possible to do this with some kind of WITH RECURSIVE magic?


texts: Invalid argument specified while applying a function on an image bands

texts:
I tried standard scaling of each band of an image as follows:
var aoi = ee.Geometry.Polygon(
        [[[-98.3640625, 36.61967042308436],
          [-98.3640625, 32.86546781020435],
          [-92.91484375000002, 32.86546781020435],
          [-92.91484375000002, 36.61967042308436]]], null, false);
Map.addLayer(aoi);

var S2 = ee.ImageCollection('COPERNICUS/S2')
                    .filterBounds(aoi)
                    .filterDate('2019-07-01','2019-07-02')
                    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',90))
                    .select('B.*').toBands();
print (S2);

var normalize = function(image) {
  var names = image.bandNames.map(function(name){
    var col = S2.select(name);
    var out = col.subtract(col.mean()).divide(col.stdDev())});
    return out.toBands()};

var normalized = ee.ImageCollection.fromImages(normalize).toBands();
print (normalized);

However, I got the following error.
Invalid argument specified for ee.List(): function normalize(_0) {
return newF.apply(this, arguments);
}
Can you help me?


texts: My raster does not line up with rest of vector shapefiles, even though both rely on NAD 83

texts:
I started creating a map in ArcMap by adding vector shapefiles. All of my vectors have data sources as:

Data Type:  Shapefile Feature Class
Geometry Type:  Line

Projected Coordinate System:    NAD_1983_StatePlane_North_Carolina_FIPS_3200
Projection: Lambert_Conformal_Conic
False_Easting:  609601.22000000
False_Northing: 0.00000000
Central_Meridian:   -79.00000000
Standard_Parallel_1:    34.33333333
Standard_Parallel_2:    36.16666667
Latitude_Of_Origin: 33.75000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

When I bring in a raster image downloaded from a North Carolina imagery earth-explorer database, the raster image looks a few hundred miles off from where it should be. The data source of the raster images is:

Spatial Reference: NAD_1983_2011_StatePlane_North_Carolina_FIPS_3200_Ft_US
Linear Unit: Foot_US (0.304801)
Angular Unit: Degree (0.01745329519943295)
Datum: D_NAD_1983_2011

When I add the raster, I get the Geographic Coordinate Systems Warning, saying that the data source of the raster ("GC_NAD_1983_2011") is different from the one used in the data frame. When I go in the Transformations options in the same dialogue box, and try to "Convert from GCS_NAD_1983_2011" to just "GCS_North_American_1983", there is no method available under the Using: dropdown, meaning I can't transform. What can I do to rectify this? How can I change the actual datum of the raster?


texts: Capturing coordinates via snap with QGIS

texts:
I would like to know if there is a way to get and extract the coordinates of a node (one of the node that conforms a polygon) in order to use these pair of coordinates to georeference a DXF file, using QGIS.

Is there a tool/procedure to do that. I mean, for instance, if there is a way to activate the snapping options for the identify tool, connect the cursor to a concrete node, and get its coordinates.

As far as I know, the Coordinate Capture plugin doesn't allow snapping options and I'm looking for something like that. Any hint or solution?

I'm using QGIS 2.10 on Ubuntu


texts: Computing unsupervised random forest classification in R?

texts:
I want to compute an unsupervised random forest classification out of a raster stack in R. The raster stack represents the same extent  in different spectral bands and as a result I want to obtain an unsupervised classification of the stack.
I am having problems with my code as my data is very huge.

Is it okay to just convert the stack into a dataframe in order to run the random forest algorithm like this?

stack_median <- stack(b1_mosaic_median, b2_mosaic_median, b3_mosaic_median, b4_mosaic_median, b5_mosaic_median, b7_mosaic_median)
stack_median_df <- as.data.frame(stack_median)

Here is the data as a csv file (https://www.dropbox.com/s/gkaryusnet46f0i/stack_median_df.csv?dl=0) - and you can read it in via:

stack_median_df<-read.csv(file="stack_median_df.csv")
stack_median_df<-stack_median_df[,-1]
stack_median_df_na <- na.omit(stack_median_df)

My next step would be the unsupervised classification:

median_rf <- randomForest(stack_median_df_na, importance=TRUE, proximity=FALSE, ntree=500, type=unsupervised, forest=NULL)

Due to my huge dataset a proximity measure can't be calculated (would need around 6000GB).

Do you know how to be able to have a look at the classification?

As predict(median_rf) and plot(median_rf) don't return anything.


texts: Converting polyline coordinates to lat/lon for use in google maps

texts:
I am given an ArcGIS polyline (or polypoint) from a WFS service (e.g. 1693669.5143, 5672605.5392), and am looking for a function (python, java, whatever) to convert it to a lat/lon so I can use it in Google Maps. I obviously have the projection as well.

I've googled this extensively but most of the results return ArcGIS's solution http://pro.arcgis.com/en/pro-app/tool-reference/data-management/add-xy-coordinates.htm - i.e. tool-specific

Is there an online tool or converter to do this, or (ideally) a function or library that I can use directly from code?


texts: How does assign_vertex_id() in pgRouting work?

texts:
I'm trying to follow this tutorial,

http://www.pgrouting.org/docs/foss4g2008/ch06.html

but on the step with assign_vertex_id() it shows an error for executing the SQL query:

SQL error:

ERROR:  column "’ways’" does not exist
LINE 1: SELECT assign_vertex_id(’ways’, 0.00001, ’the_geom’, ’gid’);
                            ^
In statement:
SELECT assign_vertex_id(’ways’, 0.00001, ’the_geom’, ’gid’);

I have not found any documentation or forums regarding how assign_vertex_id() should actually work.

The only thing that I know is that the first parameter should be a table_name and not a column name.


texts: Can I trigger a CartoDB infowindow programmatically?

texts: Is it possible to show a CartoDB infowindow for a particular point on a map programmatically? I want to open the info window when a marker is dropped on a Google Map.


texts:
I have some vector files that come without projection. So I have to define my own projection in QGIS. The information that comes with the data reads as following:

Projection: Lambert Azimuthal
Units: Meters
Datum: None
Parameters:
6370997.24063 (radius of the sphere of reference)
100 0 0.000 (longitude of center of projection)
45 0 0.000 (latitude of center of projection)
0.00000 (false easting (meters))
0.00000 (false northing (meters))

The Lambert Azimuthal is an equal area projection.
I tried to translate this into a projection definition as:

+proj=laea +lat_0=45. +lon_0=100. +x_0=0 +y_0=0 +datum=none  +units=m +a=6370997.24063 +b=6370997.24063

But my data is still wrongly projected.


texts: How to bring feature selection into print composer

texts: I have a PostGIS multipolygon layer in the map view and want to bring the feature selection into the print composer. Also, only the label for the selected feature should be visible in the map I want to create.


texts:
I used Pathfinder office 4.2 to export GeoXH data into RINEX format. It works right when I have only 1 point feature, however when I have multiple point features I get only 1 point feature in the resulting RINEX file.

I checked the RINEX file manually and it does only have 1 MARKER NAME and 1 MARKER NUMBER hence it only displays one of the features...
I also checked the points in the software and found that all the points seem to be named "point_generic" (even though they have different attributes).
Is there a way to change those names, or another way to export SSF to RINEX with multiple point features?


texts: Creating custom line symbols for Leaflet map?

texts:
I'm making an online map using leaflet and I'd like to add custom line symbols.

I am NOT trying to make custom markers. I say this because the internet is full of how to: tutorials for markers, and I've got that figured out, but i can't seem to find a tutorial for making line features.


texts: Why does ST_Area return different results for geog vs geom?

texts:
I have this simple query that returns very different values for a geom vs a geog:

SELECT ST_Area(the_geom::geography) as area1, st_area(the_geom) from neighborhood where city_id='ny-new-york-city'

+----------------+--------------------+
|area1           | area               |
+----------------+--------------------+
|1646766.08995329| 0.00017589982948304|
+----------------+--------------------+

The SRID of the_geom column is 4326. I don't get why geog column makes such a difference.

Which measurement is correct?


texts: Calculating custom HTML tables by feature for AGOL popups

texts:
I am working on a project that has all of the stops in a bus system. Each stop can serve more than one route and I have broken down ridership by route in each of the columns with a total route ridership, route boardings and route alightings. This data was an object that I transferred to a table. So there is a lot of fields that are 0 since that particular stop is not serving that route.

Each of these data fields have a standard field formula.

total route ridership = 'R_{route}'
total route boardings = 'R_{route}_BRD'
total route alightings = 'R_{route}_ALT'

In the table there is a list of routes separated by ', ' and I am using this to build a list in a for loop to create an html table. I have tried a couple of different methods to no avail to produce something that AGOL or Arcade finds acceptable.

Method 1: Arcade loop for building HTML

When Arcade outputs html it does not form a table

Method 2: Python loop for building Arcade

The first thing that I tried was building a Python script that would automate some arcade HTML and then used the field in a calculated field by just using return $feature.html

def htmlTable(routes):
list = routes.split(', ')
html = '

Route	Total Ridership	Avg Board	Avg Alight
'

for route in list:
    html += '' + route + '	"'
    html += '+' + f'$feature.R_{route}' + '+	'
    html += '+' + f'$feature.R_{route}_BRD' + '+	'
    html += '+' + f'$feature.R_{route}_ALT' + '+
'

html += '

'

return html

Method 1: Python loop for building HTML

When I put this together I am able to get the html formatting, but the popup ignores my attempt to call fields with Arcade.

def htmlTable(routes):
list = routes.split(', ')
html = '

Route	Total Ridership	Avg Board	Avg Alight
'

for route in list:
    html += '
' + route + '	'
    html += '{' + f'R_{route}' + '}	'
    html += '{' + f'R_{route}_BRD' + '}	'
    html += '{' + f'R_{route}_ALT' + '}
'

html += '
'

return html


texts:
Possible Duplicate:
  Re-projection in Geoserver

I've got a WMS service that publishes layers in WGS84 (EPSG:4326). I have several AutoDesk users that need to use the imagery in NAD83 (EPSG:2227). I have GeoServer 2.1 set up to cascade the WMS and reproject it.

Problem is that as soon as I turn on SRS handling to reproject to declared SRS the reprojected WMS layer fails. What am I doing wrong? Is GeoServer capable of doing the reprojection?


texts: Missing Geocoded Addresses in Attribute Table

texts:
I have geocoded about 13,000 addresses. I opened the Interactive Rematch Window and made adjustments where necessary. When I look at the Interactive Rematch window, I can see all of my addresses.

However, when I open the associated shapefile's attribute table, I only see 1,550 addresses.  I know that sometimes ArcMap will truncate large data files, but that's not happening here.

Also I don't have any edits to save in the Interactive Rematch Window.

Am I missing a step that refreshes the shapefile or what?


texts: In QGIS 3.4, how do you perform a does not contain filter? I know how to perform a contains filter, with the LIKE / ILIKE operators, but how do you perform a does not contain filter?


texts: Converting a pixel coordinate in GeoTIFF to a latitude and longitude?

texts:
I have a GeoTIFF with the following meta information:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]

How to I find the geographic coordinate (lat, lon) in Web Mercator of a pixel in the image (e.g. x=500, y=600)


texts: Attribute table values randomly disappear or are changed

texts:
Im using 3.10.14 version of QGIS and something very strange started happening. I performed several attraction distance and attraction reach analyses (using the Place Syntax Tool) that I saved in the same layer. Then I would usually duplicate that layer (right click -> duplicate layer) to change the visualization settings to create different maps. I have been doing this over the course of 2 months. When I opened the project recently I saw that the values for some of the analyses in the attribute table had changed either to NULL or to a mix of 1, 0 and NULL. After I redo the analysis it looks fine again, but when opening it after a couple of hours it has changed again.
Has this happened to anyone else and do you have any suggestions as to what Im doing wrong?
Update: The file was saved as geopackage. I'm attaching screenshots of most part of the attribute table and the field types.


texts: Vector layer not displaying in OpenLayers

texts:
I have this code from an Udemy course. Everything works except I can not visualize the vector layer. I can see the JSON file - wfs_url correctly but not visualize vectorSource (the data points selected in wfs_url).

<!DOCTYPE html>

    Working with Vector Data

  var view = new ol.View({
      center: [-9015705,4026797],
      zoom: 11
      });

   var map = new ol.Map({
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      })
    ],
    target: 'map',
    view: view
  });

var theurl = 'http://localhost:8050/geoserver/richland_postgres/wms';

var bk79 = new ol.layer.Image({
      source: new ol.source.ImageWMS({
        url: theurl,
        params: {'LAYERS': 'richland:bk79'},
        serverType: 'geoserver'
      })
    });
 bk79.setOpacity(0.3);

var cqlfilter = 'DWITHIN(geometry,Point(503909 3768402),5000,meters)';

var schoolcql = new ol.layer.Image({
       source: new ol.source.ImageWMS({
        url: theurl,
        params: {'LAYERS': 'richland_postgres:schools','CQL_FILTER': cqlfilter},
        serverType: 'geoserver'
      })
    });

wfs_url = 'http://localhost:8050/geoserver/richland_postgres/wfs?service=WFS&' +
          'version=1.1.0' +
          '&request=GetFeature' +
          '&typename=richland_postgres:schools&' +
          '&CQL_FILTER=' + cqlfilter + '&' +
          'outputFormat=application/json&srsname=EPSG:3857&' +
          ',EPSG:3857';
prompt('',wfs_url);

var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: function(extent) {
      return 'http://localhost:8050/geoserver/richland_postgres/wfs?service=WFS&' +
          'version=1.1.0' +
          '&request=GetFeature' +
          '&typename=richland_postgres:schools&' +
          '&CQL_FILTER=' + cqlfilter + '&' +
          'outputFormat=application/json&srsname=EPSG:3857&' +
          ',EPSG:3857';
    },
    strategy: ol.loadingstrategy.bbox
  });

  var vector = new ol.layer.Vector({
    source: vectorSource
  });

map.addLayer(bk79);
map.addLayer(schoolcql);
map.addLayer(vector);

</script>


texts: Why measure line in QGIS result in millimeter?

texts: I am using QGIS 3.4.4 I just want to measure the distance between a tree with its shadow  from a raster data, but when i use measure line the result is in millimeter. I just upgrade the software to the newest version, when I used the previous one, the measure line works, can someone help me please...


texts: How to determine Sentinel Productnames for EarthEngine/AWS/GoogleStorage

texts:
How would I know the granule name for sentinel data in Google Earth Engine? I tried to look up the granule name in the "index.csv.gz" (gs://gcp-public-data-sentinel-2/index.csv.gz) and "scihub" but cannot get the code running. Here are the field names in the index.csv.gz GRANULE_ID,PRODUCT_ID,DATATAKE_IDENTIFIER,MGRS_TILE,SENSING_TIME,TOTAL_SIZE,CLOUD_COVER,GEOMETRIC_QUALITY_FLAG,GENERATION_TIME, NORTH_LAT,SOUTH_LAT,WEST_LON,EAST_LON,BASE_URL.

On the other hand this is the sample granule ID I got from scihub: S2A_OPER_MSI_L1C_TL_EPA__20170512T230414_A002288_T50PQA_N02.04.

I am struggling how to determine sentinel ID (specifically Granule ID) in Google Earth Engine.


texts: How do get North, East, South, or West bounds from Cesium map?

texts: I see that there is an option that allows for one to set the view if you provide the coordinates for these bounds, but nowhere in the documentation do I see anything pertaining to retrieving these values as defined by the canvas container / rectangle.


texts: Tippecanoe tile-join tilejson v3 and v3.1

texts:
I have merged tiles with tile-join successfully after merging metadata of tilejson changed to v2.
is there any way to set metadata to v3 or v3.1 to tippecanoe.
or another tool which support v3 or v3.1 vector tiles schema?

github issue.
https://github.com/mapbox/tippecanoe/issues/852


texts: How to split a polygon in 9 equal parts but I want to use different orientation on QGis?

texts:
I want to split a polygon of a terrain that I generated through coordinates on QGis. But at the same time, I want to do perform such an activity by testing different orientations for the terrains being divided.
I know that there is a lot of calculus involved in irregular things like that, but I wasn't expecting it to be so hard to manipulate different orientations for the polygon divisions. In this case, I got (with some trouble), the following orientation: "|||||". And it is not even good.

However, I want the terrain to follow this orientation: "////" or maybe "\\". I want the division to happen in diagonal sections instead of that. And on top of that, I barely can make it look decent because since it is an irregular polygon and I want 9 equal parts out of it, the Plugin "Divide Polygon" creates some crazy divisions in a block like aspect.
About the 2 different sections with 10 divisions on it, these are actually the method for me in order to create "stripes". In this case, sometimes the Divide Polygon plugin creates a perfect stripe. If not, I need then to divide the area into 20 parts and work around my way in order to assemble the division in equal stripes. As you can imagine, this is a tedious task and sometimes the "Divide Polygon" plugin gives you some crazy divisions (to my surprise, the plugin also is prone to make gross mistakes if you use it intensively. Like dividing an area set to 4 parts into 3 parts).
One last thing: Although I want 9 divisions on this terrain, I only got 8 showing in this image. This is because sometimes you can get a perfect vertical cut when dividing by two a large area. Then you can try to divide the new area in two as well and such. I am doing this but I do not like to deal with QGIS in this way.
I tried also several other answers made around that appear to correlate cases but I believe that this particular case that I am asking here wasn't enough discussed on other threads. I need a simple and direct answer, and I tested several other methods and it wasn't ever mentioned my specific problem on such cases.


texts: How to create proportional circles in QGIS 3.4

texts: I am working with the polygon features and I would like the values to be represented in terms of proportional symbols (points). Values range from 0 to 30. I would like that: zero value = no circle; values equal to 1 = small circle; values equal to 2 = ..; values between 2-10 = ..; Over 10 = biggest circle. It is OK to use the same colours for the different dots. How do I do it in QGIS 3.4?


texts: Correctly scale topographic SRTM data in OSM

texts:
I am trying to create contour maps in OSM with scales of 1:50,000 and 1:250,000 using 30-meter NASA SRTM data from here. Examples of what I'm trying to achieve can be seen on LINZ NZ Contours 1:50k and NZ Contours 1:250k pages:

Using methods strongly based on the OSM Contours page, I've managed to manipulate the SRTM data to render as follows:

I then tried to get something similar to Contours 1:250k by playing around with the rendering queries, but I've noticed my contour shapes are not exactly the same, and I fear I may be going about this the wrong way. Essentially, my query is almost identical to what's in the OSM Contours query section, except that I've modified the id: contours50 query to match WHERE height::int % 110 = 0. My attempt is on the left, the LINZ 1:250k version is on the right:

The commands I'm using to get all this into the database:

# Download .hgt file from NASA Earthdata

# Store .hgt file name in variable
source_file="S37E174"

# Prepare file for import
gdal_translate -q -co TILED=YES -co COMPRESS=DEFLATE -co ZLEVEL=9 "$source_file.hgt" "$source_file.tif"
gdalbuildvrt -q -overwrite SRTM1.vrt "$source_file.tif"
gdal_translate -q -co TILED=YES -co COMPRESS=DEFLATE -co ZLEVEL=9 -co PREDICTOR=2 -projwin 174.692262 -36.780183 174.818605 -37.028786 SRTM1.vrt srtm_30m.tif
gdal_contour -i 10 -a height srtm_30m.tif srtm_30m_contours_10m

# Import into database
cd srtm_30m_contours_10m/
shp2pgsql -p -I -g way -s 4326:900913 contour.shp contour | psql -h "$my_ip" -U postgres -d gis
shp2pgsql -a -g way -s 4326:900913 contour.shp contour | psql -h "$my_ip" -U postgres -d gis

My question: is this the correct way of scaling topographic information? Is there a better way of doing this?


texts: Joining Excel table to shapefile in ArcMap?

texts:
I have a shapefile of lower layer super output areas (LSOA) for England and Wales and I have a spreadsheet of the indices of deprivation for each of the LSOA's in England. I want to join the deprivation data to the shapefile so it can be mapped in ArcMap. I tried joining the table through the Join and Relates > Join option and put in the fields for the LSOA data for both shapefile and spreadsheet (the data that matches). the attribute table added the datafields to the attribute table from the excel spreadsheet but all of the data is blank!!

I have tried looking at the field "type" as ive read that it has been an issue before but it seems to be OK.

Is there anything anyone can suggest?


texts:
The PCS and the GCS for the shapefiles, the DEM below the polygon file and the data frame are all the same. But the coordinates are completely off. When I bring newly downloaded data onto the map it lines up with the coordinates you see in the lower corner of the screen. I can't get everything else to line up.

I've tried reprojecting everything, changing Geoprocessing Environments, defining projection etc.

Could the  data be corrupted?

Should I just redownload everything from the start?


texts: Displaying multiple regions shapefiles on iOS device?

texts:
How can I display multiple shapefiles on an iOS device?

My client is asking to display all the states in a different color, and again cities in a different color and all the fields in a different color with lat/longs.

I was not able to proceed further with such a huge requirement so how should I start?

What are the requirements for moving forward?

Are there any GiS specialists needed for me to know completely about the shapefiles?


texts:
Is there any API out there which can give me Polygon value for a zipcode? I really thought that it would be fairly easy for me given the number of Maps APIs out there. Either Google Maps API or something else would be able to give me that data. But a plain Google search was not able to come back with any results which were satisfying.

The I realized that the Zipcodes can not be "strictly" defined as polygons. They belong to another type of geometrical element called "Boundaries".

  All I want it a container inside a container inside a container

However, the smallest molecule in the equation would be a zipcode. The other containers can be Polygons OR Boundaries. I am not sure if that would make a HUGE difference. However, what I want essentially is a Geometrical representation of a zipcode boundary so that my SQL Server Db can understand it and I can query it from there

Any good suggestions on how can I enter a Zipcode and get a respective Polygon/boundary value for it?


texts: QGIS record containing polygon for each lat/lon point

texts:
I have two layers in QGIS. One is a shapefile of states in the US. The other is a csv with lat/lon coordinate points, with data at each point.

I want to export a csv from QGIS that has a column that records what state each data point falls in.

For example, let's say I have 4 columns in my csv file. First is location ID, second is Latitude, third is Longitude, and fourth is count (just a number). One of the rows might look like this:

123456, 41.2, -97.3, 4

Since this point is in Nebraska (in real-life and according to my shapefile), I want to have a new column for state, so that the above row would now look like this:

123456, 41.2, -97.3, 4, Nebraska

Is this something I would use 'Join Attributes by Location' for? I'm trying that now, but its suspiciously slow.


texts: How to get spatialite_osm_net from spatialite

texts:
I am using windows 7. I downloaded spatialite-tools (precompiled) from the following site. But there is no "spatialite_osm_net" tool available. For one of my project I was trying to convert the OSM road network to spatialite road network.

http://www.gaia-gis.it/spatialite-2.3.1/binaries.html

It looks like without this tool It is not possible to do it.

For your information, I tried to make this tools from the following site according to instruction using MinGW. But It does not work when I come to $make command. I do not know is there anything missing? or how Should I do?

$ cd c:/spatialite-tool-2.4.0

$ export "CFLAGS=-I/usr/local/include"

$ export "LDFLAGS=-L/usr/local/lib"

$ export "PKG_CONFIG_PATH=/usr/local/lib/pkgconfig"

$ ./configure --target=mingw32

$ make

$ make install-strip

So now how can I get spatialite_osm_net???


texts: Making lists, looping a function through a list and controlling output names

texts:
How do you output a list of multiple featurelayers from an mxd while controlling the output names of the shapefiles and overwriting the contents of the shapefile outputs if the the file already exists? Currently I'm using both the FeatureClasstoFeatureClass_conversion() and CopyFeatures_management() functions. I can get one featurelayer to output and overwrite a shapefile in a different location. But how to do this for a list of shapefiles?

import arcpy
Mxd = arcpy.mapping.MapDocument("CURRENT")
arcpy.mapping.ListLayers(Mxd)
arcpy.env.overwriteOutput = True

sampleinputList = ["input_a", "input_b", "input_c"] #list of feature layers
sampleout_path = "T:/path" #output path
samplenames = ["input_a_name", "input_b_name", "input_c_name"]
sampleOutputnamesList = ["T:/path/input_a_name", "T:/path/input_b_name", "T:/path/input_c_name"] # each input corresponds to a respective output name
for featureclasses in sampleinputList:
    arcpy.CopyFeatures_management(sampleinputList,sampleOutputnamesList)
    #neither of these list structures worked :/
for featureclasses in sampleinputList:
    arcpy.FeatureClasstoFeatureClass_conversion(sampleinputList,sampleout_path,samplenames)

both instances give me a
    Runtime error
    Traceback (most recent call last):
    File "<string>", line 2, in <module>
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py",    line 2281, in CopyFeatures
    raise e
    RuntimeError: Object: Error in executing tool

Here's the part of the code I think should be amended:

#Filename mapping Dictionary
fileNames = {
'unique layer namea' : 'unique_shapefile_renamedA',
'unique layer nameb' : 'unique_shapefile_renamedB',
'unique layer nameC' : 'unique_shapefile_renamedC',
}

outPath = "T:\\path" #path to folder where new shapefiles will be output

for featureclass in listLayers: #this loops through the arcpy function seen on last line
    if featureclass.name.lower() in fileNames:
        outputName = '{}.shp'.format(fileNames[featureclass.name.lower()])
        outputPath = os.path.join(outPath, outputName)
        arcpy.CopyFeatures_management(featureclass, outputPath)

if the original syntax of the Copy Features management function is
    CopyFeatures_management (in_features, out_feature_class,)
then shouldn't featureclass in last line of the block be defined in the way outputPath is defined?


texts: Generating HGT Files from SRTM 90m Digital Elevation Database v4.1?

texts:
This may be a really dumb question and I apologize if it is - this is pretty new to me.

I'm using the HGT files from https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/ with a program currently and want to upgrade to more recent, better quality, data.

I downloaded the SRTM 90m Digital Elevation Database v4.1 ZIP file and see that HGT files aren't one of the formats supplied.

How can I convert what is supplied into HGT files with the same naming convention as what I currently have?

Is that possible or do I need to modify my application to suit this other format?


texts: Single WMS from single mapfile with several rasters in MapServer

texts:
I have several rasters in a single mapfile, I put them as layers like here:

LAYER
    NAME            "Layer 1"
    GROUP           "All Layers"
    TYPE            RASTER
    STATUS          ON
    DATA            "layer1.png"

    PROJECTION
        "init=epsg:4326"
    END

    METADATA
        "wms_title"         "Layer 1"
        "wms_srs"           "EPSG:4326"
        "wms_group_title"   "All layers"
        "wms_abstract"          "No abstract"
        "wms_server_version"        "1.1.1"
        "wms_format"            "image/png"
        "wms_include_items"     "all"
        "gml_include_items"     "all"
        "gml_geom_type"         "raster"
    END
END

LAYER
    NAME            "Layer 2"
    GROUP           "All Layers"
    TYPE            RASTER
    STATUS          ON
    DATA            "layer2.png"

    PROJECTION
        "init=epsg:4326"
    END

    METADATA
        "wms_title"         "Layer 2"
        "wms_srs"           "EPSG:4326"
        "wms_group_title"   "All layers"
        "wms_abstract"          "No abstract"
        "wms_server_version"        "1.1.1"
        "wms_format"            "image/png"
        "wms_include_items"     "all"
        "gml_include_items"     "all"
        "gml_geom_type"         "raster"
    END
END

LAYER
    NAME            "Layer 3"
    GROUP           "All Layers"
    TYPE            RASTER
    STATUS          ON
    DATA            "layer3.png"

    PROJECTION
        "init=epsg:4326"
    END

    METADATA
        "wms_title"         "Layer 3"
        "wms_srs"           "EPSG:4326"
        "wms_group_title"   "All layers"
        "wms_abstract"          "No abstract"
        "wms_server_version"        "1.1.1"
        "wms_format"            "image/png"
        "wms_include_items"     "all"
        "gml_include_items"     "all"
        "gml_geom_type"         "raster"
    END
END

The problem is that I want to have a single WMS layer (named as in GroupTitle) as output from this mapfile (because it is only one mapfile), instead I have four WMS layers (three coming from each LAYER and one from the GROUP LAYER).

How can I achieve what I want? I use Mapserver 6.4.1 in OL 2.11.


texts:
I have a code that is supposed to estimate the elevation (which is in percentile per pixel) where a certain value is minimal. I do this gridwise, but instead of using a grid, i want to use a reduced resolution image, with the value per pixel. I am able to find the minimum value (MOD_rp2) but am unable to link that for each pixel to the elevation band it originates from. How would I do this? Even better would be to link it to the actual elevation of that band, so the elevation for the corresponding percentile (from hgt_p)

//import Datasets
var MOD = ee.ImageCollection('MODIS/006/MOD10A1')
   .select('NDSI_Snow_Cover').filterDate('2000-10-01','2010-09-30')
   .map(function(img){return img.gte(40)});

//precentile list
var start = 5;
var end = 95;
var step = 5;
var percent = ee.List.sequence(start,end,step);

//elevation bands based on pdf
var hgt_p = hgt.reduceResolution({reducer:ee.Reducer.percentile(percent),
    maxPixels:15000,bestEffort:false}).reproject({crs:crs,scale:scale});

//create snow/no snow map (1 for land above percentile, and snow below)
var hgt_c = hgt.gt(hgt_p);
var x = function(img){return img.subtract(hgt_c).abs()};
var MOD2 = MOD.map(x);

// sum up all pixels which need to be minimized
var MOD_rp = MOD2.map(function(img){var rr = img.reduceResolution({
    reducer:ee.Reducer.sum().unweighted(),maxPixels:15000,
    bestEffort:false}).reproject({crs:crs,scale:scale});
    return ee.Image(rr)});

//minimize, to add: extract elevation from hgt_p which confirms the correct pixel)
var MOD_rp2 = MOD_rp.map(function(img){var a = img.reduce(ee.Reducer.min()).clip(C_A);
    return a;});

the full code can be found here


texts: pgRouting cannot find route at intersection

texts:
I have a simple path that is defined by two polylines:

The first line goes from point 1 to point 2, and the second from 3 to 4. After running the pgr_createTopology script the source and target's values are consistent with that. The pgr_analyzeGraph returns the following results:

ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 2
NOTICE:                            Dead ends: 4
NOTICE:  Potential gaps found near dead ends: 1
NOTICE:               Intersections detected: 1
NOTICE:                      Ring geometries: 0

My problem is when I try to run the pgr_dijkstra script. That only returns me a edge between the points (1,2) and the points (3,4). I can't get a result from 1 to 4 for example.

How can I fix this problem?


texts: Get Z value from point in LAS Dataset using ArcObjects

texts:
Using ArcObjects (C#) I am trying to build a tool that will extract the Z value from a point in a LAS Dataset when the user clicks on it (then use that value elsewhere). My intention was to use the IIdentify2 interface to extract the Z value, but that doesn't seem to support LASDatasetLayer's

I have looked at a couple of other interfaces that might be of some use like ILasPointInfo or even an event. I feel like there should be one for clicking on a point in a LAS Dataset.

Unfortunately the documentation on this whole area of LASDatasets is pretty limited. Also I am pretty new to C# and ArcObjects so its pretty hard to work out how this all fits together.

I've rewritten the code and now have a path to the desired outcome. However there is a new problem. I have tried to use the 'DisplayTransformation.ToMapPoint()' method on the arg.X and arg.Y values from the MouseDownEvent (see code below). This simply inverts the Y value of the 'device coordinates'. Because of this by the time that I try to select a LAS point, the IGeometry object that is required in the GetLasPointInfo method is still in 'device coordinates' and is being used to search for a LAS point that is in 'real-world coordinates' - therefore the error below is returned. Are there other interfaces I should be using or is there a completely different way to do this?

As you can see I've tried to set the spatial reference to that of the Scene (GDA_1994_MGA_Zone_55). I have also tried to use the Project method on the point. Most likely in vein.

If you think of anything I've left out, let me know.

Environment:
ArcScene 10.3.1, ArcObjects .NET

New code:

protected override void OnMouseDown(ESRI.ArcGIS.Desktop.AddIns.Tool.MouseEventArgs arg)
    {
        try
        {
            ITrackCancel tCancel = new CancelTrackerClass();

            ISxDocument sxdoc = ArcScene.Application.Document as ISxDocument;
            IScene scene = sxdoc.Scene as IScene;

            ILayer layer = scene.Layer[0] as ILayer;
            MessageBox.Show("Layer: " + layer.Name);

            IPoint clickPoint = new Point();
            clickPoint.PutCoords((double)arg.X, (double)arg.Y);
            MessageBox.Show(arg.X + ", " + arg.Y);

            ISpatialReference sRef = scene.SpatialReference;

            ISceneGraph sceneGraph = scene.SceneGraph;
            IScreenDisplay screenDisplay = sceneGraph as IScreenDisplay;
            IPoint scenePoint = screenDisplay.DisplayTransformation.ToMapPoint((int)clickPoint.X, (int)clickPoint.Y) as IPoint;
            MessageBox.Show("New coord sys for scene point: " + scenePoint.SpatialReference.Name);
            IPoint projScenePoint = new Point();
            scenePoint.Project(sRef);

            //IEnvelope projEnjGeoClick = geomClickPoint.Envelope as IEnvelope;
            MessageBox.Show("Projected coords: "+ scenePoint.X + ", " + scenePoint.Y);

            //ILayer layer = scene.Layer[0] as ILayer;
            MessageBox.Show("Layer: " + layer.Name);
            ILasDatasetLayer lasLayer = layer as ILasDatasetLayer;

            ILasDataset lasDataset = lasLayer as ILasDataset;
            ILasPointCloud lasPC = lasDataset as ILasPointCloud;
            IArray lasArray = lasPC.GetLasPointInfo(tCancel, scenePoint, 0) as IArray;

            ILasPointInfo lasPtInfo = lasArray.Element[0] as ILasPointInfo;
            WKSPointZ wksPtZ = lasPtInfo.get_WksLocation();
            double zValue = wksPtZ.Z;
            MessageBox.Show("Z Value: " + zValue);
        }

        catch (Exception ex)
        {
            MessageBox.Show(ex.Source + "\n\n" + ex.Message + "\n\n" + ex.InnerException + "\n\n" + ex.StackTrace);
        }

    }