texts: How can I use GDAL Merge with multiple files in one folder in python 2.7?

texts:
I have multiple files in one folder which are called, for example: 049.h11v11.hdf, 049.h11v12.hadf, 049.h11v13.hdf, 065.h11v11.hdf, 065.h12v12.hdf, and so on. The first 3 numbers are the julian day of MODIS LAI (MOD15A2) product.

The objective is to do daily mosaics with about 3 years of modis MOD13A2 products. At this point I have a folder named "LAI", that contains the years folder, e.g 2000, and so on. Those folder contains the .hdf files I named in the last paragraph. So, the idea is that the python code can run for all folders and get all same days and do the mosaic.

So far I've though about something like this:

#Libraries
import gdal
from gdalconst import *
import osgeo
import osr
import osgeo.gdal as gdal
import osgeo.osr as osr
import numpy as np
import glob
import os

#Roots:
RutaLai="/mnt/discoenred/jorge/scripts/MOD15_A2_mosaic/LAI/"

#Create dir
if not os.path.exists(RutaLai): os.mkdir(RutaLai)
.
.
#in between this are the reprojection and the multiplication by the LAI factor for each .hdf file.
.
.
#Here is were the mosaic code starts:
for raiz, carpetas, archivos in os.walk(RutaLai, topdown=False):

    arreglo = carpetas

diaTemporal = None
diaVariable = []
diaArreglo = []

for carpeta in arreglo:
    os.chdir(RutaLai+carpeta) #carpeta is defined as as the year taken from the MOD13A2 name, i.e carpeta = x[9:13]+"/" ; were x goes inside a for x in lista2, that just take the year inside the name of the file. (it is inside the code for the LAI factor, I don't put it because its really huge. If it is really needed then I can put it)
    #Lai images lists
    lista3=glob.glob("*.hdf")
    lista3.sort()

    for archivo in lista3:
        dia = archivo.split(".")

        #Definition of temporal day (diaTemporal) variable to avoid error.
        if (diaTemporal == None):
            diaTemporal = dia[0]

        if (diaTemporal == dia[0]):
            diaVariable.append( dia[0] )

        else:
            if (diaVariable != None):
                diaArreglo.append( diaVariable )

        diaTemporal = dia[0]

    #RutaDestinoMosaico = RutaMosaic + carpeta + "mosaico_" + ".hdf"
    #os.system("gdal_merge.py -o "+ RutaDestinoMosaico + " -of ENVI -n 0 '"+archivo1+"' '"+archivo2+"'")

The idea is that the code go to all e.g 049 files in one folder and do the mosaic that will be saved in a folder inside the year called "mosaic".. the files should go like "mosaic_049.hdf" (for the example given above). In order to understand it better I took a screen shot of the directory:

It is important to notice that the picture (inside year "2000") contains the files that are inputs for the mosaic. For this all I'm using a Linux machine with Ubuntu.


texts: Replicate Geodatabase and append to main SDE database

texts:
I have my ArcSDE GDB on SQL Server 2005. This we are maintaining at our end at head office.
Now we have installed ArcGIS editor license to regional offices. GIS technicians at regional offices can't edit data in main database at SQL server. thus, we have decided to provide them offline replica of the Arc SDE GDB.

Once they finish their edits, weekly/ monthly basis we will recieve GDB from each regional office.

Question is that, How we will upload their GDB into our main ArcSDE GDB without duplicating data in the server? or how we can upload only the features which are updated by the GIS technicians at regional offices?


texts: Automate - convert Labels and attach to predefined annotation class

texts:
I am developing pythonaddin tool which as a result should update/accumulate a predefined geodatabase annotation standard class. I am stucked at the point where I believe that arcpy has limitations within annotation operations.
With click on the button, I am aiming to show labels with given expresion, to set up data frame scale, convert label to annotation and append it to already predefined annotation.
In the tool, I am able to finalize steps with labels, but not to convert and update existing annotation layer using arcpy.
 if layer.supports("LABELCLASSES"):
    for lblClass in layer.labelClasses:
        lblClass.showClassLabels = True
 lblClass.expression ='"%s" & [export.DARST1] & "%s"' % ("<CLR red='255' green='0' blue='0'>", "</CLR>")
 layer.visible = True
 layer.showLabels = True
 self.df.scale = 10000

The arcpy methods arcpy.TiledLabelsToAnnotation_cartography and then arcpy.AppendAnnotation_management do not really help, because they create new annotation features and then trying to append it to predefined annotation is getting hard - makes it difficult and slow, but I will not say imposible.
Moreover, the last steps I need, would be possible to execute manually (pic below), therefore my assumption is that it can to be performed by adding some lines of code from arcobjects?!

Any help if it is possible with python for arcobjects is welcome, otherwise I assume the problem can be solved by working with ArcObjects for .NET/C#...


texts: Define original image ID for each pixel of Quality Mosaic in Google Earth Engine

texts:
I used a set of LS8 TOA images to produce a quality mosaic based on NDVI values. I wonder if there is a way do figure out a source / origin image ID for each pixel of the result mosaic. I guess it should be some kind of a thematic raster of the same resolution.

Source code I used to produce mosaic is attached below.

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');
// Get the least cloudy image in 2015.
var image = ee.Image(
  l8.filterBounds(ee.Geometry.Point(-122.2929, 37.8206))
    .filterDate('2015-01-01', '2015-12-31')
    .sort('CLOUD_COVER')
    .first()
);
Map.setCenter(-122.2929, 37.8206, 10)
// Compute the Normalized Difference Vegetation Index (NDVI).
var nir = image.select('B5');
var red = image.select('B4');
var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};
var withNDVI = l8.map(addNDVI);
// Make a "greenest" pixel composite.
var greenest = withNDVI.qualityMosaic('NDVI');
// Display the result.
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.3};
Map.addLayer(greenest, visParams, 'Greenest pixel composite');


texts:
I have a CSV file with X/Y coordinates that I would like to reproject to have long/lat coordinates.
Unfortunatelly, reprojecting or setting CRS does not work.

I have a layer of polygons which is correctly projected (EPSG: 4326), showing coordinates in degrees:

I have another layer of points, which appears in the correct CRS (EPSG: 4326), but the coordinates appear in X/Y.

This points should appear on top of the previous polygon layer:

I have tried to:

Change the CRS on the layer properties
Set the CRS of the layer
Toolbox -> Define layer projection
Toolbox -> Reproject layer

I have tried with the on memory layer, saving it as CSV, saving it as SHP and reimporting...

What am I missing?

QGIS version 3.8.0-Zanzibar


texts: Geoserver Reprojection Results in Northerly Shift Compared to GDAL

texts:
I have data in a postgresql data store in EPSG:2163 and am displaying that using OpenLayers in EPSG:3857. This results in a ~25 km northerly shift when displayed either in QGIS or OpenLayers.
To test the reprojection I used GDAL and the Geoserver reprojection console, which confirmed that the x-coordinates were reprojecting correctly but the y-coordinates were artificially high.
GDAL:
gdaltransform -s_srs EPSG:2163 -t_srs EPSG:3857
1494718.24465831 -1025289.47657808
output = -9310869.18028598 4090931.77295946 0 (this displays correctly).
Geoserver:
Reprojection Console
output =-9310869.180285975 4115158.942481106 (this displays north of expected).
Resulting in a difference of 24.2 km.
It seems to me that this a projection problem within Geoserver, any ideas what the problem is or how I can get around this? If I reproject my source data in postgresql to EPSG:3857, there are no projection problems but I want to keep my source data in EPSG:2163.


texts: How to set a definition query for an AGS layer?

texts:
Is it possible to filter out some features using a definition query set from the client on a AGS layer?

I've checked out the map server layer but it only supports find and identify. Also looked at WMS, WFS and the Feature Service but I can't find a way that works. Do I need to extend the AGS to get this functionality?

Update:
What I'm trying to do is opening an MXD in an ArcGIS Engine application. The MXD has an ArcGIS Server map server layer. I would like to filter out some of the features in one of the ArcGIS Server layers by using a definition query.


texts:
I am new in mapnik, and I found mapnik at OSM wiki page.

And I am interested in it, but I have some question.

1 the mapnik work scope

It seems that it is a library for map rendering(creation), however I wonder if the mapnik is bind to OSM data only?

I mean can mapnik work on data procesed by ArcGIS or something else?

In our production environment, we use ESRI arcgis to create,edit,render our map.

2 use mapnik to generate raster tiles based on Arcgis data source.

And we use arcgis server to generate the raster tiles, but it often cost a long time, can it only work on windows.

So I wonder if we can use mapnik to generate our raster tiles according to our ESRI based data(.gdb database)?

If it is possible, I wonder how about the style rules of the layer? Do you edit the style rules manually? Is it possible to translate the rules from Arcgis .mxd file?

3 vector tiles

BTW,we are doing some research work on the vector tiles( we tried to use it to replace the raster tile some day). And I wonder if mapnik support to create vector tiles based on our Arcgis .gdb data source?

And in the wiki page I found this:

TileStache
  Providers
      VecTiles - MVT (Mapnik Vector Tiles) and GeoJSON from PostGIS
      Vector - GeoJSON (and others) from OGR drivers
      PostGeoJSON - GeoJSON from PostGIS
  Hand-Rolled Vector Tiles – TileStache (blog)

What does the MVT (Mapnik Vector Tiles) mean?


texts:
A friend of mine gave me DEM resolutions of cells 5x5m, but there is something wrong with interpolation. When I try to make a cross-section of the terrain I get a pretty inaccurate cross-section, with weird peaks and valleys.
What should I do to get a smoother DEM?
These are the extensions:

cross-section (There is no waterfall there :-)


texts: Writing python plugin: How terminate plugin from script without closing QGIS?

texts:
I'm writing a QGIS conversion plugin in python to be able to create Obstacle-files for a model from shapefiles.

What I want to do now is just to close the plugin cleanly. This means, to close the GUI, to delete the variables and so on. I tried to use sys.exit(), but it also closes QGIS.

How do I only close the plugin?


texts: Comparing two maps with nice web interface?

texts:
Is it possible in Leaflet or OpenLayers to compare two maps with a bar to slide between the two?

Examples include:

http://www.esri.com/services/disaster-response/hurricanes/hurricane-sandy-the-aftermap
http://www.arcgis.com/apps/StorytellingSwipe/index.html?appid=c62112f6ba834227978d8f1cf2b6f5af
http://www.abc.net.au/news/specials/qld-floods/ (not a map)

I know it's possible, I don't know how to do it though. Anyone got any ideas? Examples?


texts: I'm wondering how to geocode addresses in R Studio ? Is there a link to Google Developers API ?


texts: How to connect a argis, sde version through stand alone python script. SDE.Mapping is child version of SDE.Default and I want to make child version of SDE.Mapping . I am using direct sde services


texts: WFS layer doesn't display in OpenLayers2

texts:
I have a few layers in a PostGIS database that I serve through GeoServer and that I try to display in OpenLayers.

When I had the GeoServer installed on localhost, it worked flawlessly. Now that I migrated the GeoServer on a remote server, the WFS is no longer displayed on my OpenLayers map.

Looking at Firebug, the request is sent to the server but no answer is received.

What is even weirder is that I can access the WFS from QGIS without any problem. Also, even the Geoserver examples (WMS though) work perfectly.

