texts: Raster pixel aspect ratio gets distorted when exporting from Google Earth Engine

texts:
Problem:
I have a mosaic of a Landsat 7 collection in Google Earth Engine that looks alright on display (pre-export-version) in the map window of the web app. However, when I export it, some pixels get narrower and some wider. Basically there are north-south columns with pixels that are 22 m wide instead of 30m. (see attached screenshot)
Question:
Why does this happen? Is there any way to avoid it?

Investigation:
I investigated the output in QGIS, but if I upload the exported image back into Earth Engine this difference remains, so it does not appear to be a platform dependent problem.
If I use the inspector function, the pre-export and the export-imported file have different pixel values where the overlap is distorted, so it is not just a display error.
If I do the same thing with an image from the same collection (just the raw image, no mosaic), the exported file looks identical to the original.
I have tried to export using different coordinate references such as WGS84 and different UTM-zones. None of this affects the result.
Thoughts:
First I thought that difference between the mosaic and the raw image could imply that the irregularities are a result from the mosaic operation trying to fit pixels with slightly different spatial location into the same raster. But if so, I do not really get why that is not showing in the pre-export render of the raster.
Reproducible example: (or see link: https://code.earthengine.google.com/3ecceca28e81cc3d0c304281f08b4619)
// Define study site geometry:
var box = ee.Geometry.Rectangle(14.91081234991633, 68.28371077192794, 14.965572329164377, 68.30250393637535);
Map.centerObject(box);

// Filter image collection by time and area.
var landsat_7 = ee.ImageCollection("LANDSAT/LE07/C01/T1_SR") //
  .filterBounds(box)
  .filterDate('2000-07-01','2000-08-31');

// Make median mosaic from collection:
var landsat_7_mosaic = landsat_7.median()
  .select('B.+'); // Keep only bands containing B

var visuals = {bands: ["B1","B2","B3"],
    gamma: 1,
    max: 2737.5867438745377,
    min: 2040.6889200176545,
    opacity: 1};

// Add the image to the map window
Map.addLayer(landsat_7_mosaic, visuals, 'Landsat 7 mosaic');

Export.image.toDrive({image: landsat_7_mosaic,
                       description:'imageToDriveExample',
                       scale:30,
                       region: box,
                       maxPixels:1e13});

// Import the export and add to map to compare with the pre-export


texts: Issue with webADF 931 application migrated to 10 environment within IE8

texts:
I am currently in process of migrating a webADF 931 application to Arc 10 environment and testing it within IE8. I am having this issue, the map when launched only covers partial area. I have tried setting various sizes but didn’t made any difference. I have atatched the screenshot here. if you see the area on the right and bottom these aren't filled.

Thanks and look forward to any feedback.
Jay


texts: Dynamic copyright or layer credit in QGIS 3.4?

texts:
Is it possible to add a layer credit to the bottom right corner of the map and print composer in QGIS 3.4? This fuctionality existed in 2.8, as seen here.

I can't find this option in the layer properties anymore. I realize I can do something similar manually via the Copyright Decoration, but it was nice to just check/uncheck a box when needed. The Copyright Decoration is also always on top and visible.

Perhaps someone can suggest an expression to use in the Copyright Decoration that will pull the visible layer's metadata and use that as the text string? My QGIS expression skills aren't that strong yet.


texts: Running more than 10 routes at same time in ArcGIS Network Analyst?

texts:
My team and I are trying to run food collection routes for a food bank in Mexico City but we can run more than 10 routes at a time. In order to run the eleventh we need to take one out.

Does someone know what could be the problem or where to configure this?


texts: QGIS 3.8 Batch Processing AttributeError: 'BatchInputSelectionPanel' object has no attribute 'value'

texts:
I installed the latest version of Qgis (3.8.0). When I try to run batch processing for contours, it gives me this error

Is anyone else experiencing this? I cant seem to find a quick solution to this problem.


texts: Change style of GeoJSON circle marker by feature properties

texts:
I have GeoJSON Feature Collection of points with a text property 'Proyecto' and I want to change the color of the circle marker for each point to vary by that property.

I've created a function to get the color based on the five unique values for the 'Proyecto' property:

function getColor(d) {
  return d = 'Vivienda saludable' ? '#a6cee3' :
         d = 'Programa OPAS- 1969, Prevención de conflictos, desarrollo de acuerdos y construcción de la paz en comunidades con personas internamente desplazadas en Chiapas, México'  ? '#1f78b4' :
         d = 'Familias fuertes, amor y límites'  ? '#b2df8a' :
         d = 'e-Health, acceso a servicios médicos de telemedicina en comunidades indígenas en extrema pobreza'  ? '#33a02c' :
         d = 'Casas maternas'   ? '#fb9a99' :
             '#FFEDA0';
         }

I then create the GeoJSON layer, call pointToLayer, create a circle marker and call the getColor function, and create a popup from the properties of the GeoJSON file.

var puntos = L.geoJson(puntos_chiapas, {

    pointToLayer: function(feature, latlng) {
        return new L.CircleMarker(latlng, {radius: 10, fillOpacity: 0.85, color: getColor(feature.properties.Proyecto)});
    },
        onEachFeature: function (feature, layer) {
            layer.bindPopup(
                "Proyecto: " +
                feature.properties.Proyecto +
                "
" +
                "Ubicación: " +
                feature.properties.Ubicación +
                "
" +
                "Sector: " +
                feature.properties.Sector +
                "
" +
                "Esquema: " +
                feature.properties.Esquema
            )
            }

    });

The resulting map and points all show only the color of the first 'Proyecto' text from the getColor function. What am I doing wrong?

I've included the JSFiddle here. You may need to toggle off the polygon layer to see the points.


texts: gdal.RasterizeLayer can not burn any values or use attribute from shapefile in python

texts:
I am trying to rasterize a shapefile but the raster produced by the code it has only 0 as a value. The extends and projection is correct but I am not able to set any value in the raster. I need to use the "lucode" from the shapefile. This is an integer field in the shp. I have aso tried to burn a single value to the raster with no luck.
below please find the code I am working on:
from osgeo import gdal, ogr
# Define pixel_size and NoData value of new raster
pixel_size = 5000
NoData_value = -9999

# Filename of input OGR file
vector_fn = r'D:\Africa\01_carbon\DRC\shp\LULC_join.shp'

# Filename of the raster Tiff that will be created
raster_fn =  r'D:\Africa\01_carbon\DRC\shp\test.tif'

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
source_srs = source_layer.GetSpatialRef()

# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)

target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

target_ds.SetProjection(source_srs.ExportToWkt())

#print(target_ds, [1], source_layer.name)

# Rasterize
#gdal.RasterizeLayer(target_ds,
#                    [1], source_layer,
#                    burn_values=[3]) # try to burn a contant value (3)

# Rasterize
gdal.RasterizeLayer(target_ds,
                    [1], source_layer,
                    options = ["ALL_TOUCHED=TRUE", "ATTRIBUTE=lucode"]) #using an attribute from the shapefile


texts: (d.*) of GRASS GIS does not execute from standalone script

texts:
I am in windows 10. All my module import works well and grass env for standalone script is fine. I tried to execute following commands from standalone script. Neither wx1 display gets executed nor it throws any error.

gscript.run_command('d.mon',start='wx1',resoultion='1')
gscript.run_command('d.rast',map='elevation')


texts: Does anyone know how to import Maki Icons to QGIS?


texts: GeoServer layer disappearing from Leaflet map at max zoom?

texts:
I am developing an application that uses GeoServer + Leaflet + PostGIS. I created my database-based layers within GeoServer. As a basemap I am using OSM. When the maximum zoom occurs with the active layer, the style is losing completely, but remains under the map.

I also use the "Leaflet.Wms" plugin to access the GeoServer layers. Follows the code:

var NewSource = L.WMS.Source.extend({
'showFeatureInfo': function (latlng, info) {
    this._map.openPopup(info, latlng);
}});
var wmsCorredores,
    wmsLotes,
    wmsLimites,
    wmsMunicVizinhos,
    wmsZonas;

var s = new NewSource(url, {
    "format": 'image/png',
    "transparent": true,
    "version": '1.1.1',
    "attribution": '&copy; PMSJC',
    "info_format": "text/html",
    "feature_count": 10
});

wmsMunicVizinhos = s.getLayer("geoSJC:vizinhos").addTo(map);
wmsLimites = s.getLayer("geoSJC:limites").addTo(map);
wmsPMSJC = s.getLayer("geoSJC:grupo_geoSJC").addTo(map);

After some tests I could conclude that when the zoom reaches level 19 the layer style is lost. The layer is still active, only the style is lost. The style applied to the layer is sent by GeoServer via WMS service.

I can't understand if the problem is setting the style inside GeoServer, or if I need to fix the application that consumes these layers. As I am still getting to know GeoServer, I would like to know if anyone has had this kind of problem.


texts:
I have implemented K-means and Expectation- maximization algorithm on my sentinel-1 SAR image in the SNAP toolbox.
Now I want to compare which clustering result.is better using beta -cv or silhouette coefficient or any other coefficient to test compactness of the clusters
If anyone could help me with it.


texts: PyQGIS & Qt - Revisiting the 'dialog window on top' issue

texts:
I want to keep a modeless dialog window on top of the main QGIS window. This issue has been discussed before and two solutions have been offered:

self.dlg.setWindowFlags(Qt.WindowStaysOnTopHint). This is
unsatisfactory because the dialog window stays on top of every application
window, not just the QGIS one.
Give the dialog an owner when it is created, as in self.dlg = MyCreatorDialog(self.iface.mainWindow()). This does indeed keep the
dialog on top of the application's main window, but sets the dialog as modal, so the underlying
window is inaccessible.

The processing toolbox creates modeless dialog boxes that stay on top of the application window, but not on top of other applications. How does it do it?


texts:
I am digitizing road line in India using the OpenLayers plugin. When I calculate the length of line using the field calculator it gives me the result in millimetre.

I thought it is the projection problem and defined the projection as WGS 84 UTM 44N, but the same problem exists. When I measure the line using measure line tool it gives the result in meters.

What is the problem?

I want the result to be in meters in the field calculator itself.


texts: Calculating percent cover from shapefiles

texts:
I normally do not work with shapefiles, so I am a bit lost here. I have two shapefiles each with multiple objects. The first is a set of 32 polygons (each one is a plot). The second shapefile has >10,000 objects which represent vegetation clusters of different sizes within each plot. I am trying to figure out:

1) How do I calculate percent cover of total vegetation cover within each site?

2) What percentage of each the vegetation cover is less than 5 meters in area in each plot?

This is what my data looks like in ArcGIS for a single plot. I am open to advice on how to do this in arcGIS too.

This link will take you to the shapefiles in my dropbox:
https://www.dropbox.com/sh/wyokxximppexyb3/p7VC-pfF2E


texts: Basemap not appearing: upgrade from OpenLayers 3.0 to OpenLayers 3.3

texts:
I recently updated my work's app to OL 3.3 from OL 3.0. I had a few basemap layers that looked like this:

    "basemap1": new ol.layer.Tile({
                extent: new ol.proj.transformExtent([-180.0,-90.0,180.0,90.0], "EPSG:4326", "EPSG:3857"),
                source: new ol.source.XYZ({
                        tileUrlFunction: function(coordinate) {
                        if(coordinate == null){
                                return "";
                        }
                        var z = coordinate[0];
                        var x = coordinate[1];
                        var y = (1 << z) - coordinate[2] - 1;
                        return 'http://tileurl.com/'+z+'/'+x+'/'+y+'.png';
                        }
                }),
                visible: true
    })