The page is here: http://www.henriod.info/khub2/khub2.html

Do you see any problem with this part of the code?

 var schools = new OpenLayers.Layer.Vector("Schools", {
               styleMap: style_schools,
               strategies: [new OpenLayers.Strategy.BBOX()],
               projection: new OpenLayers.Projection("EPSG:4326"),
               protocol: new OpenLayers.Protocol.WFS({
               version: "1.1.0",
                               url:  "http://thething.disco.unimib.it:12080/geoserver/wfs",
               featurePrefix: "humadat",
               featureNS: "http://thething.disco.unimib.it:12080/geoserver/khub2",
                               featureType: "schools",
               geometryName: "the_geom"
                           })
            });

Or could the problem be somewhere else?


texts: Getting the correlated FeatureDataRecord (or metadata) for a given geometry?

texts:
I'm using SharpMap 2.0, and I'm successfully loading a shapefile and parsing through it. The particular shapefile that I'm working with has about 700 polygons in it, and each of those polygons has meta data associated with it (the metadata is in the .DBF file).

How can I get at the metadata for each polygon? I'm not seeing any methods that allow access to this data...

I've searched Google, these forums and have pretty much exhausted my resources trying to figure this one out.

Any help would be greatly appreciated!


texts: compile QGIS-plugin for Windows: problem with gdal_i.lib version

texts: While compiling my plugin for QGIS 2.2 there are no errors, but it ist still not loading in QGIS because of an dll error. (function not found) I think the problem is at the gdal version. The QGIS-About says it is compiled with GDAL 1.10.1 but my plugin is compiled with GDAL 1.11. So I tried to compile my plugin with referencing the gdal_i.lib form the installed QGIS named "C:\Program Files (x86)\QGIS Valmiera\lib\gdal_i.lib" but cmake says while configuring that it is gdal version 1.11. Why? I tried it with the lib fom QGIS Lisboa - same result. Although the libs have different sizes and dates. Where is the problem?


texts: Using ArcGIS python window to change definition query on layer?

texts: I post a very simple script to the python window to change the definition query on a layer.  I save the map and load a new template then try to run the query again and nothing happens.  I have recently upgraded to 10.4 and when running 10.3 I never had this kind of issue.  I have tried deleting all variables at the end of the script.  Nothing works.  It will only run one time and I have to close out the ArcMap session before it will work again.


texts: Activate layers that are in a checkbox list

texts:
I'm trying to do a code to activate layers that are in a checkbox list, the code that I'm trying to use is this:

    MxDoc = My.ArcMap.Document
    Maps = MxDoc.Maps
    If chxActivated.SelectedItem = True Then
        GxLayer.Layer.Visible = True [COLOR="#FF0000"]HERE IS THE PROBLEM[/COLOR]
        MxDoc.UpdateContents()
        MxDoc.ActiveView.Refresh()
    End If

I add the layers with this code:

    MxDoc = My.ArcMap.Document
    Map = MxDoc.FocusMap
    'new layer
    GxLayer = New GxLayer
    GxFile = GxLayer
    'Layer file path
    GxFile.Path = "C:\Layers\Aditional Information.lyr"
    'Add the new layer to ArcMap
    Map.AddLayer(GxLayer.Layer)
    'Add the layer invisible
    GxLayer.Layer.Visible = False
    'Refresh the view
    MxDoc.ActiveView.Refresh()
End Sub


texts: How to setup leaflet map to display a WMS layer in ESRI:102012 projection?

texts:
I have a WMS server that serves layers in ESRI:102012 projection.
http://spatialreference.org/ref/esri/102012/
It accepts GetMap requests with srs=EPSG:102012 parameter.

My goal is to display a map in this projection with some WMS layers on it.

I'm using a proj4leaflet plugin (https://github.com/kartena/Proj4Leaflet) to setup a map for this projection, but I'm stuck with defining L.Transformation, scale and resolution parameters.
Can anybody tell me what to put in those parameters, or give a working example of leaflet map in this projection?


texts: Open source geometry service like ArcGIS?

texts: Does anyone know of an open source GIS server that provide similar functionality to ArcGIS' web services server? For example, the service would provide a similar geometry service that buffers feature shapes and returns the buffered shape? Or does joins or reprojections of geometry?


texts: Draw a new curve for each level curve with specific slope in QGIS

texts:
I´m new to QGIS due to a project I´m involved in.
For this project we want to draw a new curve for each level curve, but with a 0.5% positive slope for the valleys and 0.5% negative slope for the hills.
I managed to get the curves and inclination map. I also found a way to draw lines or polygons close to the curves. But I can´t find a way to perform this task on an 'automatic' way.
Could please somebody give me a clue on how to do it?


texts: How do I exclude raster layers from legend? (Silverlight, ArcGIS Server)

texts:
I'm trying to adjust my Print Widget so that raster layers do not show up in the legend. My thinking was that I could test for a layer type and not add the item to the legend if it is a raster.

Something like:

Foreach (layer in List<Layers>)
    If layer is a raster:
        do nothing;
    Else:
        Add it to the legend;

I know there's a LayerInfo thing out there but I'm not really sure how to use it. Incidentally, I'm new to C#/Silverlight/ArcGIS for Server so this simple task is a bit challenging for me. I did notice there was a GetType() method(?) that I might be able to use but I'm not sure if how if that is indeed the answer. The foreach legend item code already exists and such -- I'm looking to wrap some conditions around the legendItems.Add bit. I found the following resources, but not sure how to use them if they are helpful...
http://help.arcgis.com/en/webapi/sil...o_members.html --> the "Type" property, perhaps?

Alternatively, if I'm going about this in a difficult/wrong way, feel free to point that out! I'm a new intern with no prior C#/ArcGIS Server/Silverlight experience so it's all a bit overwhelming.

Thanks in advance!

Thread on arcgis forums: http://forums.arcgis.com/threads/83788-How-do-I-retrieve-layer-type

additional thoughts:
-- I'll somehow have to create a list of all the layers, too >_<.


texts: ESRI JS API, how to disable panning with the numeric keypad?

texts:
We have a web-based mapping application that uses the ESRI JS API. We would like to disable panning with the numeric keypad. There is a method in the API, esri.Map.disableKeyboardNavigation(), that disables panning with the dedicated arrow keys and leaves the numeric keypad enabled. We want just the opposite, i.e disable panning with the numeric keypad and leave the dedicated arrow keys enabled.

Why do we want to do this? We have a search box above the map, and if you click on the search box to give it focus and start typing, and the arrow cursor (not the vertical bar cursor, which is still in the search box) drifts back down into the map, using the numeric pad to type numbers results in panning of the map instead.

We have tried various combinations of capturing the keyup, keydown, keypress events for the map and the search box, but have not found a combination that solves our issue.

EDIT: ESRI's samples have the same behavior, so here is one of their samples in a jsfiddle: http://jsfiddle.net/Q4NVN/2/. It successfully stops the panning, but numeric input is blocked when the arrow cursor is over the map. Ideally, we should be able to enter digits whenever the vert bar cursor is in the test box, no matter where the arrow cursor is.

$("#search").keydown(function (event) {
  var x = event.which;
  if (x > 96 && x < 106 && $('#map').is(":hover")) {
    var press = jQuery.Event("keypress");
    press.ctrlKey = false;
    press.which = x - 48;
    press.keyCode = x - 48;
    press.charCode = x - 48;
    $("#search").trigger(press);
    event.preventDefault();
    return false;
  }
});

EDIT: This would appear to work, but does not take into account cases where the vert cursor is not at the end of the string (e.g. editing the search criteria):

$("#search").keydown(function (event) {
  var x = event.which;
  if (x > 96 && x < 106 && $('#map').is(":hover")) {
    $(this).val($(this).val() + (x - 96));
    event.preventDefault();
    return false;
  }
});


texts: Control drawing order of points from DBMS in MapInfo?

texts:
I have a layer of points stored in a table in a database (Oracle 11g). The table have a boolean field called Supporting. I have a thematic layer coloring the supporting points white and the non-supporting points red. Supporting points are often very close to non-supporting points. Since I am mainly interested in the non-supporting points I want to draw any red points on top of any white points. Can I accomplish this in MapInfo 11.5?

I have tried opening the table window and sorting it. That does not seem to affect the drawing order. I have also tried adding an ORDER BY column in the SQL when using "Expert mode" in the dialog "Open alternative for DBMS table". This do affect the drawing order, but not in a predictable or useful way.

(English names of options and dialogs are translated from Swedish by me, so they may not exactly match the names in the english version of the software.)


texts:
I'm trying to use tiles from an ArcGIS WMTS server with Openlayers 3 but I'm running in to a frustrating problem - optionsFromCapabilities() gives the following unhelpful error:

Uncaught (in promise) TypeError: Cannot read property 'a' of null(…)
Ie @ ol.js:68
We @ ol.js:68
Ze @ ol.js:68
(anonymous function) @ ol.js:855
(anonymous function) @ (index):18`

To keep things simple I'm just using the Openlayers WMTS Layer from Capabilities example - http://openlayers.org/en/v3.12.1/examples/wmts-layer-from-capabilities.html only changing the XML document to the layer I'd like to use.

I've verified that the XML document is being retrieved and Openlayers can parse it with WMTSCapabilities.read(). The layer and matrixSet are correctly defined. I tried using another ArcGIS server of the same version (10.22) (http://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/WMTS/1.0.0/WMTSCapabilities.xml) and that worked properly. I've thought that the problem might be projection based but it seems to be failing well before the point that projections enter in to things.

Here is my javascript code:

var parser = new ol.format.WMTSCapabilities();
var map;

fetch('http://imaps.indy.gov/arcgis/rest/services/2015Photography_Labels/MapServer/WMTS/1.0.0/WMTSCapabilities.xml').then(function(response) {
return response.text();
}).then(function(text) {
var result = parser.read(text);
var options = ol.source.WMTS.optionsFromCapabilities(result,
    {layer: '2015Photography_Labels', matrixSet: 'default028mm'});

map = new ol.Map({
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM(),
      opacity: 0.7
    }),
    new ol.layer.Tile({
      opacity: 1,
      source: new ol.source.WMTS(options)
    })
  ],
  target: 'map',
  view: new ol.View({
    center: [19412406.33, -5050500.21],
    zoom: 5
  })
});
});

Any thoughts?


texts:
I want to be able to overwrite a shapefile if it exists. I think my code needs some tweaking. I am using v10. I am able to delete the actual .shp file but the .dbf, .shx, etc still remain, so if I try to overwrite again, I get an error.

How can I remove all the files associated with the shapefile?

DoesItExist = True
geometry_type = "POLYGON"
print "Your shapefile is being generated."
out_path = raw_input("Enter in path to save shape file: ")
out_name = raw_input("Enter in name of shape file: ")

while DoesItExist == True:
    if (os.path.exists(os.path.join(out_path, out_name))):
        print "This file already exists. Do you wish to overwrite it?"
         choice = raw_input("Y/N:")
    if (choice == 'Y' or choice == 'y'):
        print "File will be over written"
        os.remove(os.path.join(out_path, out_name))
        #arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type)
        DoesItExist = False
    else:
        print "Choose new file name/location!"
        out_path = raw_input("Enter in file path: ")
        out_name = raw_input("Enter in file name: ")

else:
    print "Shapefile created in: " + str(os.path.join(out_path, out_name))
    arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type)
    DoesItExist = False


texts: QGIS not displaying dates in attribute tabl

texts:
I have a QGIS project that links to data in a PostgreSQL database. For some reason, QGIS is not displaying any data in date fields. The date fields appear to be fine in Postgres and I have performed calculations using the date fields in the database. The attribute tables in QGIS show the date field, but show the field as empty. Also, QGIS lists the date field size as -1. I'm using QGIS 2.2 and Postgres 9.0.

I've searched, but haven't seen anyone else with this problem.

Anyone have any insights?


texts: How to set arcpy code that do not automatic add to display completed layers when batch a seriels raster/shp?

texts: Maybe this is a small question.But I can't find any set method. When using arcpy to batch in arcmap, it automatic add the operated raster/shp files to new data layers. how to close it ? the set is in arcmap or in arcpy codes? I just want to carry out data using arcmap. When it loads data to a new layer, it causes routine slowly.


texts: How to detect a covered polygon by other polygon in same QGIS layer?

texts:
Identifying overlapping polygons in single layer using QGIS?

As described this question there are 3 conditions for two polygon on a same layer.

These 2 polygons can have intersection area. (How to detect-Solved)
These 2 polygons can have dublicated. (How to detect-Solved)
One polygon can cover other polygon. How to detect? (My question)

Layer: QGIS - SHP layer


texts:
I have a polygon which I know is a rectangle, although it may be at any orientation:

I want to shrink this polygon by a third of the width from two sides, so that I just get the strip in the middle - ending up with a result like the red area in this image:

I can't think how to go about doing this easily. If I try and negatively buffer the polygon then it will be shrunk on all sides. I have many tens of thousands of these polygons to do - so obviously can't do them manually.

Any ideas? I have access to ArcGIS, QGIS and am relatively experienced at Python programming - if that helps.


texts: Grouping corner-touching contiguous polygons using ArcGIS Desktop?

texts:
Using ArcGIS 10, I have a Feature Class which I converted from a Raster (GRID) (Raster To Polygon Tool).

There are approx. 180,000 polygons.
The attribute table contains a [grid_code] field which I need to preserve.
If one or more polygons are contiguous with corners touching AND they have the same [grid_code] value, then I'd like to group them (a dissolve basically) into one polygon (record).
Then I need to calculate the grouped polygon's area.

Using the Select By Location Tool, I'm able to successfully select the polygons that I need to group, but how do I then convert them into ONE polygon, while preserving the [grid_code] ?


texts: Convert hdf5 LiDAR data (from GEDI) to LAS or LAZ

texts:
I downloaded GEDI LiDAR data in hdf5 format. How can I open this file or convert it to the standard LAS format (point cloud)?

When I tried opening the hdf file with QGIS  and I got the following error:

  Invalid Layer: Raster layer Provider is not valid (provider: gdal, URL:filename)

I also tried opening it with Python:

import pandas as pd
pandas.read_hdf('filename')

I got this error

ValueError: No dataset in HDF5 file

Link to the GEDI datapool.


texts:
I want to use the LandTrendr (LT) functions in GEE and adjust them to my needs. Basically, I want to the a similar thing that this earthengine.app is doing (only NBR for a start).

I can run the analysis, compute the values. But I cannot figure out how to produce the figure.

// load the LandTrendr.js module
var ltgee = require('users/emaprlab/public:Modules/LandTrendr.js');

// define parameters
var startYear = 1985;
var endYear = 2017;
var startDay = '06-20';
var endDay = '09-20';
var aoi = ee.Geometry.Point(-122.8848, 43.7929);
var index = 'NBR';
var ftvList = ["NBR"];
var runParams = {
  maxSegments:            6,
  spikeThreshold:         0.9,
  vertexCountOvershoot:   3,
  preventOneYearRecovery: true,
  recoveryThreshold:      0.25,
  pvalThreshold:          0.05,
  bestModelProportion:    0.75,
  minObservationsNeeded:  6
};

// apply LandTrendr.js functions
var lt = ltgee.runLT(startYear, endYear, startDay, endDay, aoi, index, ftvList, runParams);
print(lt);

This should result in an array 4 x nYears dimension array. With years in (0,0,1) and calculated
values in rows 2 to 4.

Now, how to plot row 2 over years. And if I create a buffer around the aoi, do I get the values for the other pixels too?


texts: Getting QGIS 2.x python code to work in QGIS 3.x? object has no attribute 'legendInterface' problem

texts:
The following piece of code works well in QGIS 2.x, however it does not work in QGIS 3.x.

myDir = 'd:/work/output_folder/'
layers = iface.legendInterface().layers()
pipe = QgsRasterPipe()
for layer in layers:
 extent = layer.extent()
 width, height = layer.width(), layer.height()
 renderer = layer.renderer()
 provider=layer.dataProvider()
 crs = layer.crs().toWkt()
 pipe.set(provider.clone())
 pipe.set(renderer.clone())
 opts = ["COMPRESS=LZW"]
 file_writer = QgsRasterFileWriter(myDir + layer.name() + ".tif")
 file_writer.setCreateOptions(opts)
 file_writer.writeRaster(pipe,
      width,
         height,
         extent,
         layer.crs())

This is the error I get when I run the code:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
AttributeError: 'QgisInterface' object has no attribute 'legendInterface'

Does someone know what has replaced 'legendInterface' in the new version or what other changes I might have to make to get this to run in QGIS 3.0?


texts: Modify 3D numpy array using interpolation

texts:
I have a 3D array that I want to interpolate the np.nan values along the z dimension, and I just want the changes to modify my existing array. However, the changes seems not to be working. I have a test array with dimension (3,3,3) with nan values. I am accessing the z dimension and perform interpolation. I would like to have is that nan values will be replaced by interpolated values along the z dimension.

def arr_interp(array):
    arrN=np.array(array)
    ix=-np.isnan(arrN)
    idxTrue=ix.nonzero()[0]
    valueTrue=arrN[-np.isnan(arrN)]
    countFalse=np.isnan(arrN).ravel().nonzero()[0]
    arrN[np.isnan(arrN)]=np.interp(countFalse,idxTrue,valueTrue)
    # I got this function from another forum

In [177]:arrStack
Out[178]:
array([[[ 0.59016759,         nan,  0.88641936],
        [ 0.94884139,         nan,  0.89443098],
        [ 0.7730535 ,         nan,  0.86262287]],

       [[ 0.106567  ,         nan,  0.83757824],
        [ 0.21794018,         nan,  0.66128902],
        [ 0.21961069,         nan,  0.28654693]],

       [[ 0.26142901,         nan,  0.95653887],
        [ 0.19616204,  0.45      ,  0.61531475],
        [ 0.78262652,  0.76      ,  0.47786104]]])

In [181]: arrStack[0][0]
Out[182]: array([ 0.59016759,         nan,  0.88641936])

In [183]: arr_interp(arrStack[0][0])

In [184]: arrStack
Out[184]:
array([[[ 0.59016759,         nan,  0.88641936],
        [ 0.94884139,         nan,  0.89443098],
        [ 0.7730535 ,         nan,  0.86262287]],

       [[ 0.106567  ,         nan,  0.83757824],
        [ 0.21794018,         nan,  0.66128902],
        [ 0.21961069,         nan,  0.28654693]],

       [[ 0.26142901,         nan,  0.95653887],
        [ 0.19616204,  0.45      ,  0.61531475],
        [ 0.78262652,  0.76      ,  0.47786104]]])


texts:
I'm trying to create a GeoJSON REST service from a PostgreSQL connection.
First I fetch GeoJSON as a string from PostgreSQL:

import java.util.*;
import java.sql.Connection;
import java.sql.DriverManager;
import java.sql.SQLException;
import java.sql.ResultSet;
import java.sql.Statement;
import java.util.logging.Level;
import java.util.logging.Logger;
import java.sql.PreparedStatement;

/**
 *
 * @author me
 */
public class App2 {

    /**
     * @param args the command line arguments
    */

    private final String url = "jdbc:postgresql://127.0.0.1:5444/mydb";
    private final String user = "myuser";
    private final String password = "mypass";

    /**
     * Connect to the PostgreSQL database
     *
     * @return a Connection object
     */
    public String connect() {
        Connection conn = null;
        String geojson = "abc";
        try (Connection con = DriverManager.getConnection(url, user, password);
                PreparedStatement pst = con.prepareStatement("select ST_AsGeoJSON(geom) from (select ST_MakeLine(geom) as geom from (SELECT ST_MakePoint(location_longitude,location_latitude) as geom FROM locations where location_fix_time between '2019-05-13 11:08:30' and '2019-05-13 11:08:40' order by location_fix_time) a) b;");
                ResultSet rs = pst.executeQuery()) {

            while (rs.next()) {
                geojson = rs.getString(1);
                return geojson;
            }

        }  catch (SQLException e) {
            System.out.println(e.getMessage());
        }

        return geojson;

    }

    /**
     * @param args the command line arguments
     */
    public static void main(String[] args) {
        App2 app2 = new App2();
        app2.connect();
        String abc = app2.connect();
        System.out.println(abc);
    }

}

Then I'm trying to place it as GeoJSON on a REST webservice:

@Path("Mypath")
public class MyPathResource {
    @GET
    @Produces("text/plain")
    public String getText() {
        return "My First RESTful Service....";
    }
    private static GeometryJSON creategeoJSON(boolean prettyPrint) {
        App2 app2 = new App2();
        String queryres = app2.connect();
        GeometryJSON g = new GeometryJSON();
        try {
            g.readLine(queryres);
            return g;
        }
        catch (IOException e) {
            System.out.println(e.getMessage());
        }
        return g;
    }

    @GET
    @Path("przejazd/{idpoj}")
    @Produces(MediaType.APPLICATION_JSON)
    //@Produces(MediaType.TEXT_PLAIN)
    public Response przejazdfetch(@PathParam("idpoj") int idpoj)  {
    //JSONObject jsono = createJSON(true);
        GeometryJSON gjson = creategeoJSON(true);

    //return Response.ok().entity(jsono.toString()).build();
        return Response.ok().entity(gjson.toString()).build();
    }
}

This doesn't work however and when I try to load the webservice at http://localhost:8080/MyFirstREST/webresources/Mypath/przejazd/1 I get an error:

  SyntaxError: JSON.parse: unexpected character at line 1 column 1 of
  the JSON data.

Any ideas how it can be fixed?


texts: For a point data set with temperature data I have, I want to create their Thiessen Polygons in R. So far, the only algorithm I have found is that of Carson Farmer (http://carsonfarmer.com/2009/09/voronoi-polygons-with-r/). I would like to know if there is a ready to use function that creates Thiessn Polygons?


texts: ArcGIS Javascript API using CMV - Tiled Map Service Popup

texts:
Using ArcGIS Online webmaps, you can enable a popup for a tiled map service by selecting a feature layer either currently in the map, or providing a url to the feature layer service. I believe then the popup attributes for the tiled layer use the feature layer's attributes. Is there a way to do this using Javascript API? I am currently testing the CMV application but I need the popup option for a tiled layer. Here is the layer code in the CMV viewer.js file.

{
 type: 'tiled',
url: 'http://tiles.arcgis.com/tiles/cWzdqIyxbijuhPLw/arcgis/rest/services/Antelope_Areas/MapServer',
title: 'Antelope Hunt Areas',
options: {
          id: 'antelopeareas',
          opacity: 1.0,
          visible: true,
         }
}


texts:
I am doing a project that involves Web Mapping.

The project involves having a simulated communications network for a wind turbine company.

My project involves building a user front end to display the health of each turbine on the map.

I will be using using Simple Network Management Protocol (SNMP) to poll Management Information Bases (MIBs) for performance, and a few databases to enrich the data that will be then displayed on the map.

I am trying to understand the most efficient protocol i can use to display the info that is stored in the monitoring databases (e.g. HTTP, SOAP, REST, JSON, OSLC, etc) as displaying this info quickly is vital.

Would anyone have any suggestions on which to use for such circumstances?


texts: Using Tasseled Cap in ArcGIS Desktop on Landsat TM 4-5 .tif files?

texts:
I am using ArcGIS Desktop 10.1.

I am attempting to conduct a tasscap Transformation via "Raster Property Functions" on a composite file that I have created from unprocessed, raw data files (meaning tiff-files Bands 1-5 and Band7) via the "composit bands"-tool from Landsat TM 4-5.

The error message occuring when I try to apply the function is as follows:

  The function could not be inserted. The Input or Output arguments are
  invalid.


texts: How to reference a PostGIS JSON field in GeoServer SLD

texts:
I have a PostGIS table with a JSON field which I successfully published in GeoServer 2.16
The contents of the JSON field (named "display_properties") are an array and like this (some values have been omitted):

[{"name" : {"en" : "Level", "el" : ""}, "value" : {"en" : 1, "el" : 1}}, ... ]

I would like to make an SLD to symbolize the layer with an attribute-based rule from the contents of the JSON file. Is there a chance I might use something like this?

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
                       xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd"
                       xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
                       xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

  <NamedLayer>
    <Name>test_JSON_field</Name>
    <UserStyle>
      <Title>A cyan line style</Title>
      <FeatureTypeStyle>
        <Rule>
          <Title>cyan line</Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>display_properties[0]["value"]["en"]</ogc:PropertyName>
              <ogc:Literal>1</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <LineSymbolizer>
            <Stroke>
              <CssParameter name="stroke">#0099cc</CssParameter>
            </Stroke>
          </LineSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

P.S. There's an older similar question (GeoServer and PostgreSQL JSON column) but not about referencing a JSON field in SLD.


texts: Adding information from polygon to raster in R?

texts:
I would like to add information from a spatial polygon to a new raster.

r <- raster(nrow=10, ncol=10)
r[] <- runif(ncell(r)) * 10
xym = matrix(c(0,0,0,45,45,45,45,0),4,2,byrow=T) ## create coordinate list
p = Polygon(xym)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))
plot(r)
plot(sps,add=T,col="blue")

How do I create a new raster that assigns the numeric value of the polygon (FID) depending on if the polygon overlaps the raster?


texts: WMS layer not displaying in OpenLayers 3

texts:
Three of my WMS layers are not displaying on the map. ......

  trans = new ol.layer.Tile({
     name: 'Transmission Lines',
     visible : false,
     source: new ol.source.TileWMS({
         url: 'http://*******/cgi-bin/mapserv?map=/data/mapFile/trans.map',
         params: {
             'LAYERS': 'trans',
             'TRANSPARENT': true
         }
     })
  });

  energyGroup = new ol.layer.Group({
     layers: [trans, distribute, station],
     name: 'Energy'
  });

  view = new ol.View({
     center: ol.proj.fromLonLat([15.4542, -10.7587]),
     zoom: 3,
     maxZoom: 14,
     minZoom: 3
  });

  olMap = new ol.Map({
        layers: [******, energyGroup, ******],
        view: view
  });

Here is my mapfile: ...

EXTENT 14.512797 -34.507545 35.797839 -23.569398
FONTSET 'fonts.txt'
SYMBOLSET 'symbols.txt'
PROJECTION
        'proj=longlat'
        'ellps=WGS84'
        'datum=WGS84'
        'no_defs'
END

IMAGECOLOR 255 255 255
IMAGEQUALITY 95
IMAGETYPE png

OUTPUTFORMAT
        NAME png
        DRIVER 'GD/PNG'
        MIMETYPE 'image/png'
        IMAGEMODE RGBA
        EXTENSION 'png'
END
# Legend
LEGEND
    #IMAGECOLOR 255 255 255
    STATUS ON
        KEYSIZE 18 12
        LABEL
            TYPE BITMAP
            SIZE MEDIUM
            COLOR 0 0 89
        END
END

WEB
    TEMPLATE '/data/mapFile/template.html'
    IMAGEPATH '/data/mapFile/output/'
    IMAGEURL 'data/mapFile/output/'

    METADATA
        'wms_title'             'Transmission Lines'
        'wms_onlineresource'    'http://*****/cgi-bin/mapserv?map=/data/mapFile/trans.map&'
        'wms_enable_request'    '*'
        "wms_feature_info_mime_type" "text/html"
        'wms_srs'               'EPSG:4326'
        "wms_server_version"    "1.1.0"
    END
END

LAYER
    NAME 'trans'
    TYPE LINE
    DATA '/******/transmission.shp'
    STATUS ON
    DUMP TRUE
    METADATA
        "gml_include_items"     "all"
        "wms_include_items"     "all"
        "wms_feature_info_mime_type" "text/html"
        "wms_enable_request"    "*"
    END
    TRANSPARENCY 100
    PROJECTION
        'proj=longlat'
        'ellps=WGS84'
        'datum=WGS84'
        'no_defs'
    END
    CLASS
            NAME    'Transmission Lines'
            STYLE
                WIDTH   0.90
                OUTLINECOLOR 203 69 51
            END
    END
END

What am I missing here?


texts:
How rotate map using the leaflet? I was trying here, but nothing worked. I followed this tutorial, even though it was from the mapbox, I assumed it was useful:

https://www.mapbox.com/mapbox-gl-js/example/navigation/

The error is:

Uncaught TypeError: t.addTo is not a function
    at e.addControl (leaflet.js:5)
    at HTMLDocument.<anonymous> (index.html:85)

And in my code, I just add

<script src='https://api.tiles.mapbox.com/mapbox-gl-js/v0.42.2/mapbox-gl.js'></script>

And

map.addControl(new mapboxgl.NavigationControl());

I trying too http://labs.easyblog.it/maps/leaflet-compass, but it cant be use for rotate map. How do I rotate the map?


texts: Cannot display attachments in ArcGIS Pro

texts:
I'm trying to learn ArcGIS Pro 1.4, and with a little bit of pain I was progressing till today in which I cannot make my previous hyperlinks get ready in ArcGIS Pro. I think enabling attachments in first place and configure pop-ups afterwards might the be the right workflow to follow, correct?
I have reviewed the post: Work with file attachments, but it hasn't helped me so far with this issue.
These are the steps I have followed so far:
i.e.: I have a point feature class, with some of the locations having associated files which I want to convert into attachments.

Analysis tab -> Tools -> Geoprocessing pane: Enable Attachments on the feature class that I want it to have hyperlinks functionality enabled.

This operation went well:

After doing this a relationship class and a table have been created in my geodatabase, OK.

Then, If I go the Attributes panel, Attachments, click +Add, when I pick up the file, it pop-ups an error:

What am I missing?


texts: can we use SetSpatialFilter function in loop with new filter geometry on same layer in every iteration

texts:
I am trying to find shortest path between two point feature. There are many links between two points hence first i make buffer around starting point and select all feature inside that buffer. After this I select shortest link among the selected feature. Then I go to the last point of the selected linestring and make buffer around that point to select the linestring ahead of that point. I keep doing this in a loop till i reach the last point so as to obtain a sequence of link between two point but every time when i make new buffer and pass it as argument in the SetSpatialFilter function it give me the same feature as for the previous buffer.

while x:
 datasource=open_File(path="C:\Users\Hp\Desktop\ABSTRACT\SHASHANK\SHASHANK.shp")
    layer = datasource.GetLayerByIndex(0)
    layer.SetSpatialFilter(buffer)
    for feature in layer:
        featureindex += 1
        route[featureindex] = feature.GetGeometryRef().ExportToWkt()
    buffer=MakeBuffer()


texts: How are Class Extensions associated with a Feature Class in ArcGIS v10?

texts:
In ArcGIS 9.x if a Feature Class is registered with a Class Extension the CLSID of the extension is recorded in a table in the sde schema called sde.GDB_OBJECTCLASSES (EXTCLSID field).

However in ArcGIS 10 there is no such table, and I can't find anything comparable in an SDE 10 database. I assume this information must now be stored in the data server's registry, but I am not sure where to start.

I am trying to debug issues apparently caused by an upgraded (ArcGIS v10) Class Extension updating data in a Feature Class originally registered with a v9.3 Class Extension and held in a v9.3 SDE database (both versions of the extension share the same GUID), so any insight on this would be much appreciated.


texts: Creating symbology for multiple attributes within the same feature class in ArcMap 10

texts: I'm struggling to find a way to take a single point feature class that has multiple attributes that I want as feature layers with unique symbology.  Let's call these attributes Depth1 and Depth2.  A key point here is I want to apply a lyr setting to automate the symbology of each layer and it seems what may be the best way of doing this is in the properties of an output in model builder you can add a lyr file to use as symbology.  But I'm not sure how I can get these two attributes as new layers using model builder and then default to correct attribute for the symbology for each (ie. Depth1 for first layer & Depth2 for the 2nd layer).  Any ideas?  Let me know if I need to explain any of this better.


texts: Getting latitude and longitude of vertices/corner points of polygon/building from online maps such as OpenStreetMap?

texts:
I have tried using JOSM(Java OpenStreetMap) and QGIS for directly getting latitude and longitude of vertices/corner points of polygon/building but I am unable to get.

Is there any way to get?


texts: Can we use ArcGIS (sampleserver1.arcgisonline.com) API to get Households for area without any license?

texts:
In my web site I am using "http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer" API for getting households for an area on Google map.

Is this API is free? or did I use this without license?


texts: Remove polyline vertex outside polygon in Arcpy

texts:
I am trying to remove polyline vertex which are outside the buffer polygon and not creates new vertex on intersection.

I found a script which remove vertex but inside polygon. What is more it creates new vertex on intersection. It may help with this task.

import itertools, arcpy
arcpy.env.overwriteOutput = True
line_FC = arcpy.GetParameterAsText(0)
poly_FC = arcpy.GetParameterAsText(1)
output = arcpy.GetParameterAsText(2)

#clone input lines to memory for fast processing
temp = arcpy.CopyFeatures_management(line_FC, "in_memory/temp")

polygons = arcpy.CopyFeatures_management(poly_FC, arcpy.Geometry())

for poly in polygons:
    with arcpy.da.UpdateCursor(temp, ["SHAPE@"]) as uCursor:
        for line in uCursor:
            line = line[0]
            diff = line.difference(poly)
            #if the two lines are not equal, that means it intersected the polygon
            if not line.equals(diff):
                #the result of geometry.difference() is a multipart line of only those
                #parts that lie outside the polyon
                parts = diff.getPart()

                #if parts is empty that means the line is completely within the polygon
                #i.e., no difference
                if parts:
                    #We'll need to "join" the end of part1 to the beginning of part2
                    #so we'll just flatten the list of lists
                    joined = list(itertools.chain.from_iterable(parts))

                    #and create a new polyline object to update the shape
                    poly_trimmed = arcpy.Polyline(arcpy.Array(joined))
                    uCursor.updateRow([poly_trimmed])

arcpy.CopyFeatures_management(temp, output)

Is there a way to do it in arcpy? I am using ArcGIS Pro 1.3.


texts:
I have a set of  polygons out of which I want to make a KML file, my table schema is as follows :

CREATE TABLE [dbo].[utilitygeo](
    [ID] [int] NOT NULL,
    [PolygonID] [varchar](8000) NULL,
    [geog] [geography] NULL
)

I tried extracting Lat/Long from each polygon and then thought of converting them into a KML file, following is the code:

declare @g geography;

select top 1 @g =  geog from MyDB..Geoterritory

declare @i int = 1
create table ##result(
    Latitude float,
    Longitude float
)

while @i < @g.STNumPoints()
BEGIN
    insert into ##result values(
        @g.STPointN(@i).Lat,
        @g.STPointN(@i).Long
    )
END

select * from ##result

But this takes a lot of time.
Is there a better way to generate KML file from these polygons in SSMS.
Thanks


texts: OpenLayers feature of featurecollection removal issue - how to force map refresh?

texts:
I have a Features collection loaded as JSON from PostGIS, it is represent a set of vector areas over map (units responsibility alignment). I need to create new, modify and delete existing zones already loaded.

I have started with displaying existing areas and got feature removal BUT since first step is to select it with mouse click and than press "erase" icon, a feature remain in place until I click somewhere outside selected area. So basically saying it is removed from the feature collection successfully, but no map refresh events occurs until manual outside click happens.

I tried a lot of different ways but without any success. Does anyone was successful in resolving similar issue?

Some piece of code:

var selectedFeatureID;

function addSelect() {
    map.removeInteraction(draw);
    singleClick = new ol.interaction.Select();
    map.addInteraction(singleClick);

    singleClick.getFeatures().on('add', function (event) {
       var properties = event.element.getProperties();
       selectedFeatureID = properties.id;
    });
}

function removeSelectedFeature() {
    var features = vectorSource.getFeatures();
    $(features).each(function (i) {
      if (features[i].getProperties().id == selectedFeatureID) {
        vectorSource.removeFeature(features[i]);

      }
    });
}


texts: How to make distance calculated using Latitude and Longitude more accurate

texts:
I'm coding a Web App whereby users can choose their city (start) and a destination city, and find out the approximate distance between both cities, using Latitude and Longitude of both cities to calcultate. Then i thought that since this will be done using the concept of a Straight Line on a sphere, it means that i'll be getting a result of a straight line which is the shortest, but in reality (driving route especially), will have hills, curves, turns. It might be practically impossible to have an accurate distance, but how can i make the app more accurate that just a straight line distance.

I'll be most delighted by you opinions, suggestions, answers and comments.
Thanks in advance.


texts:
I am using the OSM Boundaries Map 4.0 to download the boundaries of the countries. The problem is that I want to include the disputed boundaries for countries such as Bhutan:

Is there anyway of downloading the disputed boundaries from OSM or any other accurate vendor?


texts:
What is the main difference between a

  Web Map Tile Service (WMTS)

and a

  Web Map Service (WMS)

When should you choose one over the other?

What are the (dis)advantages by choosing one over the other?

I have read the Wikipedia entries.


texts:
I have a problem converting my Raster Catalog to Raster Dataset. I have an error (see below) after 4-5 hours of process - all rasters have ~21Gb. Is there something wrong with data?

I will be really grateful for your help :)

==> Sure, here it is:


texts: I downloaded release ol3-cesium on my computer from https://github.com/openlayers/ol3-cesium/releases/ and when I open examples from the folder it doesn't show 3D globe or vector data. When I open same examples from http://openlayers.org/ol3-cesium/examples/ it's work properly. Also I try to load vector data and I can't see anything.


texts: Problem with adding new layers with EPSG:404000 in Geoserver

texts:
I'm using a Postgres database (with a PostGIS extension) connection in Geoserver to add new layers. When I add a new layer from this connection, in the Coordinate Reference Systems Section, Geoserver cannot identify correct SRS and the Declared SRS in EPSG:404000.

This error appears with most of my layers. Every layer has the same EPSG 4326 and I can confirm it by opening them in QGIS or consult it via query in PgAdmin. At first, I thought that the difference was the import tool: some were imported to the Postgres database using QGIS and worked, while the others, using the PostGIS Shapefile Importer, didn't. Later on, I tried to import the error layers using QGIS, but when I added them as a new layer in Geoserver, I still got the same error: Declared SRS in EPSG:404000.

What am I doing wrong? Is this a problem with the configuration of Geoserver, Postgis or the Shapefile SRS data?

Postgres Database is in Azure.


texts:
I'm creating a shapefile with a field type of real, but I want to lower the standard precision from 15 decimal places to 10, here is my code:

dst_filename = "Polygonized.shp"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_filename )
dst_layer = dst_ds.CreateLayer("layer", geom_type=ogr.wkbPolygon )

fd = ogr.FieldDefn( str(Polygonize_Field_Name), ogr.OFTReal )

fd.SetPrecision(10)
dst_layer.CreateField( fd )
gdal.FPolygonize( band,band1, dst_layer,0, [], callback=None )
dst_ds = None

(I'm also polygonizing, but that works fine.)

When I go to open up the file and return the precision using:

ds = ogr.Open("Polygonized.shp")
lyr = ds.GetLayer()
lyr_dfn = lyr.GetLayerDefn()
lyr_fld = lyr_dfn.GetFieldDefn(0)
lyr_fld.GetPrecision()

It clearly returns a precision of 15.

Am I doing something wrong when specifying the precision?


texts: Exporting Google basemap as .tif (non-commercial) w/QuickMapServices

texts:
I'm trying to export a Google satellite basemap as  tif to visualise a personal  project I'm working on.

I'm using QuickMapServices to fetch the Google basemap. But whenever I use "Save As" I get the error "Unsupported data source". I wondered if it was a coordinate issue, but its set to WGS84/Pseudo Mercator.

Up in Canada, it seems that my only basemap options are Bing and Google. The problem with Bing is that its too blurry to make out some features I'm trying to identify.


texts: Populate GeoExt 2 gridpanel from OpenLayers vector features

texts:
I need to populate a ExtJS datagrid with a OpenLayers features. The official example http://geoext.github.io/geoext2/examples/grid/feature-grid.html doesn't help me, because the data is loaded from a JSON file, but I need the features parsed from a WFS layer.

Could anyone help me with a source code example (jsfiddle)?


texts: QGIS reshape features to close indents in polygons

texts:
Have already run the buffer tool and do not want to make the polygons any larger with another buffer run.

Reshape features tool seems like it should work I am just not seeing how to operate it after many attempts and online searches.

Screenshot is  , showing how these large indents need to be eliminated, closed, or whatever the term is called?

Arcgis has an editor tool that does this but I need to use QGIS.


texts: How to preserve the order of fields when exporting a table from ArcMap?

texts:
ArcMap allows you to re-order fields in an attribute table, by dragging+dropping either within the layer's > Properties > Fields tab, or within the Table window.

However, the order of the fields doesn't seem to be honoured when exporting this table to a new table.

Is it possible to preserve the order of fields when exporting an attribute table?


texts: Calculating statistics per area for categorical raster using R?

texts:
I am working with the area of New York City (NYC) in ArcGIS 10.1. I have a landcover raster map with several different land cover categories. What I am trying to do is to build up a table in which, per borough, I have the percentage occupied per land cover category.

My first approach was to reclassify the raster, setting every category to NoData, keeping only one category as 1. The next step would be to perform a zonal statistics as table that provides me with the count of pixels per borough. Knowing the resolution of my raster, I can then infer the area of the borough covered by that category and then divide it by the borough's total area.

However, this process seems quite repetitive and slow, generating tons of intermediate data (up to 7 reclassify runs).

Can you think of a way to do this faster using R?


texts:
Possible Duplicate:
  Steps to Start Web Mapping

soon I will be trusted to develop a webgis system. This project will take an advantage of a currently developed business system. The business system is using pure PHP, Java and MySQL. The project plan is to embedded a mapping system so the management will be able to see all the report tables (from the business system) and when they click on some "numbers" they can view the map as well. I guess it's a standar nowadays.

What I'm hoping to get answer is what free webgis environment that you are suggesting to accomplished this? I'm not new in webgis but I'm not following the recently updated information about this webgis environment. Mostly I'm using desktop GIS.

Thank you.


texts:
I have level 2 of NOAA nc file. The data is in https://1drv.ms/u/s!AjXAWS5maHffhyqaxqd61ok-SFsY?e=mcUdnE.
I want to open the file in QGIS or GRASS GIS. But, they don't recognize it. So, I use GDAL.
I get it from gdalinfo :
$ gdalinfo noaa19_L2_d20210102_t2309066_e2319520_c20210102124945585006_iapp.nc
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: noaa19_L2_d20210102_t2309066_e2319520_c20210102124945585006_iapp.nc
Size is 512, 512
Coordinate System is `'
Metadata:
NC_GLOBAL#End_Scan_Line_DOY=0
NC_GLOBAL#End_Scan_Line_Time=0
NC_GLOBAL#End_Scan_Line_Year=0
NC_GLOBAL#IAPP_Version_Number=4
NC_GLOBAL#Satellite_Name=noaa19
NC_GLOBAL#Start_Scan_Line_DOY=0
NC_GLOBAL#Start_Scan_Line_Time=0
NC_GLOBAL#Start_Scan_Line_Year=0
Subdatasets:
SUBDATASET_1_NAME=NETCDF:"noaa19_L2_d20210102_t2309066_e2319520_c20210102124945585006_iapp.nc":Latitude
SUBDATASET_1_DESC=[0x18] Latitude (32-bit floating-point)
...
...
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