Which worked just fine in OpenLayers 3.0, but for some reason aren't appearing in OL 3.3.


texts: Export image collection animation generated by Google Earth Engine

texts:
I have created an animated thumbnail of an image collection using the following piece of code:

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

How can export the thumbnail to Google Drive as a GIF or other similar format?


texts: Creating grid layer of 2 km sq. and deciding X and Y parameters?

texts: How to create a Vector Gird polygon layer and what should put in the X and Y parameters when I would like to create a 2 km sq grid polygon


texts: QGIS Python bindings to load Cloud Optimized GeoTIFF (COG)

texts:
I am attempting to load a Cloud Optimized GeoTIFF (COG) via Python in QGIS 3.4.7. Specifically, I can manually load a COG using the tutorial posted here.
The URL for the COG is: http://oin-hotosm.s3.amazonaws.com/59d33df023c8440011d7b26d/0/b378087a-c2a5-43a0-abec-71fcfb051150.tif
I have tried using the AddRasterLayer approach:
path = "http://oin-hotosm.s3.amazonaws.com/59d33df023c8440011d7b26d/0/b378087a-c2a5-43a0-abec-71fcfb051150.tif"
self.iface.addRasterLayer(path, 'name')

and as a URI:
uri = QgsDataSourceUri()
uri.setParam('url', path)
uri.setParam("identifier", path)
rlayer = QgsRasterLayer(str(uri.encodedUri()), 'my cog layer')
QgsProject.instance().addMapLayer(rlayer)

As far as I understand, GDAL should be loading these under the hood. In the case of the former, I am seeing:
2020-08-26T22:35:45     CRITICAL    Invalid Layer : GDAL provider Cannot open GDAL dataset http://oin-hotosm.s3.amazonaws.com/59d33df023c8440011d7b26d/0/b378087a-c2a5-43a0-abec-71fcfb051150.tif:
             An error occurred while creating a virtual connection to the DAP server:Error while reading the URL: http://oin-hotosm.s3.amazonaws.com/59d33df023c8440011d7b26d/0/b378087a-c2a5-43a0-abec-71fcfb051150.tif.ver.
             The OPeNDAP server returned the following message:
             Not Found: The data source or server could not be found.
              Often this means that the OPeNDAP server is missing or needs attention;
              Please contact the server administrator.
             Raster layer Provider is not valid (provider: gdal, URI: http://oin-hotosm.s3.amazonaws.com/59d33df023c8440011d7b26d/0/b378087a-c2a5-43a0-abec-71fcfb051150.tif

In the case of the latter it fails silently.
How does one load a COG into QGIS using the Python bindings?


texts: Read layer symbology class color with ArcPy

texts:
I am using ArcPy within ArcGIS 10.1 and attempting to read symbology from a *.lyr file created by a raster dataset within ArcMap.  My goal is to be able read out the RGB (color) values assigned to the class breaks on the dataset.
It’s seems easy enough to get the lyr file into an object with:

    lyrFile = arcpy.mapping.Layer("C:/.../ReadLyrFileSymbology/conus_foresttype6.img.lyr")

However, after that I am unable to find the method or property to get symbology color (RGB) on a layer.  Any pointers would be greatly appreciated.


texts: Making work tablet integrate GNSS (GPS) on Windows 8 with QGIS?

texts:
I have tablet - Dell latitude 10, running Windows 8, with integrated Broadcom GNSS Receiver BCM47511 (GPS,Glonass...). Because it is integrated GPS and Windows 8, there is the "new" way of providing GPS data called Windows Location Provider. This "new" way means it doesn't work properly or at all in almost any navigation/localisation software at this moment - because this software expect GPS data coming from COM port. Including QGIS - at this moment it is on the road map - http://hub.qgis.org/issues/7878

Anyway after extensive search I found software which can bypass MS localisation platform. It is called Centrafuse Localizer. This software can simulate COM port and send through it GPS data so some software is able to read it.

Unfortunately not QGIS. QGIS can see the port (just to test it I changed COM 1 to 5 and 7) and it changes Serial device list in QGIS but after hitting Connect it always fails.

If you are thinking about tablet with integrated GPS think twice! From my point of view it is just too early - probably also no proper drivers yet as what I have read from Broadcom Readme file. Even when GPS "works" the refresh rate and precision is terrible comparing to my phone. And all software must be changed to work with this new interface...

Is there other way I could try to make work tablet integrated GNSS (GPS) on Windows 8 with QGIS?


texts: Trying to use SQL API from Oracle PL/SQL to CartoDB

texts:
I am trying to make an http call from within an Oracle server to CartoDB. To do so I have developed a custom PL/SQL package that will generate the SQL and make the calls. A sample
API call is below. The return status looks like a permissions problem ({"error":["permission denied for relation address"]}

Any ideas?

http://mycompany.cartodb.com/api/v2/sql?q=INSERT
INTO address(address_id,street_number,city,state,zip_code)  SELECT
'BBE20140321095101832AM5656' address_id, '3469 Shelby Ganges Road'
street_number, 'Shelby' city, 'Ohio' state, '44875' zip_code WHERE NOT EXISTS (
select 1 from ADDRESS WHERE ADDRESS_ID =
'BBE20140321095101832AM5656')&api_key=488e214c52f78b55f2fe9c05daef94074c7e3cd2


texts: How to copy or export list filter from one model to another?

texts:
I have a model tool with a variable parameter that is filtered by a string list (county names) I manually created this list filter when I created the model tool. I am now creating a new model tool and would like to use the same list of county names is there anyway to import the list filter from my other model so I do not have to manually create the list again. I am an amateur at best with python so if it is possible to keep the solution within the realm of ModelBuilder that would be best but if python is neccessary bear with me.

This is the value list I would like to reuse


texts: How to generate a raster showing min/max values from a stack of rasters? QGIS and Python

texts:
I have a stack of rasters (NDVI values over a year) and I want to generate a single raster showing the maximum difference between the layers.

So the result should be a single raster showing the difference between maximum value and minimum value for each point across the landscape. How can I do this in QGIS or Python?


texts:
I'm trying to label features using the geometry generator (QGIS 3.8). I used the following expression to build a simple circle :

buffer(centroid($geometry), 500)

In Placement i used Using Perimeter (curved)

The labeling is working pretty good except some of my features label are written clockwise and some counterclockwise. How could i control this rendering ?

Furthermore i would like the mid-point of my labeling string to be at the most northern point of the circle ...

How could i achieve that ?


texts:
I have this query :

SELECT ST_AREA(geom::geography) FROM region

The (geom) contains polygons from a shapefile imported into PostGIS with shp2pgsql.
I've appended ::geography, else the result of the function was displayed in degrees.
Now, I have this kind of answer :

544163926.22948

How can I know in which units it is displayed ? square/meters ? square/foot ?
I want it to be displayed in square/meters (m²) ou square/kilometers (km²)

How can I do that ?


texts:
In a custom layer I have a workspace variable (File GDB or SDE) that I need to persist (save/load).

Therefore I need to persist the ConnectionProperties of the workspace. Thankfully ConnectionProperties returns a PropertySet which implements IPersistStream so persisting this is trivial.
But I also need to persist the type of workspace factory to use.

Unfortunately WorkspaceFactory does not implement IPersistStream so persisting it directly does not work (returns null when reading).

Is there a way to do this?


texts: Getting nearest line to point with virtual layer using QGIS?

texts:
I have a streetmap and the position of hundred sensorlocations. I am currently trying to calculate the nearest line for each point. I have already done it with the nnjoin plugin. But now I want to do this with a virtual layer. For this, I followed the tutorial: https://pigrecoinfinito.wordpress.com/2017/04/28/qgis-minima-distanza-tra-punti-e-linee/

But unfortunately, I could not get it to work because I am very new to SpatiaLite.
To create the virtual layer, my approach looks like followed:

If I want to test this query, I am getting the following error note:

Any suggestions on how to fix this?

I have know tried it with renamed tables, but I am still getting the same error message.

Maybe something more is wrong?


texts: Calculating Pop density per sq mi within 1/4 buffer (transit route)?

texts:
I have census block data that I have clipped using my 1/4 mile buffer.
For each of my census polygons it has kept the original population, is there a way to get the population recalculated just within the buffer?


texts: Adding and visualizing unique values in raster dataset using ArcGIS Pro

texts:
I want to make a raster dataset out of points.
I have all the points already visualized on a map.
The attribute table (of the point dataset) consists of names of species (and other fields).
I want to visualize the number of different species in each raster cell.
For example: In one cell are 4 different species.
How can I count the different species and make a raster dataset that visualizes this?
The tool that I used is the point to raster.
However, I cannot select unique values and count them.
Counting only the number of points in each raster cell is not enough.


texts:
I am new to databases.

Is it possible to migrate from Oracle database to Postgres?

If yes, do you know any tested solutions?


texts: Is it possible to create Rule-Based Styling with $geometry intersecting a second layer?

texts:
I am trying to create a rule-based style based on geometries from another layer.

I have found an example that does exactly what I want, except that it uses a geometry from WKT:

  intersects( $geometry, geomFromWKT( 'POLYGON' ) )

(from here)

What I want to do is exactly the same thing, but with the second geometry based on a different layer rather than WKT.  As an example, I want to theme roads in black if they intersect a polygon layer, and grey if they do not. Is this possible?


texts:
I am trying to use the ol3 KeyboardPan, from the API documentation "this interaction is by default included in maps, the keys can only be used when browser focus is on the element to which the keyboard events are attached. By default, this is the map div".

I am not able to get this to work on my application or on any of the examples. How does one focus the browser on the map div? What am I missing?


texts: Creating raster that counts points using ArcGIS Desktop?

texts:
Is there a way to create a raster that counts points?

I have dresses (points) and I want to know how many addresses are in a square of 100m. So I also need a tool t measure one meter. I think it is quite simple, but I don't even now how this tool should be called. I'm using ArcGIS 10.4 on a Windows.


texts: Convert shp to kml in ArcMap but it is not placed correctly in Googleearth

texts:
I have a shp polygon in GGRS87/Greek Grid. I want to display the polygons on googleearth. I use "layer to kml" tool in ArcMap. However when I open the kml in googleearth the polygons are not aligned correctly (they are all moved southwest from their right position).

Has to do with the transformation of datums?

Notes:
1) in the tool I clicked on environments and as coordinate system I chose GCS_WGS_1984 and for transformation: GGRS_1987_To_WGS_1984
2) before applying the layer to kml tool, I tried to use Project tool to change the coordinate system to WGS84. But for an unknown reason I didn't have the option to use geogr. tranformation


texts: "Move feature" for points, manually moving a set of points with the mouse in QGIS

texts:
For moving a feature, the process is to edit the layer, use the "move feature" button, then drag and drop the feature.
I have a point shapefile layer, and would like to move a set of points. This video does exactly what I need to do at this timestamp: https://youtu.be/yiihopTnewg?t=50
However, in my case, there is no "Move feature" button. I can use the "Edit vertex" button, but I can only move one point at a time, and have to select in a different manner to the normal yellow selection process.
In the image below, I want to move the selected features to align with the square.


texts: ee.Algorithm.EEFlux Evapotranspiration (ET) Landsat Google Earth Engine

texts:
Does anyone know how to use ee.Algorithm.EEFux in Google Earth Engine to calculate Evapotranspiration with Landsat Images?

I am trying to use it but GEE does not recognise the algorithm.