Then I want to translate to GTiff. But, it is error.
gdal_translate -of GTiff NETCDF:"noaa19_L2_d20210102_t2309066_e2319520_c20210102124945585006_iapp.nc":Surface_Temperature test.tiff
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Input file size is 18, 0
ERROR 1: Error: -srcwin 0 0 18 0 has negative width and/or height.

Then I try to use gdalwarp. But, it is still error.
gdalwarp -t_srs EPSG:4326 NETCDF:"noaa19_L2_d20210102_t2309066_e2319520_c20210102124945585006_iapp.nc":Surface_Temperature test.tiff
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
ERROR 1: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for NETCDF:noaa19_L2_d20210102_t2309066_e2319520_c20210102124945585006_iapp.nc:Surface_Temperature. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

I also try SeaDAS. The message : No appropriate reader found.
I also try R. It is from https://www.researchgate.net/post/What_is_the_best_way_to_convert_nc_files_into_tif_files. It is error.
library(ncdf4)
library(raster)
fname <- file.choose()
nc<-nc_open(fname)
PCP <- brick(fname, varname="Surface_Temperature")
writeRaster(PCP, filename="sst.tif", format="GTiff",overwrite=TRUE)

Error in the step:
PCP <- brick(fname, varname="Surface_Temperature")
[1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named Across_Track BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
[1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named Along_Track BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
Error in .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
cells are not equally spaced; you should extract values as points
In addition: Warning messages:
1: In min(rs) : no non-missing arguments to min; returning Inf
2: In max(rs) : no non-missing arguments to max; returning -Inf
3: In min(rs) : no non-missing arguments to min; returning Inf

So, how to export to GTiff, so I can open in GRASS or QGIS ?


texts:
I performed a "join" going to the properties option in a shapefile (in this case is just a table), in order to join the attributes from another shapefile that has a field in common, specifically "OBJECT_ID".

Some features appear with the right numbers and others with "Null".

However, in the shapefile I wish to join the information is correct and there is no "Null" in the trable.

Any idea why this happens?


texts: Grib file information to create QGIS projection

texts:
I'm trying to create a specific projection with QGIS
(something like that: +proj=lcc +lat_1=69.75 +lat_2=52.87 +lat_0=61.31 +lon_0=7.665 +x_0=0 +y_0=720 +ellps=GRS80 +units=m +no_defs).
I have some grib files which contain these informations:

Input geometry:

Data comes from:         233         255          38
Projection     :           3
Nlon,Nlat,Nlev :         720         800          65
South,West     :   52.870000000000005        7.6650000000000000
North,East     :   69.756418673948701        45.483402182515576
Dlon,Dlat      :   2500.0000000000000        2500.0000000000000
Projlon,Projlat(1,2)    :   20.000000000000000 62.399999999999999        62.399999999999999

Scan: WE:SN winds(grid)
Lambert Conf: Lat1 52.870000 Lon1 7.665000 Lov 20.000000
Latin1 62.400000 Latin2 62.400000 LatSP 0.000000 LonSP 0.000000
North Pole (720 x 800) Dx 2.500000 Dy 2.500000

And I have my image (800 x 720) with top left corner at (0,0). The aim of this personnal projection is to put images in Lat Lon in WGS84.

How can I create this projection?


texts: 35km difference between Google Earth and Map coordinates when manual input

texts:
I have GPS coordinates. If i enter them manually in the Google Maps Engine, there is a 35 km distance between the points and the same coordinate in Google Earth.

However, if i import the Earth's KMZ into MAPS, they are located correctly but cannot see the coordinate.

Is there some behind the scenes conversion made by Maps?


texts: Data Frame Property options have changed

texts: I have Googled and searched ESRI help files and experimented until I am blue, not to mention reloading ArcGIS Desktop10 for Repair and adding the SP5 to it. I no longer have the options of border or background for data frame in layout view. I used to -- and now don't. I don't know what changed it yet it is impossible for me to combine data frame layers with the data frame borders visible. This is what I see when I right click on the data frame properties Both on the data frame in the layout view and on the data frame in the TOC in layout view:   What happened and how do I fix it? Thanks.


texts:
using arcmap 10.3

Trying to create a 'live' draw of lines between two coordinate points. So far I have been able to easily create the lines into a new shapefile using the XY to Line tool however I would like to make this layer only display and not create and actual shapefile each time.

So each time the SQL query is pushed, the connecting lines are updated.

This question has been previously asked but I am not familiar with ModelBuilder or python (willing to learn with guidance).


texts: Make toolbar icon visible in the target machine without changing the path

texts:
I had done a simple standalone application that loads a shapefile and does some basic actions like zoom in, zoom out, pan etc in the viewer. In my Python file, I had set the icon path for all these (zoom in, zoom out, pan) as directory which exist in my desktop. Now, when I send the application to the user, he is unable to view the toolbar icons. The only way is to change the icon path of the target machine and run the file or is there any other way to solve this?

# Create the actions for our tools and connect each to the appropriate
    # method
    self.actionAddLayer = QAction(QIcon("....../Standalone_External_QGIS_Application/AddVectorlayer.png"), \
                                  "Add Layer", self.frame)
    self.connect(self.actionAddLayer, SIGNAL("activated()"),self.addLayer)
    self.actionZoomIn =QAction(QIcon("....../Standalone_External_QGIS_Application/Zoomin.png"), \
                               "Zoom In", self.frame)
    self.connect(self.actionZoomIn, SIGNAL("activated()"), self.zoomIn)
    self.actionZoomOut = QAction(QIcon("....../Standalone_External_QGIS_Application/Zoomout.png"), \
                                 "Zoom Out", self.frame)
    self.connect(self.actionZoomOut, SIGNAL("activated()"), self.zoomOut)
    self.actionPan = QAction(QIcon("....../Standalone_External_QGIS_Application/pan.png"), \
                             "Pan", self.frame)
    self.connect(self.actionPan, SIGNAL("activated()"), self.pan)
    self.actionExplore = QAction(QIcon("....../Standalone_External_QGIS_Application/Identity.png"), \
                                 "Explore", self.frame)
    self.connect(self.actionExplore, SIGNAL("activated()"), self.Explore)
    self.actionZoomFull = QAction(QIcon("....../Standalone_External_QGIS_Application/ZoomtoFull.png"), \
                                 "Zoom Full Extent", self.frame)
    self.connect(self.actionZoomFull, SIGNAL("activated()"), self.zoomFull)


texts:
I am developing a script to calculate the values of two fields (C and D) based on value of field A (It is a Point feature):

A C D
0.2 1 1
0.2 1 2
0.2 1 3
0.4 2 1
0.5 3 1
0.5 3 2

So the idea is

n = 1
k =0
while i < endrow
{
if a[i]==a[i+1]
c[i]=n
d[i] = k+1
if a[i] < a[i+1]
i=i+1
k=0
c[i]=n+1
d[i]= k+1
}

Not sure - cursor can help me. But is there any way for such comparison and update the column?

Thanks

adding another complication- field E which has the maximum D for every unique A or C; so it will be:

 A C D  E
0.2 1 1 3
0.2 1 2 3
0.2 1 3 3
0.4 2 1 1
0.5 3 1 2
0.5 3 2 2

Because my ultimately goal is to keep only D ==1 and label the remained points with C and E (say: there are 3 points at location 1)

A C D  E
0.2 1 1 3
0.4 2 1 1
0.5 3 1 2

Right now appears, I have to run another combination of search and update cursor after calculation of C and D but if there is any other concise and sleek approach! If I will figure it out - I will post it here!


texts: Accessing HTML using jQuery within LeafletJS popup.

texts:
So I want to add some HTML code to a Leaflet popup, but I also want to be able to to perform some jQuery functions when this HTML is clicked, like so:

.bindPopup( "<button id='clickHere'> <strong> stuff </strong> </button>" ))

But I'm unable to access this using jQuery selectors, because it's inside  tags.

Is there any way to solve this problem?


texts: Generating series of ellipses on line using ArcGIS Desktop?

texts:
I am trying to create a series of ellipses (mocked up in green below) of defined dimensions along a line. I have got as far as generating regularly spaced points along the line that are attributed with the major (axis perpendicular to line direction) radius r1 and minor (axis parallel to line direction) radius r2.

The aim here is to build a model so that I can go back and easily create new polygons for tweaked values of r1, r2 and varying point geometries (points will always be on a line, and r1, r2 will always have same geometric relationship to line.

I have access to ArcMap and extensions.


texts: Building a Definition Query using Sub-queries with multiple conditions

texts:
Background

Since layer selections can't be be published as a map service, I need to duplicate N layers and write definition queries for each to satisfy my needs.  I can't create new standalone layers because of multiple relationships this feature class participates in (multiple M-N and 1-M relationships) otherwise I would have created a python script to take care of this.  Here is a flowchart of my tables and feature classes.  Tables/features classes are stored in a file geo-database (ArcGIS 10.4) for this particular example.

My question will focus on selecting features from the Parcel feature class based on records in the Activity Record ID Table.

Question

I’m trying to create a definition query using sub-queries with multiple conditions but running into some issues and not sure if SQL statment approach is correct.

I can create a subquery that selects parcels based on the parcel ID's in the activity record ID Table table class.  This basically selects all the parcels with Parcel ID's that are present in the Activity Record ID Table.

I want to refine the sub-query further by only selecting parcels from the activity table for a specific activity in that table. This is what I have so far but arcmap crashes as soon as apply/test the query.

I'm not sure if the SQL statement is correct but there has to be a way to create sub queries for with multiple conditions.


texts:
Using ArcGIS Javsacript API on one server to grab a PDF from a custom print service published on a different server. We have used the exact same script successfully on other apps/servers before with no problems. On a new server, however, we've run in to CORS issues, so we've tried both accepted approaches to dealing with this:

(1) Used working proxy (confirmed with other geoprocessing service in same app)

(2) esri.config.defaults.io.corsEnabledServers.push( [www.server.here] )

When using the first method (1), we initially get an esriJobFailed response. If we hit fire button again, it seems to submit the job as normal... seems strange.

If we use the second method, it submits the job fine.

At this point we are at the same step... PrintTask has been submitted and our web app is waiting for the response.

But we keep getting an "Unable to complete operation." error. When we look closely, it seems as if the PrintTask has successfully ran, but the result is not being delivered back to our app properly.

This screen shot shows a status 200 reponse with an "esriJobSucceeded" value:

But if we look at the next request, it shows that we didn't make it to the finish line:

Neither of the callbacks functions for the PrintTask (success, error) fire.

But if we listen for each xhr request and find the one with the "esriJobSucceeded" value, we can construct the correct URL and load our PDF output:

var url = "http://server.com/arcgis/rest/directories/arcgisjobs/folder/printservice/" + resp.jobId + "/scratch/web_map.pdf";
window.open(url);

This works, but it's obviously not what we want. Any thoughts on this?

UPDATE

Pretty basic callbacks. But the script doesn't even reach them.

// PRINT SERVICE
gp1 = new PrintTask("http://server.com/arcgis/rest/services/App_Print_Service", {async: true});
gp1_params = new PrintParameters();
gp1_params.map = map;

function gp1_function() {
    gp1.execute(gp1_params, gp1_complete, gp1_error);
}

// GP Success
function gp1_complete(result) {
    window.open(result.url);
}

// GP Error
function gp1_error(error_info) {
    console.log(error_info);
    alert('Unexpected error');
}


texts:
The concept of composite hillshade is to combine multiple hillshades with different sunlight orientation to avoid blind spot. It requires to merge 3 hillshades images with 315⁰ (NW sunlight, default layer), +355⁰ layer, +275⁰ sunlights, respectively. This increase details and elegance, as shown below.A full non-GIS description is here.

Using any srtm file, the 3 initial hillshades can be done via :

gdaldem hillshade input.tif hillshades_A.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges
gdaldem hillshade input.tif hillshades_B.tmp.tif -s 111120 -z 5 -az 355 -alt 60 -compute_edges
gdaldem hillshade input.tif hillshades_C.tmp.tif -s 111120 -z 5 -az 275 -alt 60 -compute_edges


texts: Running Python script in GRASS every time map extent is altered

texts:
I am looking to develop a python script that triggers every time the user zooms or pans to a new map extent. The script I have in mind will clip a raster to the visible extent thus reclassifying the colour ramp and hopefully revealing more subtle anomalies in, for example, a raster of elevation data.

My question is how to go about 'triggering' the script when the view extent is altered by the user. My intention was to access a variable that describes the extent of the currently visible map display in GRASS 7.x. A while loop could be executed that would continuously check this variable for changes and trigger the rest of the script.

I thought the variable could be accessed by the;

d.info -g

command, however this appears to be returning the co-ordinates of the current grass region as opposed to the visible extent of the map.

Problem: I cannot access the visible extent of the map in the monitor display
Reproduce: Try to create a python object that contains the visible extent of the map in the monitor display.

After launching grass with the wxGUI I have so far tried launching and selecting a new monitor using the following commands:

import grass.script as grass
grass.run_command("d.mon", start='wx0')
grass.run_command("d.mon", select='wx0')
grass.run_command("d.rast", map='elevation')

This open as new map display window and displays an elevation raster. Then I retrieve some information re the the new display using;

info = grass.read_command("d.info", flags='r')
print(info)
>> screen rectangle: 0.000000 640.000000 0.000000 480.000000

Then I resize the display window manually and rerun;

info = grass.read_command("d.info", flags='r')

The output remains the same, despite me having manually resized the display window. The same problem occurs when I run d.info with any of its flags though (I'm not entirely sure what each flag does).


texts: Correct installation of postgresql and postgis

texts: Goodmorning everybody! I'm new in this forum and I hope you would help me in solving my problem. I'm practicing the use of Qgis and i need to setup postgresql and postgis on my mac yosemite 10.10. Does anybody of you have a correct step-to-step guide which shows how to setup correctly that kind of local server?


texts: QGIS using project variables in filter queries

texts:
In QGIS 2.14.11 LTR I have to apply the same (or similiar) filter query to various layers. I have the road name to filter for in a project variable and try to apply the following query:

"road_name" = @current_road

This works as expected in a 'select by expression', but returns 0 row as a filter query, where I obviously have to use

"road_name" = 'the road name'

As variables are obviously not evaluated in filter queries, is there a deeper reason for this that I don't see, and is there another way/syntax to use variables in filter queries?


texts: Using ArcPy Insert Cursor with NEXTVAL of Oracle sequence?

texts:
I have the need to be able to set a value of a gdb table to the NEXTVAL of an Oracle sequence. I have to be able to do this through an ArcPy function that will also handle setting the ArcPy OBJECTID in the table (such as InsertCursor will automatically increment that field).

I have tried the following with no luck:

my_insert = arcpy.InsertCursor(My_SDE_Table_Path)
new_record = my_insert.newRow()

new_record.setValue('SEQUENCED_ID', 'MYSEQUENCE.NEXTVAL')

my_insert.insertRow(new_record)

This just attempts to write the string "MYSEQUENCE.NEXTVAL" to that field.

Anybody know of a solution to use sequences with ArcPy? I've scoured the Esri documentation but it's like a labyrinth.


texts: How to find layer ID and name from Identify query using ArcGIS Web ADF?

texts:
I need to retrieve the layer ID of every layer returned by the IQueryFunctionality.Identify function (I am using ADF with .NET)

Identify function spits out a DataTable array for layers which have answered to the query. Since there's no way to know what layers have not answered, you need to check the TableName (which is the layer name) for every DataTable in the array.

But what if you have two layers with the same name? You can't know what layer answered and what layer didn't

The closest answer I could find on the internet is a link at esri forums but the topic got disbanded before reaching a solution (which may not exist perhaps).

This sounds like a problem to me, can someone please help me?


texts: Reducing number of polygons in a shapefile to increase the render speed for R Leaflet?

texts:
I'm trying to increase the speed of my R Leaflet map by reducing the size of the polygons.

I found a great solution on stack exchange here, but for some reason, I get the following error when I try to plot the shapes. "Error in pgons@Polygons[[index]] : subscript out of bounds"

I'm guessing the issue stems from the getSmallPolys function, but I don't know how to fix it. If this would work, it would increase the speed of my Leaflet map a lot.

library(sp)
library(rgeos)
library(raster)
library(leaflet)
library(readr)
# Source of getSmallPolys function: "https://gis.stackexchange.com/questions/62292/how-to-speed-up-the-plotting-of-polygons-in-r"

# Get the main polygons, will determine by area.
getSmallPolys <- function(poly, minarea=0.01) {
  # Get the areas
  areas <- lapply(poly@polygons,
                  function(x) sapply(x@Polygons, function(y) y@area))

  # Quick summary of the areas
  print(quantile(unlist(areas)))

  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))

  # Get only the big polygons and extract them
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]] >= 1){
      poly@polygons[[i]]@Polygons <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:length(poly@polygons[[i]]@Polygons)
    }
  }
  return(poly)
}

AUS <- getData("GADM", country = "AUS", level = 1)
CAN <- getData("GADM", country = "CAN", level = 1)
USA <- getData("GADM", country = "USA", level = 1)
GBR <- getData("GADM", country = "GBR", level = 1)
shape <-
  rbind(AUS, CAN, USA, GBR)

# From 195 mb to 78.7 mb (32.1 mb if filtered)
shapesSimple <-
  gSimplify(shape, tol=0.1, topologyPreserve=TRUE) %>%
  SpatialPolygonsDataFrame(., data = shape@data)

# This would be ideal for render speed but I get an error when I try to plot it
# From 195 mb to 1.9 mb (1 mb if filtered)
shapesSimple <-
  gSimplify(shape, tol=0.1, topologyPreserve=TRUE) %>%
  getSmallPolys() %>%
  SpatialPolygonsDataFrame(., data = shape@data)

# Testing the speed
start_time <- Sys.time()
leaflet() %>%
  addTiles() %>%
  addPolygons(data = shapesSimple)
end_time <- Sys.time()
end_time - start_time


texts: Python script to perform the same operation over multiple different point shapefiles

texts:
From classes at uni, I have just some basic knowledge of Python. After an introduction to variables, datatypes, lists etc., we were using GDAL/OGR and despite some troubles, I was able to write some simple scripts performing geo-operations. We were also only very briefly introduced to ArcPy.