texts:
I created a small (4x4) cost raster with a 10m cell resolution, and used the ArcGIS Cost Distance tool to calculate the Output Distance raster and the Output backlink raster, assuming a single source point in the bottom left (SW) corner. Then I used Dijkstra's algorithm to check the output distance raster values; they are all correct. However when checking the backlink raster, it appears to be incorrect. The cell value one east and 2 north in the backlink raster produced is 5, but should be 4 so it goes diagonally. This makes me wonder how often incorrect least cost paths are produced.  Am I missing something, or does ArcGIS often produce incorrect backlinks and therefore incorrect least cost paths?  I'm using ArcGIS for Desktop 10.3.

The cost raster (from Chang's textbook 7th edition) is:

  1 2 2 1
  1 4 5 1
  2 3 7 6
  1 3 4 4

The (correct) Output Distance raster produced by ArcGIS is:

  40.0  51.2  71.2  86.2
  30.0  55.0  84.9  92.4
  15.0  28.3  78.3  125.7
  0.0  20.0  55.0  95.0

The (incorrect) backlink raster produced by ArcGIS is:

  3 4 5 5
  3 5 4 6  <-- 5 should be a 4
  3 4 5 4
  0 5 5 5

You can see from the output distance raster that the correct direction would be 4: from 55 to 15 rather than to 30. Heck, even a 3 would go to 20 and be better than 30. Or is ArcGIS doing something I don't understand?


texts: How to create a new service area analysis?

texts: I need to create and run a new service area analysis. I'll need to do an analysis of how far you are coming by car from one hospital in 30 min, 60 min and 90 min. Is there anyone who can help me with a model of this in ModelBuilder?


texts: Why aren't 3D Globes used as the basis of GIS systems?

texts:
Maps have been around for thousands of years and have been widely accepted as the means of representing spatial information. We all know about several map projections,each with their merits and demerits.

However, I often wonder, with all the progress GIS has made, why not make good 3-d globes as the basis of GIS Systems. Most of the outputs of GIS are intended to be presented digitally and not printed. Having near accuarate representation of the curved earth will solve many problems regarding choosing the correct projection.

Is ArcGlobe and GoogleEarth (and the likes) the way forward to represent and analyse information upon? What do you think?


texts: How to reference the ArcGIS Server JS API using RequireJS

texts:
When using RequireJS to load a JavaScript library from a CDN, an example syntax is shown at https://github.com/requirejs/example-jquery-cdn:

requirejs.config({
    "paths": {
      "jquery": "//ajax.googleapis.com/ajax/libs/jquery/2.0.0/jquery.min"
    }
});

How should the ArcGIS Server JS API be referenced? The page at https://developers.arcgis.com/javascript/ gives the syntax for a standard web page:

<script src="http://js.arcgis.com/3.12/"></script>

but when using RequireJS:

requirejs.config({
    'paths': {
      'arcgis': '//js.arcgis.com/3.12/'
    }
});

this returns NetworkError: 404 Not Found - http://js.arcgis.com/3.12/.js

What is the correct syntax for referencing the ArcGIS Server JS API when using RequireJS?


texts: Using custom dash pattern for polygon outline (simple fill) in QGIS

texts:
I am trying to use a custom dash pattern for a polygon outline with a simple fill but I can't use anything other than the presets (see the screenshot)

I am using QGIS but looking for something similar to the "Cartographic Line Symbol" from ArcMap.

I know it is possible with a Simple Line, but is that achievable with the outline of a Simple Fill?


texts: Changing date on ESRI Basemap, to display different basemap

texts:
I am trying to create a map using the ESRI Basemap on ArcGIS 10.3.
The Problem I am having is that the area I want to map is covered by cloud, making details impossible to see. Is there a way to change the date on the ESRI Basemap, to a previous date so that I can eliminate the cloud cover?


texts: ogr2ogr / GDAL remove features in FileGDB while in use by ArcGIS Server map service

texts:
I'm using ogr2ogr to pull features out of Oracle and put them into a File GeoDatabase. The FGDB is being used by an active ArcGIS Server map service. This is all done for performance / reducing DB load. The DB data is refreshed nightly and then the ogr2ogr process runs so I don't worry about stale FGDB records.

The data can both grow and shrink overnight during the update so I need my FGDB refresh to both remove and add features.

I know this is possible if I use ogr2ogr's -overwrite option, e.g.

ogr2ogr --config FGDB_BULK_LOAD YES -f "FileGDB" <path to>.gdb OCI:user/pass@SID -overwrite -nlt LINESTRING -a_srs EPSG:4326 -nln layer_name -sql "select * from src_table"

But -overwrite fails if I don't first stop the ArcGIS Server map service.

I can use -update -append to update FGDB records while the service is running but this is 'additive' - records are added but nothing is removed.

ogr2ogr --config FGDB_BULK_LOAD YES -f "FileGDB" <path to>.gdb OCI:user/pass@SID -update -append -nlt POINT -a_srs EPSG:4326 -nln layer_name -sql "select * from src_table"

I'm looking for a way to either

Truncate the FGDB layer without removing it prior to an -update -append or
Change my ogr2ogr command to remove features not pulled from Oracle

I also saw this post this page on programmatically stopping / starting AGS map services but option 1 or 2 is my preference as I want an approach that doesn't require authentication / tokens over HTTP. If necessary I can bring in Python and perhaps the arcpy module.

Any input on GDAL / ogr2ogr parameters I haven't seen yet that might help would be much appreciated (or alternative approaches of course).


texts: Add a character every nth character in QGIS Field Calculator

texts:
I can't find an expression to add a character every nth character from the right of a string, in field calculator.

As an example, I have a list of UK postcodes which I need to add a space to 3 characters from the right:

From: N165KP To: N16 5KP

Note this has to be done from the right as the first lot of characters (before the space) can vary in length.


texts: Automatic generation of equidistant points inside a polygon in QGIS?

texts: I'm trying to find a way to generate a layer of points inside a given polygon of, for instance, 100 m x 600 m. Every points must be at least 20 meters away from the others and from the polygon borders, and the polygon must contain the maximum amount of points.


texts: Unsupervised classification using ISODATA clustering in R

texts: Is there a way (e.g. package) to perform an unsupervised classification using the ISODATA clustering algorithm in R?


texts:
I´m working currently through the QGIS training manual. However, the data for the spatial analysis exercise in chapter 9.4 (http://docs.qgis.org/2.2/en/docs/training_manual/complete_analysis/analysis_exercise.html) is not available under the link provided (can´t find the archive "exercise_data/more_analysis").

Somebody has the data at hand or can tell me where to download the archive.


texts: Drawing polygon using mouse click in Cesium?

texts:
I need to draw a polygon dynamically using my mouse. I can do this very easy in OpenLayers but I'm very new in Cesium and can't find any tutorials to do it.

Many examples only show how to create a polygon using hard coded geometries.

This is how to add a polygon using hard coded geometries:

var redPolygon = viewer.entities.add({
    name : 'Red polygon on surface',
    polygon : {
        hierarchy : Cesium.Cartesian3.fromDegreesArray([-115.0, 37.0,
                                                        -115.0, 32.0,
                                                        -107.0, 33.0,
                                                        -102.0, 31.0,
                                                        -102.0, 35.0]),
        material : Cesium.Color.RED
    }
});


texts: Joining attributes based on equal attributes and nearest neighbour using PyQGIS?

texts:
I'm relatively new to QGIS and Python.

Starting with two layers, one containing polylines (streets) and one containing points. I'm trying to automate some labeling of points. Therefore I want to add a field containing an angle of the street where the point is situated. This angle already exists in my line-layer (segments). Now I want to join the data on condition that the street-name is the same. So that I can automate the labeling (parallel to the street).


texts: Definition Query then Zoom to Layer using ArcPy?

texts:
I'm trying to create something in here. I want to have an automated process of doing a DEFINITION QUERY then ZOOM TO LAYER with one step.

the concept is like that, a dialog comes up with one form that asks me to put a unique id for the specific parcel that I'm trying to look up, then I type the parcel numder then hit OK then it brings me to that specific parcel directly.

Dialog pop-up
Empty form to be filled
Hit ok
It zooms directly to that parcel

How can I do this using ArcPy?

An automated way to do the definition query then zoom to layer.


texts: Problem with converting rasters to arrays and adding them to main array, then converting to raster via NumPyArrayToRaster : output is NoData values

texts:
I have a process that aims to loop through a list of rasters (testing with only 3)
covert each raster to an array, add the array to an existing 'template' array,  then save the final array back to a raster. I'm not getting errors, but all of the test output is NoData Values. Thanks for your help!

#import modules
import arcpy, numpy, os, time

#set workspace
arcpy.env.workspace = r"C:/GIS/project/test_rasters"

#set snap raster environment
template = "C:/GIS/project/rast01_0"
arcpy.env.snapRaster = template

#define in raster
template_raster = arcpy.Raster(template)

#set output location
outpath = r"C:/GIS/project/out_rasters/rast01_x"

#create an array to serve as the frame for additional arrays.
template_array = arcpy.RasterToNumPyArray(template_raster)
template_array = numpy.float64(template_array)
main_array = template_array * 0

#check array type
print main_array.dtype.name
#start counter
count = 0

#get a list of rasters in workspace
rasters = arcpy.ListRasters()
#for each raster in the workspace
for raster in rasters:
    print "processing: " + raster + " " + time.asctime()
    #make each raster an array, add it to main raster array
    array = arcpy.RasterToNumPyArray(raster)
    main_array = main_array + array
    #count number of raster/arrays added to the main array
    count = count + 1
    print count

#print final count of summed arrays
print count
#save array of sums to raster
part_sum_raster = arcpy.NumPyArrayToRaster(main_array)
part_sum_raster.save(outpath)

print"script completed:  " + time.asctime()


texts:
I am just wondering what the exact function of skewX and skewY in raster constructors is. According to the documentation linked:

  ST_MakeEmptyRaster — Returns an empty raster (having no bands) of
  given dimensions (width & height), upperleft X and Y, pixel size and
  rotation (scalex, scaley, skewx & skewy) and reference system (srid).

It seems that altogether (scalex, scaley, skewx & skewy) defines rotation. But obviously, scalex and scaley also defines the scale operation. This makes me wonder what exactly do these four parameters define.

Is the defined function/transformation equivalent to an affine transformation?

If it's just rotation and scaling, wouldn't it be more efficient to use scalex, scaley and the angle of rotation, e.g. theta?


texts: Highlight roads that are within a buffer or overlapping buffer

texts: So I have points with a mile buffer. Is there a way to "highlight" the roads that fall under the mile buffer? My original plan is to clip the roads with the buffer, but this will take some time because I have a lot of facilities. I was wondering if there's another tool in the Arctoolbox I can maybe use to speed the process. I can only use tools that come with ArcGIS and I'm on version 10.2. I guess I can also write a python script for my original plan.


texts: WMS MapServer OpenLayers 4 Output is blank

texts:
Tested the MapServer WMS Server with GetMap and GetCapabilities, but OpenLayers shows up blank. Using QGIS as WMS client also fails saying error downloading. The OSM layer in the OpenLayers works fine when the map projection is set to EPSG:3857, and doesn't work when EPSG:4326 is used. The local WMS service doesn't get displayed in either case. The OpenLayers template works fine. (But that's not even a WMS Service that's just the browse mode)

WMS Mapfile:

    MAP
    NAME "distt"
    IMAGETYPE PNG
    EXTENT 68 6 98 38
    SHAPEPATH "/ms4w/Apache/htdocs/data/"
    SIZE 600 600
    IMAGECOLOR 255 255 255
    PROJECTION
         "init=epsg:4326"
    END
    WEB
      IMAGEPATH "D:/ms4w/Apache/htdocs/data"
      IMAGEURL "D:/ms4w/Apache/htdocs/data"
      METADATA
        "wms_title"           "distt"
        "wms_onlineresource"  "http://127.0.0.1/ms4w/apache/cgi-in/mapserv.exe?map=D:\ms4w\apache\htdocs\data\MAP1.map"
        "wms_srs"             "EPSG:4326"
        "wms_format" "image/png"
        "wms_enable_request" "*"
        "ows_enable_request" "*"
      END
    END

    LAYER
      NAME "distt"
      TYPE POLYGON
      STATUS ON
      DATA India_districts_shp
      METADATA
      "wms_title" "distt"
      "wms_srs" "EPSG:4326"
      "wms_server_version" "1.1.1"
      "wms_enable_request" "*"
      "ows_enable_request" "*"
    END
    PROJECTION
         "init=epsg:4326"
    END
      CLASS
      NAME "Countries"
      STYLE
        COLOR 200 0 200
        OUTLINECOLOR 0 0 0
        END #STYLE
      END # CLASS
    END #LAYER

  END # MAP

HTML Openlayers:

    <!doctype html>

        <link rel="stylesheet" href="https://openlayers.org/en/v4.1.1/css/ol.css" type="text/css">

        <title>MNCFC CROP Map</title>

Map


texts: Calculating the error in a cut and fill between two surface models

texts:
I have a LiDAR suface model from 2006 and a LiDAR surface model from 2013.  These two surfaces have about a 0.4 foot (0.12m) difference between the Z values on an adjacent road (2013 is 0.4 feet higher).  I have a nearby agricultural field that has had soil removed from a 92,000 sq foot (0.86 hectares) area.

Using a drone, Metashape, and a Trimble 7X (to generate ground control) I was able to create a new surface of the excavated area.  My drone surface is about 0.2 ft higher than the 2013 LiDAR surface for those areas not excavated.

Using ArcGIS 10.7 and the Cut and Fill tool in Spatial Analyst I was able to determine that about 9,658 cubic yards (~7,400 cubic meters) were removed. I have not altered or adjusted the surfaces to try to make one fit the other.

I understand that there are uncertainties in any spatial dataset including LiDAR and I am well versed in Structure from Motion reconstructions so I am positive that my drone surface is a good as it gets.

Given the differences and uncertainties between the two LiDAR datasets and the drone surface, how do I estimate the error in my 9,658 cubic yards total?


texts:
I know Tilemill runs under the BSD license, but what about the "base data" that is included with it? I couldn't find any legal information in the guide.
The best thing I could find is this:

"The Mapbox GeoData Library is a collection of free datasets that have been optimized to work well with TileMill. The library is browsable from the TileMill interface through the Add Layer dialog - click on the ‘browse’ button next to the file field, then click the ‘Mapbox’ button with the cloud icon.

All text within the data has been converted to UTF-8 which is the preferred encoding for working with data in TileMill.

All the data has been reprojected to Google Mercator projection, and is properly auto-detected as such by TileMill. All geometry data has also been cropped (where necessary) to fit within the Google Mercator square Earth bounding box of 180.0° east/west and 85.05112877980659° north/south.

All shapefiles have been indexed using Mapnik’s shapeindex tool. (The article Speed Optimization: Shapefile Indexes explains the benefits of this and how you can make sure your own shapefiles are indexed.)"

https://tilemill-project.github.io/tilemill/docs/manual/mapbox-geodata/

But there is no indication or links to licenses, etc. since I intend to use maps commercially (in YouTube videos, which have ads before them), this is quite important, because for commercial application the term "free" isn't enough.


texts: Dividing country into regions based on some initial center point using ArcMap

texts:
I am new to GIS.
I have some lat-long (say 5) spreads across a country and I want to make polygons around each point such that all the polygons cover the whole area of the country and within the polygon, its aerial distance is minimum compared to other points.
A rough illustration is shown in image.


texts:
I am attempting to use XML files, for example Google places API extracted XML files, in ArcGIS 10.4/10.5, and have yet to find an easy way to do so. I do not know what type of XML file it is exactly, I do know that it contains longitude/latitude data. I have attached a pic of the start of the code, if that helps to identify.

So far I have succeeded to convert an XML file to a shapefile through the Add X,Y feature in ArcGIS which turns the file into an Excel file, and then using the FME converting software to convert the Excel file into a working shapefile. This is a rather long way, and takes some time if I want to work with a large number of XML files.

Is there an easier way to do so? Otherwise, is there an easier way to work on spatial features with XML files than ArcGIS?


texts:
I am trying to create a correlation chart between NDVI and LST. Everything working fine except the error message saying Error generating chart: FeatureCollection.randomPoints: Polygon too large to be randomly sampled. Must be smaller than a hemisphere. I am also posting the custom script. Please check the script and tell me where I am wrong.
https://code.earthengine.google.com/de9a4508a6bb69365e3bf63cb455b61c


texts: Simple way to add a legend in QGIS 3.8.2?

texts:
I'm making a simple map using QGIS 3.8.2 and data from OSM.  The maps looks good, but I'd like to add a legend.  It seems like there used to be an "add new legend" button in older versions of QGIS, but I can find nothing like that in 3.8.2.  Is is still possible to add a legend without writing python code?

Addendum: Ok, after seeing ahmadhamb's answer I am very confused.  His screen and menus look very different than mine.  I am attaching a screenshot of what I have, with the opening screen in the background and the "About QGIS" screen in the foreground.  Where do I get the QGIS that ahmadhamb has???


texts: Is it possible to export geospatial PDF from QGIS via GDAL/OGR?

texts: Is there any plugin for QGIS which would help me in exporting the map document as a geospatial PDF (similar to Terrago's GeoPDF)?


texts: Geopackage on S3 to PostGIS RDS using a python Lambda

texts:
I am attempting to write/invoke a Lambda written in python (with a geospatial layer that includes gdal and the usual python geospatial gang: https://github.com/developmentseed/geolambda). This lambda would grab a geopackage from an S3 bucket and immediately upload it to an existing Postgres DB on RDS(POSTGIS enabled) without modification. These files can be quite large (GBs) and I imagine I'll need to stream the body of the S3 object...
I have previously loaded geopackages from my machine to this same database using ogr2ogr. What strategies/libraries are available to perform these operations without calling a subprocess to use ogr2ogr? I could probably do this still, I just couldn't get the subprocess to work initially.


texts: I am completely new at gis but I have to make an add in to work with arcmap to make it zoom to a place I want it to go. (I use a textfield in a php website to tell where to go in arcmap) Is there a zoom function I can use in my case?


texts: Remove the temporary layers from the parameter selection in the processing tool

texts:
I wrote a progressing algorithm using Python. The tool uses multiple other processing tools in the toolbox and creates temporary outputs. Those outputs seem somehow saved while I was running the script while developing. Now, when I call my tool from the toolbox and want to select layer as an input, they appear in the selection even though they are created a while ago and they are not even loaded with QGIS Project that I am working.

This is how the input parameter in my script

self.addParameter(
    QgsProcessingParameterVectorLayer(
        self.GEOLOGIE_LAYER,
        self.tr('Pfad zur geologischen Karte'),
        defaultValue=0,
        optional=False
    )
)

Any idea how I can get rid of them?

I think this part of the script is producing those temporary layers

res2 = processing.run("native:clip", {
'INPUT': vlyr_ger,
'OVERLAY': vlyr_ezgStats,
'OUTPUT': 'TEMPORARY_OUTPUT'})
vlyr_gerClp = res2['OUTPUT']
vlyr_gerClp.setName('vlyr_gerClp')
removeBinaryFields(vlyr_gerClp)
lyrList.append(vlyr_gerClp)


texts: Geoserver GeoTiff Raster Layer Turns Shade Pink (Not Pink Error Tiles)


texts: MySQL select linestrings into a multilinestring

texts:
I have a mysql column of linestring type with 2000 linestrings (high water mark values to be more specific).

I want to merge them into one multilinestring so i can see what the shape of these values are in MySQL workbench spatial viewer or somewhere else.

Is there a mysql function that can join all the linestrings into a linestring or any other way to do this?


texts: Calculating daily precipitation data from CHIRPS over several years using Google Earth Engine

texts:
I am still new to the Google Earth Engine and try to calculate the daily precipitation based on the CHIRPS. I have created a code to calculate monthly precipitation.
Could you advise me on how I change this code to daily precipitation?
var chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY"),
geometry = /* color: #d63000 */ee.Geometry.MultiPoint();

// set start and end year
var startyear = 2018;
var endyear = 2020;

// make a date object
var startdate = ee.Date.fromYMD(startyear, 1, 1);
var enddate = ee.Date.fromYMD(endyear + 1, 1, 1);

// make a list with years
var years = ee.List.sequence(startyear, endyear);
// make a list with months
var months = ee.List.sequence(1, 12);

// Specify Country names
var worldcountries = ee.FeatureCollection('FAO/GAUL/2015/level2');

var filterCountry = ee.Filter.eq('ADM2_NAME', 'Kurgantepa district');

var country = worldcountries.filter(filterCountry);

Map.addLayer(country);
Map.centerObject(country, 6);

// Get the geometry
var country = country.geometry();

var monthlyPrecip =  ee.ImageCollection.fromImages(
  years.map(function (y) {
    return months.map(function(m) {
      var w = chirps.filter(ee.Filter.calendarRange(y, y, 'year'))
                    .filter(ee.Filter.calendarRange(m, m, 'month'))
                    .sum();
      return w.set('year', y)
              .set('month', m)
              .set('system:time_start', ee.Date.fromYMD(y, m, 1));

    });
  }).flatten()
);

var meanMonthlyP =  ee.ImageCollection.fromImages(
  months.map(function (m) {
    var w = monthlyPrecip.filter(ee.Filter.eq('month', m)).mean();
    return w.set('month', m)
            .set('system:time_start',ee.Date.fromYMD(2018, m, 1));
  }).flatten()
);

var title = {
  title: 'Monthly precipitation',
  hAxis: {title: 'Time'},
  vAxis: {title: 'Precipitation (mm)'},
};

var chartMonthly = ui.Chart.image.seriesByRegion({
  imageCollection: meanMonthlyP,
  regions: country,
  reducer: ee.Reducer.mean(),
  band: 'precipitation',
  scale: 5000,
  xProperty: 'system:time_start',
  seriesProperty: 'SITE'
}).setOptions(title)
  .setChartType('ColumnChart');

print(chartMonthly);


texts: Folium and timestamped GeoJSON issue , not reading the data correctly

texts:
I am having issues with Folium and using TimestampedGeoJson. I have the following dataframe structure below. I am trying to display this data with Folium and a time slider to be used on the date field:
    Date        Geometry
0   2016-03-09  POLYGON ((10.00000 25.00000, 10.00000 25.00000, 10.00000
                          26.00000, 12.00000 26.00000, 18.00000 25.00000,
                          17.00000 25.00000))

I am trying to follow this tutorial below and still not working:
https://towardsdatascience.com/visualizing-nyc-bike-data-on-interactive-and-animated-maps-with-folium-plugins-c2d7645cd19b
def create_geojson_features(df):

    features = []
    for _, row in df.iterrows():
        features = []
        feature = {
            'type': 'Features',
            'geometry': {
                'type':'MultiPolygon',
                'coordinates':(row['Geometry'])
            },
            'properties': {
                'times':(pd.to_datetime(row['Date'])),
            }
        }
        features.append(feature)

    return features

df1 = create_geojson_features(df)
world_map = folium.Map(location=[10,0],tiles="cartodbposition",zoom_start=2,max_zoom=6,min_zoom=3)

    TimestampedGeoJson(
          df1
        , period='P1D'
        , add_last_point=True
        , auto_play=False
        , loop=False
        , max_speed=1
        , loop_button=True
        , date_options='YYYY-MM-DD'
        , time_slider_drag_update=True
    ).add_to(world_map)

The geometry and times are not being read correctly via the function and/or the TimestampedGeoJson. Also, in the coordinates field is a shapely object once the function create_geojson_features is executed.


texts: OSM riverbank relations : unable to display them

texts:
I am new in the GIS community and I am having trouble obtaining riverbanks in QGIS. It seems that when a relation is not complete, it is simply not displayed.

And since all riverbanks are now part of a huge relation and since all the nodes involved in this relation are NOT in my extract : QGIS is not displaying it!

Could someone help me? (currently, the only solution I have is to make a recursive request on averpass turbo .... but I can't get it from the planet-latest.osm.pbf file ..... )


texts: Geotools filterfactory not found exception

texts:
Can anyone help me out of this exception.

In a geotools query program i am using

Filter filter = CQL.toFilter("CONTAINS(the_geom,POINT(-115.56,67.567))");

and when am running the program using spark(sbt) ..it is throwing exception as follows can anyone help me....

Exception in thread "main" org.geotools.factory.FactoryNotFoundException: No factory of kind "FilterFactory" found.
        at org.geotools.factory.FactoryRegistry.getServiceProvider(FactoryRegistry.java:375)
        at org.geotools.factory.FactoryCreator.getServiceProvider(FactoryCreator.java:145)
        at org.geotools.factory.CommonFactoryFinder.lookup(CommonFactoryFinder.java:346)
        at org.geotools.factory.CommonFactoryFinder.getFilterFactory(CommonFactoryFinder.java:300)
        at geocode.feature.FeatureExtraction.findShapefile(FeatureExtraction.java:64)
at geocode.batch.GeocodeMain$.main(GeocodeMain.scala:33)
at geocode.batch.GeocodeMain.main(GeocodeMain.scala)
at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
at java.lang.reflect.Method.invoke(Method.java:497)
at org.apache.spark.deploy.SparkSubmit$.org$apache$spark$deploy$SparkSubmit$$runMain(SparkSubmit.scala:672)
at org.apache.spark.deploy.SparkSubmit$.doRunMain$1(SparkSubmit.scala:180)

at org.apache.spark.deploy.SparkSubmit$.submit(SparkSubmit.scala:205)
at org.apache.spark.deploy.SparkSubmit$.main(SparkSubmit.scala:120)
at org.apache.spark.deploy.SparkSubmit.main(SparkSubmit.scala)


texts: distanceTo is giving wrong result in OpenLayers?

texts:
I am using distanceTo from OpenLayers to calculate the distance between two locations

These are map details:

  Geographic  = new OpenLayers.Projection("EPSG:4326");
  Mercator = new OpenLayers.Projection("EPSG:900913");

   map = new OpenLayers.Map({
   div: "map_canvas",
   projection: Mercator,
   displayProjection: Geographic,
   center: new OpenLayers.LonLat(0, 0),
   minResolution: "auto",
   minExtent: new OpenLayers.Bounds(-1, -1, 1, 1),
   maxResolution: "auto",
   maxExtent: new OpenLayers.Bounds(-180, -90, 180, 90),
   units : 'km'
});

these are my lat and lon values

lat = 47.162494, lon = 19.503304
lat = 46.2525684, lon = 20.147060499999952

google maps function result is 111861.5962110549

function distanceBetweenMarkers(latlng1, latlng2){
    return google.maps.geometry.spherical.computeDistanceBetween(latlng1, latlng2);
}

openlayers function result 6237980.800397318

function distanceBetweenMarkers(latlng1, latlng2){
        var point1 = new OpenLayers.Geometry.Point(latlng1.lon, latlng1.lat);
        var point2 = new OpenLayers.Geometry.Point(latlng2.lon, latlng2.lat);
        return point1.distanceTo(point2);
    }

Why is the difference between google maps computeDistanceBetween and openlayers distanceTo?


texts: Testing if geometry is inside another using OGR in Python?

texts:
I need to tell whether geometries from a shapefile are inside other geometries or not.
For example, I want to test if there are trees (represented as points in trees.shp) inside the urban areas (reprensented as polygons in area.shp)

What's the easiest way to proceed?

I loop over the points, and within this loop, I have a nested loop over the polygons.
Then, I test if my point geometry is inside the polygon with the Within(Geometry geo) method from OGR.
It's very long, not optimized at all.

What should I do to fasten the computation?

def PointsInPolygons(self, pointsLayer, polysLayer):
    nbInside = 0
    total = pointsLayer.GetFeatureCount()
    for i in range(total):
        pointFeat = pointsLayer.GetFeature(i)
        pointGeo = pointFeat.GetGeometryRef()

        for j in range(polysLayer.GetFeatureCount()):
            polyFeat = polysLayer.GetFeature(j)
            polyGeo = polyFeat.GetGeometryRef()

            if pointGeo.Within(polyGeo):
                nbInside+=1
                break


texts: Modelling habitat heterogeneity in space: R function

texts:
I'm working on spatial landscape simulations and I would like to simulate a landscape where habitat quality varies in space.

I did already some simulations assuming a landscape with higher habitat amount in the middle and lower habitat amount on the borders of the landscape. In this case I used the gaussian.kernel function in the SpatialEco R package.

Now, I would like to simulate a landscape with higher habitat amount in the borders and lower habitat amount in the middle.

My idea was to use a bivariate truncated bimodal distribution (for instance beta distribution, with U-shape). But I do not know how.
Do you have any suggestion?

I work in R and I prefer to do it in-there.


texts: Finding image overlap in Landsat images - Google Earth Engine

texts:
I am trying to construct a time series of NDVI derived from Landsat 5, 7 and 8 using the Google Earth Engine's Python API. To ensure the validity of this TS, I am harmonising the TM/ETM+ imagery to OLI using the coefficients found by Roy et al. (2016, https://doi.org/10.1016/j.rse.2015.12.024), following the method outlined on https://developers.google.com/earth-engine/tutorials/community/landsat-etm-to-oli-harmonization. As the NDVI TS is pretty central to my research question, I would like to quantitatively assess the performance of this harmonisation using the OLS regression approach used by Roy et al. However, this requires overlapping pixels from different paths at similar dates to be extracted.
Is there a way of extracting this in GEE?


texts: Can overlays have different resolutions from the base layer?

texts:
I am trying to show a regional map on top of base map but without success. While debugging I noticed that my regional map has a different set of resolutions from the base. Is this allowed in OpenLayers?

See http://jsfiddle.net/Sq9hE/15/

If you change the map options to use yorkLayer's resolutions array then the regional tiles will load but no the base map's.


texts:
I have a static list that I would like to pull values from and use in a query. I can populate the list dynamically from the data but the set is quite large and I think having a static list from a csv would increase performatance. I cannot get Arc to accept the return from the input to the select by attribute clause.

class ComboBoxClass1(object):
    """Implementation for addin2_addin.combobox (ComboBox)"""

    def __init__(self):
        self.editable = True
        self.enabled = True
        self.dropdownWidth = 'W' * 16
        self.width = 'W' * 16
        mxd = arcpy.mapping.MapDocument("CURRENT")
        self.mxd = arcpy.mapping.MapDocument("CURRENT")
        self.df = arcpy.mapping.ListDataFrames(mxd)[0]

    def onSelChange(self, selection):
        selection = str(selection).replace("'[",'').replace("]'",'')
        layer = arcpy.mapping.ListLayers(self.mxd, "block_cards", self.df)[0]
        a = arcpy.SelectLayerByAttribute_management(layer, "NEW_SELECTION", "town = '{}'".format(selection))
        self.df.extent = layer.getSelectedExtent()
        arcpy.RefreshActiveView()

    def onFocus(self, focused):
        fldName = "town"
        with open(r'towns_table.csv', 'rb') as f:
            reader = csv.reader(f)
            your_list = list(reader)

        layer = arcpy.mapping.ListLayers(self.mxd, "block_cards", self.df)[0]
        arcpy.SelectLayerByAttribute_management(layer, "CLEAR_SELECTION")
        self.items = sorted(list(your_list))

I have tried several things to strip the brackets from the return (in python everything comes back '[denver]') with no luck. i think thats the sticking point but i could be off on that.


texts: Import a custom function from the python console

texts:
I'm trying to import a custom function from the python console, but it return the error:

'QgsExpressionFunction' object is not callable

My example function (saved as function.py file in /.qgis2/python/ directory) is:

from qgis.core import *
from qgis.utils import qgsfunction
from qgis.core import QgsExpression

@qgsfunction(args="auto", group='Custom')
def TriangleArea(b,h,feature,parent):
    area = b*h/2.0
    return area
    QgsExpression.registerFunction(TriangleArea)

To import and execute the function i use:

from function import *

A = TriangleArea(3, 5)
print A

what's wrong?


texts: How to identify and highlight invalid geometries?

texts: I'm am trying to clean/edit shape files that were created by algorithm in Qgis. My problem is when I try to split some polygons,they have invalid geometry. Is a way to identify these error maybe in a different colour for me to fix them?


texts:
We have built an OSM Tile Server following the instructions at

http://switch2osm.org/serving-tiles/building-a-tile-server-from-packages/

We have imported planet_latest.osm.pbf.
Our server renders very slowly at low zoom levels (especially 1 -10) and increases in speed as the zoom levels increase.

We believe this is to be expected, but could someone please explain why the rendering speeds vastly differ depending on the zoom level?


texts: Adding a geopackage (gpkg) to GeoServer through REST API

texts:
I'm trying to use the REST API to upload a GPKG file (SQLite with vector data) into GeoServer with curl but I always get Error auto-configuring coverage.
The (wrong) command I'm currently using is
curl -v -u <USR>:<PASSWORD> -XPUT -H "Content-type: application/geopackage+vnd.sqlite3" --data-binary @<PATH&FILENAME.GPKG> "http://<IP>:<PORT>/geoserver/rest/workspaces/<MYWORKSPACE>/coveragestores/<FILENAME>/file.geopackage%20(mosaic)"
Thanks to Ian Turton (see the comments below), I noticed that I'm using the wrong endpoint since "geopackage%20(mosaic)" is for raster data only. Which is the right endpoint I have to use with a vector gpkg? I tried searching online but I found nothing...


texts: Selecting Minimum values from in_memory dataset!

texts: I am trying to select the minimum value from a data set being stored in the in_memory workspace. This prevents me from using sub-queries. I'm hoping that someone has found a work around for this issue or something similar.


texts:
I have a PostgreSQL database with one table for each US state. I am trying to make a program which iterates through each of the tables and manipulates the geometry. To do this, at the beginning of the loop I am creating a topology with the line
SELECT topology.CreateTopology('state_topo', find_srid('schema', 'stateTable', ''geom''));
Since I'm iterating I want to delete the topology before restarting, so I've included the line:
SELECT topology.DropTopology(''state_topo'');
The program works fine with just one state table, however when there are multiple I get the following error:
ERROR:  cannot DROP TABLE "edge_data" because it has pending trigger events CONTEXT:  SQL statement "DROP SCHEMA state_topo CASCADE"
How do I fix this bug?


texts: creating many shapefiles from one multi-polygonal shapefile

texts: Does anyone know of a way to create a standalone shapefile from every entry (row) in the attribute table of a large set of polygons (4909 in total)? preferably using python in a way that the results could be stored as a list of values which could be assigned to a variable in an equation later.


texts: Move TEXT_ELEMENT using arcpy.mapping and data driven pages

texts:
I'm customizing an automated mapping solution using arcpy.mapping and data driven pages.
Using a looping statement, I'm hardcoding the position for the mapping elements (legend, north arrow, inset map, etc. so that for each map, the mapping elements do not interfere with the area of interest.  My code works for all mapping elements except for TEXT_ELEMENTS. Text element stays in the first position assigned for all of the maps.

Not quite sure why this is the case, cause I have been able to get the north arrow, a neat line, legend, scale bar, and inset map (dataframe) to move position for 16 different maps.

Sample code below showing how I'm moving the inset map and one of the text elements is shown below:

import arcpy, os
# Create an output directory variable
outDir = r"C:\Desktop\Test"

# Create a new, empty pdf document in the specified output directory
finalpdf_filename = outDir + r"\FinalMapBook.pdf"
if os.path.exists(finalpdf_filename): # Check to see if file already exists, delete if it does
  os.remove(finalpdf_filename)
finalPdf = arcpy.mapping.PDFDocumentCreate(finalpdf_filename)
# Create a Data Driven Pages object from the mxd
mxdPath = r"C:\Desktop\Test\States.mxd"
tempMap = arcpy.mapping.MapDocument(mxdPath)
tempDDP = tempMap.dataDrivenPages
# Create objects for the layout elements that will be moving, e.g., inset data frame, scale text
dataFrame = arcpy.mapping.ListDataFrames(tempMap, "LocatorDF")[0]

for pgIndex in range(1, tempDDP.pageCount + 1, 1):

  # Create a name for the pdf file you will create for each page
  tempMap.dataDrivenPages.currentPageID = pgIndex
  pageName = tempMap.dataDrivenPages.pageRow.Name
  temp_filename = r"C:\Desktop\Test\temp_pdfs\Land_Cover_" + \
                            str(pgIndex) + pageName + ".pdf"
  if os.path.exists(temp_filename):
    os.remove(temp_filename)

  #Code to modify mapping elements for each page

  # Code for Map1
  if (pgIndex == 1):

    #Set inset map position
    inset = arcpy.mapping.ListDataFrames(tempMap, "LocatorDF")[0]
    inset.ElementPositionX = 3.25
    inset.ElementPositionY = 0
    # Move text elements
    elm = arcpy.mapping.ListLayoutElements(tempMap, "TEXT_ELEMENT", "Westport")[0]
    elm.ElementPositionX = 0.3125
    elm.ElementPositionY = 3.0227

  # Code for Map2
  if (pgIndex == 2):

    #Set inset map position
    inset = arcpy.mapping.ListDataFrames(tempMap, "LocatorDF")[0]
    inset.ElementPositionX = 0
    inset.ElementPositionY = 3.25
    # Move text elements
    elm = arcpy.mapping.ListLayoutElements(tempMap, "TEXT_ELEMENT", "Westport")[0]
    elm.ElementPositionX = 2.3125
    elm.ElementPositionY = 2.0227

  # Code for Map3
  if (pgIndex == 3):
    #Set inset map position
    inset = arcpy.mapping.ListDataFrames(tempMap, "LocatorDF")[0]
    inset.ElementPositionX = 3.25
    inset.ElementPositionY = 3.25
    # Move text elements
    elm = arcpy.mapping.ListLayoutElements(tempMap, "TEXT_ELEMENT", "Westport")[0]
    elm.ElementPositionX = 1.3125
    elm.ElementPositionY = 2.0227

  # Code to export current page and add it to mapbook
  tempDDP.exportToPDF(temp_filename, "RANGE", pgIndex)
  finalPdf.appendPages(temp_filename)

  arcpy.RefreshActiveView()

# Clean up
del tempMap

# Update the properties of the final pdf
finalPdf.updateDocProperties(pdf_open_view="USE_THUMBS",
                             pdf_layout="SINGLE_PAGE")

# Save your result
finalPdf.saveAndClose()


texts: I downloaded NDVI greeness data from the USGS Integrated Vegetation Mapping Viewer.  In the metadata it says that visual greenness values range from 0 to 100.  My data is from the period ending 3/5/12 in the NYC area my data table values only range from 0 to 68 which makes sense since it's winter.  I need to display the data in terms of percent, which seems easy enough, but has me hung up.  When I load the raster into the dataframe you can see that the values in the table of contents range from 105 to 255 (the BinValues).  I think I am having trouble understanding the actual data within the raster and how to manipulate it in order to obtain %.  The data is also in GCS (NAD83) and I will have to project it into Albers Equal Area to match my other datasets assigning the cell size to that of another dataset as well (30m) so I have to use a bilinear sampling.  I am a newbie at this so I hope this makes sense.  Thanks.


texts: Get state but want country in GeoPy and Nominatim

texts:
I use the geopy library in Python to find the location of countries. I use country names as input,  like:
geolocator = Nominatim(user_agent="anonymous@gmail.com")
geolocator.geocode("Georgia")

The output: (32.3293809, -83.1137366)
But this is de longitude and latitude of the state Georgia and not the country Georgia.
I tried: geolocator.geocode(country="Georgia"), but that doesn't seem to work; no output.
How do I get the country instead of the state? Am I missing a parameter?


texts: Sentinel-2 Level 2A Product Archive in Google Earth Engine

texts:
My study target region is Eastern Africa (from 2017) and am interested in Sentinel-2 level 2A products from Google Earth Engine (GEE). However, the earliest imagery I could find in GEE for this region is from December 2018. But for France (my comparative study area), the archive dates back to March 2017.

How can I access the Sentinel-2 Level 2A products in Google Earth Engine before December 2018 for this target region?


texts: Truncated time series table output over multiple polygons using Google Earth Engine

texts:
I'm facing issue to export a CSV table of a time series from multiple polygons in GEE. Using few polygons it works well, but when I'm adding more polygons, the output is truncated. My CSV has Polygon's names as column and Time as row names.  The example code can be found here

var S2 = ee.ImageCollection('MODIS/006/MOD11A2')

.filterDate('2018-09-01', '2018-10-31')

.filterBounds(geometry)

.select ('LST_Day_1km');

var S2cell = S2.map(function(img){
  return img. multiply(0.02)
  .subtract(273.15)
  .copyProperties(img,['system:time_start','system:time_end']);
});

var triplets = S2cell.map(function(image) {
  return image.reduceRegions({
    collection: geometry.select(['Name']),
    reducer: ee.Reducer.mean(),
    scale: 1000
  }).filter(ee.Filter.neq('mean', null))
    .map(function(f) {
      return f.set('imageId', image.id()).setGeometry(null);
    });
}).flatten();
print(triplets.first());

var format = function(table, rowId, colId) {
  var rows = table.distinct(rowId);

  var joined = ee.Join.saveAll('matches').apply({
    primary: rows,
    secondary: table,
    condition: ee.Filter.equals({
      leftField: rowId,
      rightField: rowId
    })
  });

  return joined.map(function(row) {
      var values = ee.List(row.get('matches'))

        .map(function(feature) {
          feature = ee.Feature(feature);
          return [feature.get(colId), feature.get('mean')];
        });
      return row.select([rowId]).set(ee.Dictionary(values.flatten()));
    });
};

var table1 = format(triplets, 'imageId', 'Name');

var desc1 = 'table_easimes';
Export.table.toDrive({
  collection: table1,
  description: desc1,
  fileNamePrefix: desc1,
  fileFormat: 'CSV'
});

print(geometry)

As you can see my Shapefile asset has 15 polygons; but the output only shows 8 polygons. The worst is that I would to compute it over 1600 polygons.

What's wrong ?


texts: Showing change of lot division/ownership over time?

texts:
I am new to GIS programs and using ArcGIS Pro.

I am mapping a neighborhood and looking at creating a detailed database of the ownership history. Currently I have the lots configured to have an attribute field for each individual time the lot ownership legally changed, called an Instrument.

This is pretty basic and an ugly solution as it requires me to have 16 fields for lots that have much lower ownership changes. I've been looking into temporal options but it seems to be just for changing the Attributes. I am looking for a solution that allows the polygons in the map to change according to the time slider resulting in something like this in map form:


texts: Do feature classes need to be versioned for editing with ArcGIS Collector/feature service?

texts:
I have an Oracle 18c/10.7.1 geodatabase with feature classes that have the following properties:

Unversioned
SDO_GEOMETRY datatype
Registered with the GDB

My organization plans to serve up the feature classes to a feature service/ArcGIS Server for on-premise/live-editing via ArcGIS Collector (not ArcGIS Online).

Do the feature classes need to be versioned in order for users to edit them via Collector?


texts:
When I want to load a layer style, I can open the layer properties and in symbology I get a list of styles I can choose from. On the left of each style I get a small picture that shows me a preview on how the style looks like. With one klick on the style I get a bigger preview on the top right.

I would love to have a similar drop down menu for my labels. Currently, I have a list of label styles but I can't see how they look like without opening them and applying on my layer.

Is that somehow possible?


texts:
I'm trying to add a WFS/vector layer with OpenLayers 3. So far I've tried

var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: function(extent) {

        var newExt = ol.proj.transformExtent(extent,'EPSG:3857','EPSG:25832');
      return 'http://localhost:8080/geoserver/wfs?service=WFS&' +
          'version=1.1.0&request=GetFeature&typename=postgre:polygon_test&' +
          'outputFormat=application/json&srsname=EPSG:25832&' +
          'bbox=' + newExt.join(',') + ',EPSG:25832';
    },
    strategy: ol.loadingstrategy.bbox
  });

  var vector = new ol.layer.Vector({
    source: vectorSource,
    style: new ol.style.Style({
        fill: new ol.style.Fill({
            color: 'black'
        }),
        stroke: new ol.style.Stroke({
            color: 'blue',
            width: 2
        })
    })
  });

var map = new ol.Map({
target: 'map1',
layers: [
    osmLayer,
    vector
],
view: new ol.View({
center: ol.proj.fromLonLat([10.09, 56.45]),
    maxExtent: new ol.extent.boundingExtent(384015, 5846682, 516991, 5931295),
    zoom: 15
    })
});

If I debug the page with Chrome and look at the network tab I get a reasonable GeoJSON response but nothing is displayed. What am I missing?

I found the reason for the problem. First of all I've missed projection in the format definition. I looks like this now:

format: new ol.format.GeoJSON({defaultDataProjection:'EPSG:25832'})

Furthermore GeoServer adds this:

"crs" : {"type" : "name","properties" : {"name" : "urn:ogc:def:crs:EPSG::25832"}}

at the end of the GeoJSON data. If I download the GeoJSON data and put it in a file on the server and strips the crs-part, the features are added to the map.

Either I need to stop GeoServer from adding CRS to the output or I have to configure OpenLayers to accept it.

Any suggestions?


texts: QGIS crashes when trying to open ecw files

texts:
I added Ubuntugis unstable repo and installed Qgis, libgdal-ecw-src, libecwj2-3.3 but still does not work.

When I attempt to open an ecw file with Qgis it crashes and tells me:
'std::length_error': basic_string::_S_create

I am using ubuntu 13.04 and when I do:

gdalinfo --formats|grep -i ecw

it seems to work.

ECW (rw): ERDAS Compressed Wavelets (SDK 3.x)
JP2ECW (rw+v): ERDAS JPEG2000 (SDK 3.x)

What is the problem?


texts: How do I disable the Add Spatial Database Connection function in ArcGis?

texts: Is it possible to remove/disable the Add Spatial Database Connection and Add OLE DB Connection functionality from ArcMap and ArcCatalog?  Ideally the end user could hit the Add Data button, and the add connections buttons would not show up in the add data dialog.  Any pointers in the right direction are much appreciated.


texts: Is it possible to rename a table column with an expression in the field calcultor of QGIS?

texts:
I want to be able to rename a few column names with an expression in the field calculator in order to be able reproduce this action with the QGIS graphical modeler.

I tried to do this with the refactor field tool of the modeler which works but the information of my columns disappears. Any suggestions ?


texts: Search MXDs with specific dataSource layer using Python

texts:
I work with ArcGIS 10.3 and I try to search a specific layer (with dataSource named "D:\desktop\Project\layers\1.jpg" ) in hundred of mxd's that spread in folder "D:\PROJECTS" and it divided to hundred of sub folders, using python 2.7.8:

import arcpy,os,sys
import arcpy.mapping
from arcpy import env

env.workspace = r"D:\PROJECTS"
for mxdname in arcpy.ListFiles("*.mxd"):
    mxd = arcpy.mapping.MapDocument(r"D:\desktop\Project\\" + mxdname)
    dfList = arcpy.mapping.ListDataFrames(mxd, "*")
    for df in dfList:
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if lyr.isGroupLayer == True: continue
            if lyr.dataSource == r"D:\desktop\Project\layers\1.jpg":
                print mxdname, mxdname.pathway
    mxd.save()
del mxd

Finally, I want that python will print all the mxd's source name that contain the specific layer I search for.
When I run the code I get error:

project1.mxd

Traceback (most recent call last):
  File "C:\Users\yaron.KAYAMOT\Desktop\idle.pyw", line 13, in <module>
    print mxdname, mxdname.pathway
AttributeError: 'unicode' object has no attribute 'pathway'


texts: ArcObjects 10.1 IPixelFilter (IClipFilter) not filtering for IMosaicDataset

texts:
I am attempting to extract an irregular polygon from a Mosaic Dataset using ArcObjects in C#. This sample shows how to open a mosaic dataset in ArcObjects. Referencing existing (VB) code, this snippet is from a working tool which uses IPixelOperation to reduce the raster to only the portion covered by a polygon:

Dim pRaster As IRaster = pRasLayer.Raster
Dim pPixOp As IPixelOperation = pRaster
Dim pShapeWS As IWorkspace = pShapeWSfact.OpenFromFile(pFolderPath, 0)
Dim pShapeFws As IFeatureWorkspace = pShapeWS
Dim pShapeFC As IFeatureClass = pShapeFws.OpenFeatureClass(pFileName & ".shp")
Dim pFeatCur As IFeatureCursor = pShapeFC.Search(Nothing, True)
Dim pFeat As IFeature = pFeatCur.NextFeature
Dim pClipFilter As IClipFilter = New ClipFilterClass()
Dim pGeom As IGeometry

Do Until pFeat Is Nothing
  pGeom = pFeat.ShapeCopy
  If pGeom.SpatialReference.FactoryCode <> pMap.SpatialReference.FactoryCode Then pGeom.Project(pMap.SpatialReference)
  pClipFilter.Add(pGeom)

  pFeat = pFeatCur.NextFeature
Loop
pPixOp.PixelFilter = pClipFilter

Which illustrates how an IClipFilter can be used to restrict the display of a raster... I am attempting to use IClipFilter as IPixelFilter to return the irregular portion of a raster within a polygon:

IMosaicWorkspaceExtensionHelper mosaicExtHelper = new MosaicWorkspaceExtensionHelperClass();
IMosaicWorkspaceExtension mosaicExt = mosaicExtHelper.FindExtension(gMosDSsrcWorkspace);
IMosaicDataset pMosDS = mosaicExt.OpenMosaicDataset(gMosDSname);
IRasterDataset3 pInDS = (IRasterDataset3)pMosDS;
IRaster pInRaster = pInDS.CreateFullRaster();
IPixelOperation pPixOp = (IPixelOperation)pInRaster;

IClipFilter pClipFilter = new ClipFilterClass();
// add some geometries to the clip filter, in the same spatial reference as the raster
pPixOp.PixelFilter = (IPixelFilter) pClipFilter;

// then read/write the raster
IPnt ReadBlockSize = new PntClass();
IPnt SrcUpperLeft = new PntClass();
IPnt DestUpperLeft = new PntClass();

// create the pixel block for read/write
ReadBlockSize.SetCoords(pDestCols, RowsToRead);
IPixelBlock pReadBlock = pInRaster.CreatePixelBlock(ReadBlockSize);
IRasterEdit pWriteEdit = (IRasterEdit)pWriteRaster;
for (int RowMin = 0; RowMin < pDestRows; RowMin += RowsToRead)
{
    SrcUpperLeft.SetCoords(pSrcColOffset, pSrcRowOffset + RowMin);
    DestUpperLeft.SetCoords(0, RowMin);
    pInRaster.Read(SrcUpperLeft, pReadBlock);
    //(pClipFilter as IPixelFilter).Filter(pReadBlock);// tried this too, but it still didn't work
    pWriteEdit.Write(DestUpperLeft, pReadBlock);
}

The result:

The image is the size of the dataset which is correct but the pixels are unfiltered. Does IPixelOperation, IClipFilter only work in ArcMap or with IRasterLayer? Do I need to supply the geometries in cell coordinates? Is there an updated/different interface required to work with an IMosaicDataset?... I've tried several different combinations of readers/writers/filterers and none of them seem to have any affect - What am I missing?


texts:
I have the UTM code 36S and want to convert it to EPSG.
Somebody told me that if the UTM code contains 'S', then the first three digits of the EPSG code are 327, and if the UTM code is 'N' then use 326 instead.  And then you just append the UTM number (in this case 36) to get '327' + '36' => '32736'.
Does this always work?  And should I be using some python library to do this for me?


texts: Programmatically rematch unmatched addresses using a different composite locator

texts:
I have a table of about 1500 addresses that updates and geocodes nightly using task scheduler and a couple of python scripts. The majority of these addresses are in the county for which I work. We have a local composite locator that can get all these addresses no problem, but the locator fails to match the addresses outside the county because it is built on county-wide datasets. The addresses that do occur outside the county range all over the United States.

So, instead I have it set up to connect to and use the ESRI World Geocoding service to geocode the addresses each night, which works like a charm. However, with 1500+ addresses geocoding nightly, this is very credit-intensive, and since we can geocode the majority of them without using the service, I am looking for a way to use our local locator first, and then rematch the unmatched (out of county) addresses using the ESRI World Geocoding service. Any thoughts?


texts:
I have a GeoJSON dataset of the following format (LineString features):

{
  "type": "FeatureCollection",
  "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

  "features": [
    { "type": "Feature", "properties": { "id": 433775, "clazz": 72, "flags": 7, "source": 26697, "target": 48416, "kmh": 10, "time": 1 }, "geometry": { "type": "LineString", "coordinates": [ [ 13.3782144, 52.516451399999973, 0.0 ], [ 13.3782407, 52.5163161, 0.0 ] ] } },
    { /* [... around 13k features ...] */ },
    { "type": "Feature", "properties": { "id": 377768, "clazz": 51, "flags": 7, "source": 270569, "target": 270570, "kmh": 5, "time": 600 }, "geometry": { "type": "LineString", "coordinates": [ [ 13.388153799999985, 52.492538299999971, 0.0 ], [ 13.3881033, 52.4924441, 0.0 ] ] } }
  ]
}

For a research project I need this featureset to be cut into tiles of different zoom levels.

without generalization / simplification of geometries
without caching and writing results to new GeoJSON files on disk

I had a look at geojson-vt by mapbox and it seams it only generates the tiles on the fly, but I need them on the disk for further processing and investigation.

Is that possible? If so, how?

Updates:

I'm on ArchLinux, command line tools like ogr2ogr or any C/C++ tools would be the best, regarding speed.
I already tried tilestache but it is totally broken and not well maintainted. Couldn't get it running.
Input aswell as output should be GeoJSON (or any vector format).
I'm talking about tiles a la Google Maps. The zoom factor for the resulting tiles should be 12, 13 or 14.

The result should be something like this:

./12/2200/1343.json
./12/2200/1344.json
[...]
./14/8800/5372.json
./14/8800/5373.json

A collection of GeoJSON tiles cut from the original vector data.


texts: Unable to truly polygonize raster in QGIS

texts:
I've used QGIS's Polygonize to convert the following raster into outlines. The polygonize tool does produce a vector output, but it isn't in plain polygon form but some thousands of single pixel polygons.

I understand that this is the expected behavior of the Polygonize tool and one can now simplify these weird squares for it to just be two polygons - one for brown and one for purple. I tried setting a really low threshold of 0.000001 for v.generalize (in the Processing toolbox) for these shapes to be converted into a line and yet there seems to be no improvement. The documentation of v.generalize did not seem too useful.

How does one go about making this into two plain polygons?


texts: Loading vector tile layer in Leaflet map?

texts:
I need to load a vector tile layer in a Leaflet map.

The vector tile is the vector tile layer about Mapillary sequences (look at https://a.mapillary.com/#vector-tiles ...), and the tile URL pattern is:

https://d2munx5tg0hw47.cloudfront.net/tiles/{z}/{x}/{y}.mapbox

the vector tile use the Mapbox vector tile format.

I've searched on the net but I haven't found a sample: I've seen that this could be done using Mapbox, but if it's possible I'd like to use only Leaflet


texts: Transpose FeatureCollection export Google Earth Engine

texts:
I'm exporting the pixel values from a set of points (n = 10,000) across the Landsat collection. At the moment the code exports each point as a new row; however, I would like each point to be a new column. Is is possible to transpose the FeatureCollection prior to the export

Link to the export code

The main code snippet is

var mapfunc = function(feat) {
  var geom = feat.geometry();
  var addProp = function(img, f) {
    var newf = ee.Feature(f);
    var date = img.date().format('yyyy-MM-dd');
    var bands = img.bandNames();
    var value = img.select(bands)
    .reduceRegion(ee.Reducer.first(), geom, 30);

    var fun = function(name, ini) {
      var ini = ee.Feature(ini)
      var val = value.get(name);
      var newname = ee.String(date).cat('-').cat(name);
      return ee.Feature(ee.Algorithms.If(val,
                        ini.set(newname, ee.String(val)),
                        ini.set(newname, ee.String('No data'))));
    };

    return ee.Feature(bands.iterate(fun, newf));

  };
  var newfeat = ee.Feature(overall_collection.iterate(addProp, feat));
  return newfeat;
};

var site_data = stable_points.map(mapfunc);


texts: Clipping raster with mask in QGIS gives AAIGrid error?

texts:
I am using QGIS 2.18.4 Las Palmas on Windows 10 (32 bits).

I need to clip raster (TIF) 1 km resolution imagery files downloaded from MODIS using a shapefile mask of Argentina. To do so I tried to use the Clipper tool in the GDAL plugin. All of my layers are in the same coordinate reference system and I have the "Enable on the fly CRS transformation" option disabled.

After doing Raster >> Extraction >> Clipper >> Clipping mode: mask layer, I get the following error message:

  Output driver 'AAIGrid' not recognised or does not support direct
  output file creation. The following format drivers are configured and
  support direct output: ...

followed by a list of 45 supported outputs.

After a thorough online search of related problems and potential solutions, I feel like I have exhausted all possibilities.

Could someone point me in the right direction please?


texts: Update virtual field on demand (not permanently) (QGIS)

texts:
I quite often use virtual fields. Some of them are relatively complex aggregates. Since virtual fields updates on every input, by moving the map, open the attribute table and so on, it can cost a lot of CPU-Resources.
Is there a way to suppress the permanent update of a virtual field or just update it on demand? F.e. it would be fine, if virtual fields updates just by saving the layer or ending the edit mode.


texts: How to select two points on a map using QGIS?

texts: I need to select two points on a map in QGIS.  I selected the element selection tool in menu view and I want to select two points.  How can I accomplish this? Are there other tools or methods to select points in QGIS?


texts: Label placement for cadastral layer in QGIS

texts:
Using QGIS 2.14 and 2.18

Is it possible to place the "label" for the cadastral polygons based on the property shape? Meaning, based on the objects geometry.

The default placement is either horizontal or vertical (mostly I'm not sure how else to do this)

Ideally, I would like to place the labels as shown below

Is this possible given so much of improvement in the QGIS labeling engine in the recent releases.

If possible, can someone please provide the details?


texts: Re-projecting base map (WGS 84 EPSG: 4326)

texts:
I am trying to change my base map from the open street plug in to Google Street Map from WGS 84 EPSG: 3857 to WGS 84 EPSG: 4326.

My spatialite polygons (originally mapped on United States) projects into the ocean.
Is there a way to project the spatialite polygons (mapped in WGS 84 EPSG: 4326 setting) to a the new basemap in WGS 84 EPSG: 4326 setting?

Here is the properties image for reference:


texts: Attribute table in QGIS won't save custom values

texts: I am trying to edit a point shapefile in QGIS and I have editing on for that layer. I have added a new column for the custom data, I have set it to text and have 8 characters allowed in the column, but every time I save and toggle editing off all of the values go back to null.


texts: Use of Grass and R w.r.t. Land Use Prediction

texts:
In my data mining term project I need to conduct a work on geospatial data related to a big city to show possible urbanization fields in next years. I've used GRASS lately on other projects so I'm a bit familiar with GRASS.

After some research I've seen that R can be used for spatial statistics issues and I've read about forecast, raster, sp, spatstat, geospt etc libraries of R and it seems it can be helpful for my project. But on the other hand it has a very steep learning curve. Though I'm not sure that GRASS is whether sufficient or not at that point, learning process of R is seems scaring to me.

I'm coding for 4 years with C, Java, C#, Python etc but R's syntax, working mechanism is like an alien to me. What do you suggest for this problem. Do you think GRASS can be used solely or do i need to invest on R?

Thanks in advance.

EDIT: I have etm+ images produced at 3 different years(2010, 2011, 2012) for same location.
I will use first 2 for training and last one for test.
I will define parameters such as city centers, roads, urbanized places etc. to predict areas which are probable to be urbanized.
There are some researches using Bayesian Belief Networks or Artificial Neural Networks to achieve this. I will use 2 or 3 parameters(roads, rivers) to estimate urbanization model of a city. Using vector data and raster data I think I can predict urbanization model more or less. So I need to evaluate raster cells against vector data to calculate their urbanization probability.


texts: RuleBasedRenderer with two different symbols using PyQGIS

texts:
I am trying to create a RuleBasedRenderer with two layers. The first layer is used to fill the polygons and the other should create a hatching effect with lines. However, the code below, creates the layers needed, but both line color and fill color are the same. Do you know how can solve this?

        symbol = QgsSymbol.defaultSymbol(layer.geometryType())
        symbol.deleteSymbolLayer(0)
        symbolLayer2 = QgsSimpleFillSymbolLayer()
        symbolLayer1 = QgsLinePatternFillSymbolLayer()
        symbol.appendSymbolLayer(symbolLayer1)
        symbol.appendSymbolLayer(symbolLayer2)
        renderer = QgsRuleBasedRenderer(symbol)
        # get the "root" rule
        root_rule = renderer.rootRule()

        for label, expression, color_name in rules:
            # create a clone (i.e. a copy) of the default rule
            rule = root_rule.children()[0].clone()
            # set the label, expression and color
            rule.setLabel(label)
            rule.setFilterExpression(expression)
            if color_name == 'orange':
                print(len(rule.symbols()))
                rule.symbol().setColor(QColor(color_name))

            else: rule.symbol().setColor(QColor(color_name))
            root_rule.appendChild(rule)

        # delete the default rule
        root_rule.removeChildAt(0)

        # apply the renderer to the layer
        layer.setRenderer(renderer)
        layer.triggerRepaint()

P.S. The if statement is for testing purposes.


texts: $area calculation doesn't work or wrong projection?

texts:
I've been reading the QGIS manual and I've encounter a problem I don't know what to do.
I downloaded data from OpenStreetMap and the raster from http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp
The thing is that when I try to calculate the area it turns in REALLY small numbers such as 6.18e-06
I've tried with a couple of projections and it doesn't seem to work.

Any ideas?


texts: Creating network for pgRouting from set of multiline strings?

texts:
I have a table which is a set of multiline strings.

I wish to create a pgrouting topology from that table so I can do shortest-path problems.

I have tried pgr_createTopology but it failed.

Is there a good tutorial on this?


texts: Upload KML / GML / GeoJson to Geoserver through Rest API

texts:
What is the best option, if any, to upload a vector format like GeoJson, and or GML to Geoserver?

The environment from which to upload will be a client web-application, so i cannot use curl. In the backend, i am using a C#-Proxy to connect to Geoserver.

Do I have to convert the file to shape (*.shp) first in some way, which i am able to upload via REST, or is there an option to load the data directly?

On the docs, there is an option for uploading property files, but i cannot find any documentation on it and am not shure about how to use that or if it can help solving my problems.


texts: Attributes encoding problem with rgdal::writeOGR

texts:
I performed an intersection of two SpatialPolygonsDataFrame on R (3.5.1) using raster::intersect.
Here's the structure of the resulting SPDF attributes table

> str(shp@data)
'data.frame':   944 obs. of  7 variables:
 $ COD5      : chr  "1122" "1122" "1122" "1122" ...
 $ LIV_1     : chr  "1" "1" "1" "1" ...
 $ LIV_2     : chr  "1" "1" "1" "1" ...
 $ LIV_3     : chr  "2" "2" "2" "2" ...
 $ LIV_4     : chr  "2" "2" "2" "2" ...
 $ LIV_5     : chr  NA NA NA NA ...
 $ SHAPE_AREA: num  23499 10474 5227 13752 71721 ...

I exported the SPDF with rgdal::writeOGR

writeOGR(shp, 'C:/foo', 'shp', 'ESRI Shapefile', overwrite_layer=T)

but QGIS 3.4 LTR was unable to display attributes correctly

Is it an enconding problem with rgdal::writeOGR? I read the documentation of the function but I didn't figure out how to solve it.
I tried with

writeOGR(shp, 'C:/foo', 'shp', 'ESRI Shapefile', layer_options="ENCODING=UTF-8", overwrite_layer=T)

but the result didn't change. I also tried with the updated rgdal package version (1.4-4).


texts: Export to .dbf changes my data in ArcMap 10.2.2

texts: I have an excel table with parcel ID's (PINs).  The PINs are in "number" format in excel.  Iimported this excel table into ArcMap with the goal of joining them with the county parcel shapefile.  I right click on the excel table>export>navigate to myfolder and save as a .dbf.  I add the .dbf to my map and the data is not correct.  When I export, my numbers change. Does anyone know how to fix this?


texts:
I have a polyline shapefile containing a number of lines, some of which overlap. By overlap, I mean that they follow exactly the same path for some of their route, but also have bits that are not shared. As it's difficult to explain in words, the following images should help:

In the example above, you can see how the two selected lines overlap, and I want to merge them into one line feature.

Is there a way to do this with a built-in ArcGIS tool? If not, is there an easy way to check for overlapping in Python? (and then an easy way to merge two lines in python?)


texts: Filter features based on their attributes using Python

texts:
How to get features by their attributes (similar to 'IQueryFilter' in ArcObjects) in QGIS using Python? Instead of getting all features and filter it manually, is there any option to use where clause to filter it out?
Example: I have a field name called "Counties". It has more than fifty thousand of features i.e not possible to fetch all features and filter it because of time consuming. So i can query it by using iqueryfilter.whereclause = 'Counties = Norwich' in ArcObjects. Similar thing I need in PyQGIS.


texts: Using QGIS raster calculator for expression of height with RGB?

texts:
I'm using QGIS 3.4.1, and downloaded the SRTM TIFF data for the area of my interest, and added the TIFF format SRTM file to my project as a raster layer.

Here, I wish to make my own color gradient for SRTM raster layer using RGB values.

Red and green will describe the color every 255m. Basic equation for R and G is height = R + (G*255). Here, if height is 500m, this will be 245 + (1*255), so RGB value will be (245, 1, 24). If height is 1000m, this will be 235 + (3*255), so RGB value will be (235, 3, 40). Value for blue is explained below.
Blue will be used as some sort of "category" of height. B increase 8 for every 255m from -512m (it starts from -512m), up to B=255. For example, From -512m~-256m height, B=0. From -255m~0m height, B=8. From 1m~255m, B=16. From 256m~512m, B=24... Accordingly, above example of 500m should be (245, 1, 24), and 1000m should be (235, 3, 40)

This way, I wish to create my own color gradient scheme, to express height using specific RGB numbers. I'm trying raster calculator now, but as a beginner, making equation is not easy for me.

What should be the variable names for R, G, and B values?

I have one raster layer name "SRTM@1". I was thinking of equations like "SRTM@1" = "Red" + 255*"Green", but of course it is not working. I'm not sure if this is correct expression for RGB, and if my logic behind this equation is correct from raster calculator.


texts: Read and then write rasterio GeoTIFF file without loading all data into memory

texts:
I work with large drone image raster files, on the order of 40,000 x 40,000 and above. I have a large uncompressed GeoTIFF file and I want to use rasterio to rewrite the file in compressed format. I can do this by loading all data into memory, but is there a way to execute this writing without loading everything to memory?
For my original file I can open it with:
 dat = rasterio.open("grid_001.tif")

Then to rewrite the file with compression I tried:
profile = dat.profile.copy()
profile.update(
        compress='lzw')

with rasterio.open("grid_001_compressed.tif", 'w', **profile) as dst:
        dst.write(dat)

This will give me an error:
 ValueError: Source shape (44134, 44162) is inconsistent with given indexes 3

This is an expected error, because when I open the dataset it creates an iterator or lazy object without actually accessing the data. Now, if I did a command like:
 dat = dat.read()

This will load all of the data from the file into memory, and I can an array of dimensions [3, 44134, 44162] that I can write. BUT, this takes a lot of memory to implement.
Hence, is there a way to perform the same operation, but without loading everything into memory? I am not sure if windowed reads would help in this case or anything.


texts:
Some weather providers like Yahoo give a visibility percentage while others like OpenWeather.org give cloudiness percentage.

To normalize would it be acceptable to say that visibility = 100 - cloudiness?