Now, I am trying to process tasks for my bachelor thesis where I have 11 shapefiles of geocoded data. I need to perform the same operation over each of the point shapefiles (they contain very similar data - addresses of students' studying at our faculty during the given academical year). I am sure it must be possible to write a simple script that would iterate through each of the shapefiles and perform desired operation but I am unable to put it together.

Could anyone at least give me a hint how to do it?

I am aware that this is not a place where I ask a question and you do everything instead of me.

I suppose I should be using ArcPy to do this, am I right?

What I specifically want to do is extract (like select and put them in a new shapefile) those points which are inside buffers I created. I have 5 concentric rings with given perimeters.

You probably understand that clicking 5×11 times very similar operations is not what is particularly scholarly.

Which part of 'help' should I explore?


texts:
I have an issue on the available fonts on the installed GeoServer.

is a GeoServer with available fonts correctly installed.

is the GeoServer with the error.
How can I correct this issue?


texts: Clip/Intersect tool is failing to process all data in my shapefile

texts:
I am using ArcMap 10.3.
I have this polygon which contains information about different types of soil. I want to clip it, so that only the area inside the other polygon (black border) will be in the output.

However, when I use the clip tool, this is the output:
There is large areas that have not been implementet in the output file, even though it is inside the clip polygon. There is also an area in the lower right corner that is outside of the clip polygon that is in the output file.

And if I use the intersect, this is the output (the clip output is over the intersect layer, showing the difference):
The intersect tool managed to get some more areas in the output, but also here, some areas outside of the polygon are inside the output file.

I used only standard options when performing these tasks, but it seems like ArcMap has troubles with some of the polygons.
Any ideas of what I am doing wrong?


texts: I recently reformatted my hard drive, so I downloaded QGIS so I could install it again. But now that I have, I seem to be missing all of my geoprocessing tools. Normally there were quite a few available from the drop-down Vector menu. Now there are only a handful. Also, I can't find the processing toolbox that I used to use, not sure what that was called but it seems to be missing now. It was where you could access all of the geoprocessing commands including GRASS.


texts: Setting raster symbology to standard deviations in ArcPy

texts: I have an ArcPy routine and I would like to set the symbology of a derived raster to standard deviations, so that I can then reclassify on that basis and convert to polygons. Is this possible? I looked at the documentation and it appears the options in symbology are limited.


texts: Atlas - Export one image for each street of multiple areas

texts:
I have 2 layers:

OSM roads polyline layer
various area polygon layer

I need to export each area with every road highlighted.

Here are 2 exemples of export of the same area with 2 differents roads:

I selected the road manually and saved the image manually.

The problem is that there are about 10 areas with more than 100 roads.

Is it possible to automate the image export for each areas and roads ?

With the atlas tool I can export automatically each areas with this window:

Exemples of 2 areas:

How can I highlight every road one by one for each area and export it as image automatically ?


texts:
So I'm writing some script for QGIS 1.8 and I'm unable to resolve this error.

from qgis.core import *
pallyr = QgsPalLayerSettings()

The second line gives me error that QgsPalLayerSettings is not defined. When I do it in QGIS 2 or higher then i don't get any errors and it imports this class. Somehow I have to stick to QGIS 1.8. How can this be resolved?


texts: Can't start web services (GeoServer) in OSGeo-live 11 Lubuntu

texts:
I have created a virtual machine using Osgeo-live 11 in Lubuntu and when I try to start the GeoServer or any other web services they start correctly but when I run it on web I get errors like in the picture below:

I have used this previous package of OSGeo live like 6, 8, 9 and everything was OK, in this release I don't why the web services don't work.


texts:
I'm getting a Lenovo thinkpad w520 with 8gb of ram/2000m quadro GPU/2760qm Corei7 to use as my secondary work laptop.

I mostly deal with large amounts of ortho photos and some GPS software.

Do you guys think this is overkill or any of these specs useless? I'm just hoping this was a good purchase, I didn't think a mid consumer laptop would be able to run GIS very well (without crashing)....given that I do not do much 3D work, just alot of different GIS tasks, is the Quadro 2000m going to benefit me more than say a Nvidia 540m mid level card?


texts:
I use a map control and set a MapClick event for it.

And i have a tool that ClientAction set to point.

Event fired when i clicked on the map.

But args.ScreenPoint was empty.

What's wrong ?

I want get x and y in map click.


texts: Close dialog when algorithm is finished (QGIS plugin)

texts:
I've read many threads on this but I'm not able to figure out the solution (I'm a beginner to the plugin creation...).

What I what to do is to close the dialog window when the algorithm has finished its job.

That's the piece of code of the alg.dialog.py:

    class CreateModelDialog(QDialog, FORM_CLASS):
        def __init__(self, iface):
            QDialog.__init__(self)
            self.iface = iface
            self.setupUi(self)
            self.buttonBox.rejected.connect(self.reject)
            self.buttonBox.button(QDialogButtonBox.Ok).clicked.connect(self.createModel)

            self.manageGui()

........

    def createModel(self):
        modelName = self.bxDBname.text()
        modelType = self.cmbBxModelType.currentText()
        lengthString = self.cmbBxLengthUnit.currentText()
        timeString = self.cmbBxTimeUnit.currentText()
        workingDir = self.OutFilePath
        print 'Direttori recuperata ..... ' , workingDir
        isChild = 1.0
        createModel(modelName, workingDir, modelType, isChild, lengthString, timeString)

So self.buttonBox.button(QDialogButtonBox.Ok).clicked.connect(self.createModel), when clicked runs the createModel method.

But the dialog remains opended, where should I edit the code?

This is the __init__ function:

def __init__(self, iface):
    QDialog.__init__(self)
    self.iface = iface
    self.setupUi(self)
    self.buttonBox.rejected.connect(self.reject)
    self.buttonBox.button(QDialogButtonBox.Ok).clicked.connect(self.createModel)

    self.manageGui()

This is the function in the main plugin.py file:

def runCreateModel(self):
    """Run method that performs all the real work"""
    # Create the dialog (after translation) and keep reference
    self.dlgMD = mdDialog.CreateModelDialog(self.iface)
    # show the dialog
    self.dlgMD.show()
    # Run the dialog event loop
    self.dlgMD.exec_()


texts: Can Geoserver return the raster value of a lat/lon point

texts:
I have a GeoTIFF published by Geoserver 2.2 and being used by Leaflet via WMS. How can I get Geoserver to return the raster value at a lat/long point?

I'm guessing it has to do with using WCS, but no examples seem to exist for doing this!

This will be similar to querying PostGIS with ST_Value(raster, ST_SetSRID(ST_Point(lon,lat),4326)). I dont really want to import the same raster into PostGIS just to find the raster value at a point. Or is this the recommended solution?


texts: Geoserver CSS Styling - Bug using env var twice?

texts:
Currently working on a css style for one of my layers. My plan is to use an env var from the query string to drive some of the styling. I was going to use the env var as a flag string. Basically, if 'A' is in the string, that will trigger labels to show, if 'B' is in the string, icons will show, etc, but I'm running into an issue where it seems like the env var is causing my conditionals to evaluate incorrectly.

My style

/* @title label */
[strIndexOf(@style_flags, 'A') > -1] {
    label:          [strIndexOf(@style_flags, 'A')];
    font-family:    "Trebuchet MS Bold";
    font-fill:      white;
    font-style:     normal;
    font-size:      12;
    halo-radius:    0.5;
    halo-color:     black;
}

I open the layer in geoserver's OL previewer, and append '&env=style_flags:A' to the end.

But when I change style_flags to 'B', I get this.

If strIndexOf(@style_flags, 'A') evaluates to -1, the rule should evaluate to false, and no label should show, right?

EDIT:

I was appending this '&env=flags:ABC' to the OL URL and using subsets of 'ABC' to test the toggle functionality. I played around with setting defaults on the env var, because it appeared that my flags var string wasn't overriding the default. I used labels to test if my variable was being read in, and it was, but it didn't seem to affect how the rules were evaluated.

my full style:

@mode "Flat";

/* @title 1-1000 */
[strIndexOf(@flags("C"), 'C') <> -1] [speed BETWEEN 1 AND 1000] {
    mark:       url('http://localhost:8080/geoserver/www/styles/images/arrow2.png');
    mark-size:  [speed*0.25];
    mark-rotation: [course];
}

/* @title Over 1000 */
[strIndexOf(@flags("C"), 'C') <> -1] [speed > 1000] {
    mark: url('http://localhost:8080/geoserver/www/styles/images/arrow2.png');
    mark-size:  125;
    mark-rotation: [course];
}

/* @title PLANE BLUE */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'PLANE'] [color like 'BLUE']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/PLANE_BLUE.png');
}

/* @title PLANE RED */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'PLANE'] [color = 'RED']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/PLANE_RED.png');
}

/* @title PLANE GREEN */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'PLANE'] [color = 'GREEN']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/PLANE_GREEN.png');
}

/* @title CAR BLUE */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'CAR'] [color like 'BLUE']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/CAR_BLUE.png');
}

/* @title CAR RED */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'CAR'] [color = 'RED']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/CAR_RED.png');
}

/* @title CAR GREEN */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'CAR'] [color = 'GREEN']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/CAR_GREEN.png');
}

/* @title SHIP BLUE */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'SHIP'] [color like 'BLUE']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/SHIP_BLUE.png');
}

/* @title SHIP RED */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'SHIP'] [color = 'RED']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/SHIP_RED.png');
}

/* @title SHIP GREEN */
[strIndexOf(@flags("B"), 'B') <> -1] [type = 'SHIP'] [color = 'GREEN']{
    mark: url('http://localhost:8080/geoserver/www/styles/images/SHIP_GREEN.png');
}

/* @title label */
[strIndexOf(@flags("A"), 'A') <> -1] {
    label:          [name];
    font-family:    "Trebuchet MS Bold";
    font-fill:      white;
    font-style:     normal;
    font-size:      12;
    halo-radius:    0.5;
    halo-color:     black;
}


texts: Georeferencing an image to a Progressive JPEG

texts:
I am using the Georeferencing tool in ArcGIS 10.2 to georeference a normal JPEG image.
I want to produce a georeferenced Progressive JPEG image from that image, but there is no JPEG Progressive option in the Rectify window dialog:

It seems ArcGIS has the capabilities to produce Progressive JPEG, just not through the georeferencing toolbar. Rather, it can produce a Progressive JPEG by exporting the map to an image through File - Export Map:

Exporting my georeferenced image to a Progressive JPEG through Export Map would be cumbersome and maybe cause loss of image quality.

Do you know any workaround to this?


texts: How to maintain data in attribute table when importing ASCII file to ArcMap

texts:
I am very new to using ArcGIS. I have an Excel sheet with data in it where each cell represents a grid from a model output and contains a value (eg. amount of precipitation in that grid cell from the model output). I have taken the data and put it into ASCII format in order to import it to ArcMap. The data appears to import and display correctly on the map, but the attribute table does not maintain the original model out put values (eg. precipitation and location). Is there a way to change the attribute table to better reflect the data and get it in a format similar to what I started with? I will use these numbers later on for further analysis.

Below is a screen shot containing the map view, attribute table, and original text file that I imported.


texts: r.Viewshed GRASS function returns NaN insteat of 0 and 1 for 1 in -b Boolean Mode

texts:
r.Viewshed GRASS function returns NaN insteat of 0 and 1 for 1,
i am not not able to recode NaN Values to 0 with r.recode function.
How to reclassify or recode raster from 0 to 1? pl. Help


texts:
Is there a generic issue with the QgsPluginLayer class in QGIS 2.0.1? When I take the basic Watermark sample from here:
https://github.com/sourcepole/qgis-watermark-plugin

and make a couple of mods so that the plugin runs on 2.0.1 (all that appeared to be required was to add a metadata.txt file) it fails with an error "Could not draw Watermark plugin layer because:" and thats it. It seems like there is a similar issue active with an OpenLayers plugin as reported here:
http://hub.qgis.org/issues/8679

but I am working on a plugin using the QgsPluginLayer class which is also now broken so I am wondering whether anyone can suggest a workaround or have any idea about when a fix may be available?

Thanks in advance!


texts:
I converted my shapefile to autocad.
How can I see my fields in autocad.


texts: Error installing Shapely package using pip command

texts:
I'm using Windows 10 and python 2.7. Even though I updated pip it's getting error:

ERROR: Command errored out with exit status 1:
     command: 'c:\python27\python.exe' -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'c:\\users\\pc\\appdata\\local\\temp\\pip-install-eztkma\\shapely\\setup.py'"'"';
__file__='"'"'c:\\users\\pc\\appdata\\local\\temp\\pip-install-eztkma\\shapely\\setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' egg_info --egg-base pip-egg-info
             cwd: c:\users\pc\appdata\local\temp\pip-install-eztkma\shapely\
        Complete output (9 lines):
        Traceback (most recent call last):
          File "<string>", line 1, in <module>
          File "c:\users\pc\appdata\local\temp\pip-install-eztkma\shapely\setup.py", line 80, in <module>
            from shapely._buildcfg import geos_version_string, geos_version, \
          File "shapely\_buildcfg.py", line 200, in <module>
            lgeos = CDLL("geos_c.dll")
          File "c:\python27\lib\ctypes\__init__.py", line 362, in __init__
            self._handle = _dlopen(self._name, mode)
    WindowsError: [Error 126] The specified module could not be found
    ----------------------------------------
ERROR: Command errored out with exit status 1: python setup.py egg_info Check the logs for full command output.


texts:
I would like to override the default contextual menu in QGIS 3 in the ToC.
I've created a class that inherit from QgsLayerTreeViewMenuProvider.
from qgis.gui import QgsLayerTreeViewMenuProvider
from qgis.PyQt.QtWidgets import QMenu

class LayerTreeViewMenuProvider(QgsLayerTreeViewMenuProvider):

    def __init__(self):
        super().__init__()

    def createContextMenu(self):
        menu = QMenu()
        menu.addAction('Test', lambda : print('test'))
        return menu

This code works well but I would like to get in all cases (right click in ToC or to a vector layer, raster layer, ...) the default contextual menu to adapt it. So I've made :
def createContextMenu(self):
    menu = super().createContextMenu()
    return menu

That doesn't work, I've got a RecursionError: maximum recursion depth exceeded


texts: Why memory map file is larger than original slide

texts:
I am having a tiff slide with the dimension of 90112 x 100352. The slide size is 1.3GB . To process this slide i am creating a numpy memmap file with above dimension like below.

If I check the size of the "test.memmap", its showing nearly 27GB. I dont know why its this much bigger than original tiff file.

I am reading the original slide patch by patch and copying its RGB values with "memMapRef". After reading whole slide i am drawing a rectangle on "memMapRef" at some X,Y coordinates. After drawing the rectangle, now I am trying to create image using below code from "memMapRef" numpy array. Here i am getting the Error "size does not fit in an int".

I am using the "imsave" method from scipy.misc.

My Entire code is below:

import cv2
import numpy as np
from scipy.misc import imread, imresize, imsave
import openslide
level = 0
read_size = (10000, 10000) # reading size per patch
tumor_file = "Test_Tumor.tif" # 1.3 GB file
tumor_slide = openslide.open_slide(tumor_file)
max_OX,max_OY = tumor_slide.level_dimensions[level] # max_OX = 90112 max_OY = 100352
memMapRef = np.memmap('test.memmap', dtype='uint8', mode='w+', shape=(max_OY,max_OX,3))
def process_slide():
    for start_OX in xrange(0,max_OX,read_size[0]):
        for start_OY in xrange(0,max_OY,read_size[1]):
            size_OX=min(read_size[0],max_OX-start_OX)
            size_OY=min(read_size[1],max_OY-start_OY)
            tile = np.array(tumor_slide.read_region((start_OX,start_OY),level,(size_OX,size_OY)), dtype="uint8")[:,:,:3]
            yield tile, start_OX, start_OY

for img_tile, X, Y in process_slide():
    h,w,d = img_tile.shape
    memMapRef[Y: Y+h, X: X+w] = img_tile

x1 = 500
y1 = 1000
x2 = x1+256
y2 = y1+256
cv2.rectangle(memMapRef,(x1,y1),(x2,y2),(0,255,0),10)
imsave("output.png", memMapRef) # Error "size does not fit in an int".

Where I did mistake and How to fix it?