texts: Exporting training data for deep learning tool in ArcGIS Pro?

texts:
I am trying to export training data for deep learning so I used ArcGIS Pro.

I use sentinel imagery and it exports just 2 or 4 images and these images not open.

I try to use 3 band composite google earth imagery also.

ArcGIS Pro exports more images but all of them are black and do not show the object.

What is the problem?


texts: Raster Buffer error with package updates

texts:
I have just update my r packages and downloaded the newest version of R, and now the buffer function in the raster package no longer works as before.

An example of my issue:
Running this example code: with updated packages and R version 3.4.3 (2017-11-30) -- "Kite-Eating Tree"
library(raster) #version 2.6-7
The Rcpp package was also updated to 0.12.14

    r <- raster(ncols=150, nrows=150, xmn=0)
    r[] <- NA
    r[11050:11075]<-1
    plot(r)
    b<-buffer(r,width=600000)
    plot(b)

Running this same piece of code on another computer with R version 3.3.3 (2017-03-06) -- "Another Canoe"
library(raster) #version 2.5-8 I get the result I am looking for...


texts:
I have created Cache imagery using the Manage Tile Cache tool of ArcGIS Desktop. This tool creates cache imagery from the data provided. After creating this cache, i can open that in ArcMap easily from ArcCatalog, just like a raster imagery.

How can I open that cache created dataset in QGIS?

I want to open that Cache dataset in QGIS, and view it.


texts:
How do i create a new mapset from QGIS 3.2 for GRASS?

I can´t find the Plugin GRASS in manage and install plugins...


texts: Converting Shapefile elevation field to height field

texts:
Have been trying to update the elevation field in my shapefile (vector file *.shp) into a height field. Have tried to rasterize the shapefile then deduct the height using Raster Calculator with a reference dted (raster map). Then I filled the empty space with nodata. After which I polygonize.

I have managed to get some resemblance of the height but however, I have found that the height in the shapefile is rounded off and the processing process actually cost me to lose some features.

All of this is done in  QGIS 2.2

Is there a better way I can do this?

One of the way I am trying to explore is to use field Calculators in vector shapefile properties. Can I reference an external raster map layer with this? If so, how can I do this?


texts: Troubleshooting with Tiger for d3 -- is this the wrong projection? Why's it look like this?

texts:
UPDATE: Solved, sort of … I downloaded the file afresh and tried feeding it through ogr2ogr and topojson again … and wouldn't you know … after a while it worked. I think I was possibly not setting an id property.  This was a big help: https://stackoverflow.com/questions/14565963/topojson-for-congressional-districts

ORIGINAL POST:
So, for the friends who have used Tiger shapefiles, clue me in on my mistake here?

I suspect it's something with wrong projection or wrong conversion from shapefiles?  The .xml suggests the projection is WGS_1984_Web_Mercator_Auxiliary_Sphere

I converted with ogr2ogr then to topojson, a la the Mike Bostock canonical example: http://bost.ocks.org/mike/map/
(FWIW, I also tried with mapshaper.org and also got a bunch of shards.)

UPDATE: I wasn't clear abt mapshaper.org -- it displays this file perfectly online, but when I download topo from it, yea, I get shards.

This should be a map of two Georgia counties, an area maybe 100 miles across. The shapefile came from U.S. Census Tiger.

I got another shapefile of the same area from our state mapping agency and it mapped perfectly.

Any ideas?

UPDATE: here is the shapefile I've been using: https://www.dropbox.com/s/kux63jdo66h87uk/reference_map_shape.zip?dl=0

d3.json("050_002.json", function(error, counties) {
  if (error) return console.error(error);
  console.log("ok counties:", counties)

  var projection = d3.geo.mercator()
    .scale(100)
    .rotate([84.2, -33.9, 0]);

var path = d3.geo.path()
  .projection(projection);

svg.append("path")
  .datum(topojson.feature(counties, counties.objects['050_00']))
  .attr("d", path)


texts:
For ArcGIS for Desktop (Concurrent) users how do you add multiple computers so that they can also use ArcGIS for Desktop as well? Is it through the License Manager, ArcGIS Administrator, or is there some other way of establishing a network connection?

Any help or suggestions would be beneficial.

Thanks.


texts: How to improve the accuracy of geolocation on openlayers 5

texts:
I wonder if it is possible to improve the accuracy of geolocation on openlayers 5, currently it is off by 5000 meters although I enabled highaccuracy option.

   var geolocation = new ol.Geolocation({
  // take the projection to use from the map's view
    trackingOptions: {
          enableHighAccuracy: true
        },
    projection: view.getProjection()
});


texts:
I created a raster file using interpolation (bspline and idw). The original vector data were GPS elevation data I collected on the ground (in Greece). After interpolation, the file contains information such as this:

north: 38.2179846
south: 38.2173841
east: 22.1286094
west: 22.1281249
rows: 499
cols: 403
21.031 21.032 21.033 21.034 21.035 21.036 21.037 21.038 21.039 21.039 21.04 21.041 21.042 21.043 21.043 21.044 21.045 21.045 21.046 21.046 21.047 21.047 21.048 21.048 21.049 21.049 21.05 21.05 21.05 21.051 21.051 21.051 21.051 21.051 21.051 21.052 21.052 21.052 21.052 21.052 21.052 21.051 21.051 21.051 21.051

... etc.

So, everything in the file looks ok and it displays well in 2D:

The problem is when I change the view from 2D to 3D view nothing shows at all. I've tried zooming in/out, full extent, but nothing happens. Could it be a problem of scale? And if so, how one would go about solving it? Thank you for any pointers.

This is the output of computational region:


texts: Using WEKA J48 classifier for satellite image classifications?

texts:
I am trying to build a tool from WEKA software using the J48 classifier for automatic satellite image classification to build landuse land cover maps.

How do I input the satellite image in WEKA after converting it into ASCII format?


texts: Do I need a database for generating WFS-T with QGIS-Server?

texts:
I have hard difficulties in understanding: I have a QGIS-Project (in QGIS) and I want to create a WFS-T by using QGIS-Server (and the QGIS-Server I have is connected with e.g. Windows-Server).

Is the Data for the WFS-T stored in the QGIS-Server?
Or do I need for WFS-T a DB, where QGIS-Server grabs the data?

The QGIS-Documentation said nothing about a DB, instead it says "QGIS as Dataserver". But this Blog and others always talk about QGIS-Server in connection with a DB.


texts: Converting (Polyline Data) JSON to GeoDataFrame using GeoPandas

texts:
I have this particular json and I want to make the JSON into a GeoDataFrame (with geometry and its attributes); How to set the geometry and attributes to a geodataframe?
I have tried using geopandas.read_file() and it returned an error: "expected string or bytes-like object";

I also tried geopandas.GeoDataFrame.from_features() but it returns this table (which is the problem):

I want it to look like something like this (a geometry column, and the other attribute columns):

I am aware of Esri's ArcGIS API for Python but, in this case, I can't use ArcGIS API or ArcPy.


texts:
In QGIS 2.18.16 I'm using the following code (and it works).

import requests
service_url = "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}"
service_uri = "type=xyz&zmin=0&zmax=21&url="+requests.utils.quote(service_url)
tms_layer = iface.addRasterLayer(service_uri, "Google Sat", "wms")

How should I modify it to work in the QGIS 2.99?


texts: Hiding NoData values of a raster in MapServer

texts:
In a MapServer NoData value pixels can be hide by using the following code inlcusion in a mapfile:

`PROCESSING 'BANDS=3,2,1'`

Now I have raster with multiple combinations of NoData Values other than "RGB:000".
I have tried with combinations like OFFSITE 15 20 106 but its not reflecting in the output. What is the way to offsite the NoData Values of multiple combinations.


texts: Using r.mapcalculator in Python console?

texts:
I'm trying to use the r.mapcalculator algorithim in a python script, but I'm having problems in using the r.mapcalculator to tranform a raster layer with float values into integers.

Using the processing toolbox of QGIS, the algorithm works flawlessly. However, using the python code below, the values don't get rounded, instead, I only get NaN values.

processing.runalg("grass:r.mapcalculator","D:/OneDrive/reclassified maps/teste.tif", None, None, None, None, None, "round(A)", None, None, "D:/OneDrive/reclassified maps/teste_arredondado.tif")

EDIT: Just solved the problem. The last None parameters cannot be "None".
So, for the GRASS_REGION_PARAMETER I used the extent() comand to create a string.

ext = teste.extent()
a = str(ext.xMinimum())
b = str(ext.xMaximum())
c = str(ext.yMinimum())
d = str(ext.yMaximum())
size = a + ',' + b + ',' + c + ',' + d

The last parameter (GRASS_REGION_CELLSIZE_PARAMETER) just input 0. The final form is

processing.runalg("grass:r.mapcalculator",
teste, None, None, None, None, None,
"round(A)", size, '0',
"D:/OneDrive/Doutoramento/GIS/Aveiro Coelho (2005)/reclassified maps/teste_arredondado.tif")


texts: Distance to object, taking into account size of object

texts:
I have a number of polygon features on the edge of a lake feature. I also have sample points throughout the lake feature. I am looking to add some sort of indicator for these points as to how influences they are by a nearby wetland. For example, being 100m away from a large wetland should give a different value compared to 100m from a small wetland. In addition, 100m away from two small wetlands should give the same as being 100m away from a single 2x as large wetland.

Is there a single tool for this or do I need to do a calculation on my own. E.g. find distance, find size of nearest wetland and multiply.

Thinking maybe to write a program that takes distW1 * areaW1 + distW2 * areaW2 etc...


texts: Create a shapefile from inMemory workspace

texts:
I am trying to convert a raster dataset into a shapefile using IConversionOp.RasterDataToPolygonFeatureData API.

See the code below:

IConversionOp conversionOp = new RasterConversionOp() as IConversionOp;
IWorkspaceFactory inMemoryWorkSpaceFactory = new InMemoryWorkspaceFactory();
IWorkspaceName inMemoryWorkspaceName = inMemoryWorkSpaceFactory.Create("", "MyWorkspace", null, 0);
IName name = (IName)inMemoryWorkspaceName;
IWorkspace inMemoryWorkspace = (IWorkspace)name.Open();

conversionOp.RasterDataToPolygonFeatureData(rasterDataset as IGeoDataset,
inMemoryWorkspace, inMemoryWorkspace.PathName, false);

I am using InMemory as RasterDataToPolygonFeatuerData fails when the path argument has spaces and I do have fair chances to have spaces in the output path.

So I tried to create the FeatureClass using InMemory workspace. The FeatureClass gets created successfully and now I need to create the Shapefile on disk.

Can someone help me in creating the Shapefile on a disk from InMemory workspace?


texts: Is there any API in C# for Zooming selected feature?

texts:
I am new to ArcGIS and want to implement the Zoom to Selected Features functionality available in Selection menu option (Selection-> Zoom to Selected Features), using C# API.

Is any API available in C# ?


texts: QGIS - Moving Line Features Large Distances

texts:
I have two layers on my map: reserve boundaries (polygons) and reserve trails (lines). I am using WGS 84 (EPSG:4326) and projected my trails from WGS 84 / UTM Zone 17N (EPSG:32617) to match my polygon and project CRS.
However, my trails layer is not aligned with my reserve boundaries layer and the distance between them is too great to be moved manually using the editor tools. Following up online I found the "Move Features" plugin, but I cannot get it to perform the transformations. Does anyone have a solution?


texts:
What I want to do appears simple. However, there seems to be a bug.

I have an image from GE which I want to georeference and set as a background layer. I thus set some placemarks on GE, wrote down the coordinates, created a picture and uploaded it to the georeferencer. I then wanted to put in the coordinates for the placemarks manually, however, they do not fit. In order to see which would be kinf of fitting coordinates, I chose the "from map canvas option" and got totally different coordinates for which I can't make out the unit at all. I have specified WGS 84 everywhere as coordinate system and I also have some other layers on QGIS, which display the correct lat/ long values when asked by the information function.

Do you know what mistake I could have made?


texts: Exporting Raster with Applied Color Ramp with ArcPy?

texts:
There is NDVI image of a small agricultural area. I was also given separate shapefiles for all fields (37 fields) inside that area. I have written a script that clipped NDVI image for each field using the shapefile. When a clipped NDVI image is saved, it is exported in tiff format and grayscale.

Shape_Files_Name = ['TX1_10_poly.shp', 'TX1_11N_poly.shp', 'TX1_11S_poly.shp'...........]

for n in Shape_Files_Name:

    #mask out each field
    Clipped_File_Name = random.randint(1,101)
    arcpy.gp.ExtractByMask_sa (NDVI_Layer, n, Clipped_File_Name )

    #save the raster
    arcpy.CopyRaster_management (Clipped_File_Name, NDVI_Save_Location_and_Extension)

I was trying to find a way to export them to a known colour ramp.  Thus,  a dummy layer is created and applied the desired colour ramp to it. Later, for each NDVI field image, just applying the colour ramp by the following codes

arcpy.ApplySymbologyFromLayer_management(in raster, color ramp layer)

It works. But I would like to export those NDVI images in the colour ramp I applied.

Is there any way to do that in ArcPy?


texts: How can i detect the number of polygon is drawed in a edit session in QGIS

texts:
I am working with qgis 1.8 and i want to know how many polygon drawed in one edit sessions
 is there any signal emitted with qgis ???


texts: Seeking detailed description of qgs file structure?

texts:
Is there a detailed description on the .qgs file structure somewhere?

I checked at http://mrcc.com/qgis.dtd and got

  The requested URL /qgis.dtd was not found on this server.

and at https://svn.osgeo.org/qgis/trunk/qgis/qgis.dtd which wasn't very detailed.


texts: Is it possible to change the year range of time manager in QGIS?

texts:
I have a shapefile with a start and end time field, properly formatted. However my dates range from 1864-1995. When I try to use the shape file with the QGIS time manager plugin I keep getting this message:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/PROGRA~2/Quantum GIS Wroclaw/apps/qgis/./python/plugins\TimeManager
   \timemanagerguicontrol.py", line 342, in updateTimeExtents
    self.dock.horizontalTimeSlider.setMinimum(mktime(timeExtents[0].timetuple()))
ValueError: year out of range

I found that one limitation of time manager is: The range of manageable years is limited due to limitations in time.mktime. The exact range is platform dependent.

Is there a way to find out what the exact range is or change it?

When I tested a sample of my shapefile for 1977 and 1984 they worked fine.


texts:
I managed to scrape an .mvt (mapbox vector tile) file from an online interactive map: here. I procured a landuse development layer. What I want to do is, to reproject this to EPSG:4326 WGS 84 so that I can combine it with other mumbai shapefiles to create a landuse map with transport lines etc. From the site where I scraped, I am unable to get any projection information. Internet tells me that mapbox vector tiles usually use the EPSG:3857 projection, which I have tried to set in QGIS, but that appears to be incorrect too. I also cannot find any resources online to deal with converting the mvt geometry into coordinate systems.
Any tips?
I am new at this.  I have attached a sample .mvt file from the site here for reference.


texts: Remove duplicate features based on attributes using ModelBuilder?

texts: I need to create a model in the ArcGIS ModelBuilder. I am working with a feature class of zip-codes for a state. I have got to the point where the zip-codes that have duplicates are selected from the feature class. There is a field called "Spots" within this feature class that has numeric values from 0-10. Of the zip-codes that have duplicates, I want to delete the feature that has the smaller number in the "Spots" field. Furthurmore, if the spots field is the same for the duplicates, I'd like to keep only the bottom row. I hope this makes sense! Any help would be great!


texts: Install osgeo live package in existing Ubuntu?

texts:
I am using Ubuntu. In my secondary computer I use Lubuntu based OSGEO-LIVE 12.0. OSGEO-LIVE is a great software when it comes to any GIS application and reduces tedious work of individual package installation. But I want to set up all the packages from OSGEO-live into my existing Ubuntu 18.4.

Is there any way of doing this in bulk in short time?


texts: How to clip a Oracle query layer directly in a Python script?

texts:
I have a query layer directly running from ORACLE database in a MXD. I want to clip this layer  with a layer being created in another MXD in a Python based tool. My question is to do this clip analysis, do I have to make a layer of this query layer locally or I can clip it directly from the query layer? I have success doing it with MakeFeatureLayer_management but then creating a layer of 6-7 million records every time this tool runs does not appear efficient and fails the whole purpose of creating query layer. SDE was never implemented for this project and setup was done before I joined the project.
Alternatively I tried to run the model in query layer MXD but it will work stand alone and results in error with Python.


texts:
I've used a Landsat 8 image for vegetation analysis in an urban area. I've got a large maximum number of EVI value as shown in this image:

Is this normal in calculating EVI and NDVI? Because I know that value range in NDVI and EVI are between -1 and 1, although only 1 DN has max value.

And 2nd question is, what are the best vegetation indices for analyzing in urban areas with Landsat 8 images? I use ENVI apps to calculate band math.


texts: Sample for MS Visual Studio C++ using GDAL library

texts:
I need to create a sample console application on c++ that will use *.hdf file as an input and will read it.

As I got, I need a GDAL library. I've found a GDAL tutorial on official GDAL website, but it's splitted on unended parts, so I can't find a good working example of use.

Can someone provide it to me?

Also I have some troubles with the GDAL library. As I found, I need the stable library from gisinternals, yep?


texts: ExtractByMask with one mask and multiple layers ArcMap?

texts:
I have a multiple cost distance rasters and multiple feature layers. I need to find and join the feature layers (or associated shapefiles) that lie within a threshold of costs on the cost rasters. For this, I have a loop which goes:

def joinfun(row, start, end):
    FID=str(row.getValue("FID"))
    where = '"FID" =' + FID
    out_fl='out_fl{0}'.format(FID)
    for FID in range (start,finish):
        arcpy.MakeFeatureLayer_management(fc,out_fl,where) #i want acodist_{0} to cycle through all the codist layers, while out_fl stays the same
        if arcpy.GetRasterProperties_management((arcpy.sa.ExtractByMask('codist_{0}'.format(range(start,finish)), out_fl)), "MINIMUM") <= threshold:
            arcpy.SpatialJoin_analysis('out_fl{0}'.format(), out_fl, '"joined_{0}".format(FID) "and" out_fl') ##the format() should be the number of the codist_{0} which passed the condition, so that the out_fl is the one which is considered joinable
            arcpy.Delete_management('out_fl{0}'.format(FID))
            arcpy.Delete_management('out_fl')

the issue here is the if statement. I need to find a way to see if the out_fl lies over an area of cost which is <=threshold, no matter the cost layer. Therefore, in the following statement, something needs to be done about the codist_{0}.format(range()) part:

if arcpy.GetRasterProperties_management((arcpy.sa.ExtractByMask('codist_{0}'.format(range(start,finish)), out_fl)), "MINIMUM") <= threshold:
                arcpy.SpatialJoin_analysis('out_fl{0}'.format(), out_fl, '"joined_{0}".format(FID) "and" out_fl')

I need it to use the all codist_{0} rasters I have (from 4910 of them, from 0=0 to 0=4909) and check each for each out_fl (also 4910 of them in total). Does anybody have experience in or advice on this matter?

the SpatialJoin_analysis of the statement is the most important, because in the end, it is the goal of this process. I need to join all out_fl{0} that fall within the cost distance threshold of the next. the cost distance rasters are the accumulative cost at each cell leaving from out_fl{0} outwards.


texts:
Previously I've asked an question and Micha gave me a helpful advice:How to find the watershed area for each point in the river network using DEM？

After getting the useful DEM file, I began to setup GRASS and use it to solve my problem, I downloaded it from its official site:Index of /grass64/binary/mswindows/native, choosing this one:

  WinGRASS-6.4.2-2-Setup.exe

to download.

I have downloaded the main program and two datasets, and then launch GRASS after fininshing installation:

I create an new folder and put downloaded datasets in it, in my case,

  C:\grassdata

But as I choose this folder as my GIS data directory, the window soon vanished and nothing came up.

I tried many times but the result was just all the same, so I want to setup GRASS correctly, and then use these two useful modules:

r.watershed elev=DEM thresh= drain=flow_dir_grid
r.water.outlet drain=flow_dir_grid basin=basin1 easting=X1 northing=Y1

to solve my problem.

How to setup GRASS and could someone provide step by step guidance?Thanks.


texts: Extract by Mask tool not clipping to mask extent?

texts:
I am using ArcMap 10.4.1.

I am using the Extract by Mask tool to clip a raster using a polygon shapefile. When the tool finishes running, the output raster does not fill the entire mask extent (see attached image). I have made sure that my raster and shapefile are using the same coordinate system.

Any idea what's going on?


texts: Understanding the Map Tile creation properties

texts:
I am trying to create a tile base map.
I have an image that I want to create the tile from and I am using GDAL.

I managed to create the tiles but when I view the result it seems that the tiled map is:

It is flipped
The proportions are not right

This is the map overlay image:

This is the result of the tiling:

I checked many times that the coordinates are good by comparing the image corners on an orthophoto.

An other thing is that in google earth with the same coordinates the map image fits precisely to google map.

EDITED

  gdalinfo map_overlay.png

with this output:

  Driver: PNG/Portable Network Graphics Files: map_overlay.png Size is
  3527, 2494 Coordinate System is `' Image Structure Metadata:
  INTERLEAVE=PIXEL Corner Coordinates: Upper Left  (    0.0,    0.0)
  Lower Left  (    0.0, 2494.0) Upper Right ( 3527.0,    0.0) Lower
  Right ( 3527.0, 2494.0) Center      ( 1763.5, 1247.0) Band 1
  Block=3527x1 Type=Byte, ColorInterp=Red   Mask Flags: PER_DATASET
  ALPHA  Band 2 Block=3527x1 Type=Byte, ColorInterp=Green   Mask Flags:
  PER_DATASET ALPHA  Band 3 Block=3527x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA  Band 4 Block=3527x1 Type=Byte,
  ColorInterp=Alpha

Then translate:

  gdal_translate -of VRT -a_srs EPSG:4326 -gcp 0 0 31.7431761644 35.1680410195 -gcp 3527 0 31.7493769674 35.1784535489 -gcp 3527 2492 31.7431011291 35.1784951643 -gcp 0 2494 31.7431761644 35.1680410195 map_overlay.png map_overlay.vrt

with this output:

  Input file size is 3527, 2494

and lastly:

  gdal2tiles.py -z 14-21  map_overlay.vrt

The overlay is a simple .png image without any special data.

I will appreciate any help.

Shani


texts: Identify duplicate attributes and keep value of earliest record

texts:
I have rows in a table that have duplicate ID values in one field. I would like to identify these duplicates, keep the ID value of the earliest entry, and null out the ID field for the other duplicates. Since there is no date of entry field, I am relying on the ObjectID to identify the earliest entry (the lowest number). I am using a file geodatabase table, and operating in Arc 10.3.

For example:

999999 repeats twice, therefore I want to keep the ID for ObjectID 1, but null out the fm_uniqid for ObjectID 2. The same goes for 999998, which is repeated 3 times. I did the following manually.

from this:

to this:

The following code (source), helps me identify the duplicates, but now I'm wondering about my next steps.

import arcpy

fc = r"C:\Users\YCZ\testgdb.gdb/test_table"
checkfield = "fm_uniqid"

def checkDuplicate(fc, checkfield):
    vals = [row[0] for row in arcpy.da.SearchCursor(fc, checkfield)]
    dct = {x:vals.count(x) for x in vals}
    for k, v in dct.iteritems():
        if v > 1:
            print k, v

if __name__ == '__main__':
    checkDuplicate(fc, checkfield)

I am looking for advice on how to continue before I attempt any code. I believe the best way may be to put these identified duplicates into a list, then iterate through each value individually: Select by attribute for duplicated fm_uniqueIDs, identify the lowest objectID value, and set the fm_uniqid of all rows with an objectID greater than the lowest value, to null.

This seems quite tedious for a large dataset (40,000+ records)--are there recommendations of other approaches?


texts: Natural Earth: How to get the position of each country?

texts:
I've just split the 1:10m Natural Earth countries into individual JSON files which I plan to render and partially stitch back together using RaphaelJS. The 1:10m map is very big (10mb or so) so loading it all in one go in the browser is not really feasible.

But how do I get the position of each country? Opening the 1:10m map in ArcMap's Attribute Table shows no signs of any latitude, longitude, X or Y or whatever.

Maybe this info needs to be extracted from an exported JSON of the full map using some JS virtual canvas / positioning system? I.e. for each of the individual countries that need to be stitched together, some server-side script would query the full JSON to get their paths?

Thanks!


texts: Does ArcPy have spatial search function for geometry?

texts:
My problem:
Using ArcPy I want to loop over buffers and select geometry_features inside each buffer and do something (update) only the objects found within that specific buffer.
The code below explains what I want to do (more or less):

def _update_connections_inside_buffers(self):
    buffers = arcpy.SearchCursor(self.__buffer_class_name)
    in_layer = "connections"
    for i_buffer in buffers:
        shape = i_buffer.shape
        # can not use a geometry to do a selection.. very inconvenient!!
        connections = arcpy.SelectLayerByLocation_management(in_layer, "WITHIN", shape)
        self._update_connections(connections)

However: this will not work because SelectLayerByLocation_management() does not accept a geometry, "shape" as an argument but expects a feature_class_name (layer name). Is there a arcpy method that can do a search using a spatial predicate. I could not find how to do this from the ESRI manual.


texts: PostGIS: How to find Nearest Neighbour but only within a direction range?

texts:
I am attempting to find the nearest neighbour to a geometry in a certain direction but am becoming stuck at the planning stage.
I understand the basic tools that I need but am failing to find anything that can limit on a direction range.
I am finding the nearest neighbour and getting the coordinates of the two closest points as so:
SELECT
    tb1.degrees
    --taking ST_StartPoint and ST_EndPoint to get my closest points
    ,ST_ShortestLine(
        tb1.geom
        ,(
            SELECT geom
            FROM table2 tb2
            ORDER BY tb1.geom <-> tb2.geom ASC
            LIMIT 1
        )
    )
FROM table1 tb1;

But I do not want to find the nearest neighbour to the 'left' of this geometry, only right. When I say left, I actually mean (degrees-90) of the bearing for tb1.geom.
So how do I find the nearest neighbour but only in directions n°+180? Or, is this possible?


texts: Where do I find the ogr2ogr command for Mac OS X Lion?

texts: I am looking for the equivalent of FWTools for Mac OS X, Lion.  Is there anything I can easily install and get that command?


texts: Error when adding WMS layer to GeoNetwork Map

texts:
I am setting up local stack of pgadmin4 (v2), geoserver (2.12.1) and geonetwork (3.4)(on tomcat 9).

I've been testing to see if the layers i've uploaded to pgadmin and added to geoserver work, and they do. I can easily add them as a wms service with QGIS.

However, I am trying to add them also to map of geonetwork, and it doesn't seem to work. No matter which layer I use, it won't work. Even other WMS services that I know work.

I get this error:

  Something went wrong while loading tile 'http://localhost:8080/geoserve ...' for layer 'river_ciliwung[3857]'.

Very strange, because geonetork is able to do a WMSrequest  GetCapabilities no problem.

See picture, it retrieves all the layers:

furthermore, I used chrome (f12) to inspect the requests, and this is also very strange. GeoNetwork does a bunch of requests and they all get canceled.

How do I add a wms map layer to the geonetwork map?


texts:
I'm new in QGIS and I’m having some troubles. I have produced a surface plot in python and I would like to import it in QGIS but I don't know how do that, can anyone help me?

P.S. I also have the long lat coordinates of the plot for the geo referencing


texts:
I'm trying to implement a standard table partitioning scheme in postgres but I end up with this error.

Using postgres 9.4 and postgis 2.1 on centos 7

Tables:

create table experiments (
    id serial,
    value int
);

create table nada_neurons (
    id serial,
    experiment_id int,
    geometry geometry(MultiPointZ)
);

functions:

CREATE OR REPLACE FUNCTION
public.create_partitions()
RETURNS TRIGGER AS
$BODY$
BEGIN

EXECUTE 'CREATE TABLE public.nada_neurons_' || LASTVAL() || ' (CHECK ( experiment_id = ' || LASTVAL() || ')) INHERITS (public.nada_neurons)';
EXECUTE 'CREATE INDEX neurons_experiment_id_' || LASTVAL() || ' on public.nada_neurons_' || LASTVAL() || ' USING btree(experiment_id)';
EXECUTE 'CREATE INDEX neurons_geometry_' || LASTVAL() || ' on public.nada_neurons_' || LASTVAL() || ' USING gist(geometry gist_geometry_ops_nd)';

RETURN NULL;
END;
$BODY$
LANGUAGE plpgsql;

CREATE TRIGGER add_experiment_trigger
AFTER INSERT ON public.experiments
FOR EACH ROW EXECUTE PROCEDURE public.create_partitions();

CREATE OR REPLACE FUNCTION neurons_insert_function()
RETURNS TRIGGER AS $$
BEGIN
    EXECUTE 'INSERT INTO nada_neurons_' || NEW.experiment_id || ' VALUES ' || NEW.* || '';
    RETURN NULL;
END;
$$
LANGUAGE plpgsql;

CREATE TRIGGER insert_neurons_trigger
    BEFORE INSERT ON nada_neurons
    FOR EACH ROW EXECUTE PROCEDURE neurons_insert_function();

The insert:

postgres=# insert into nada_neurons(experiment_id, geometry) values(5, 'MULTIPOINT(0 0 0, 5 5 5)');
ERROR:  column "geometry" is of type geometry but expression is of type numeric
LINE 1: INSERT INTO nada_neurons_5 VALUES (8,5,010400008002000000010...
                                               ^
HINT:  You will need to rewrite or cast the expression.
QUERY:  INSERT INTO nada_neurons_5 VALUES (8,5,01040000800200000001010000800000000000000000000000000000000000000000000000000101000080000000000000144000000000000014400000000000001440)
CONTEXT:  PL/pgSQL function neurons_insert_function() line 3 at EXECUTE statement

EDIT: looking at kttii's response below, this blog post by Paul Ramsey and this stackoverflow question I believe the core issue is that the EXECUTE statement causes the PostGIS constructor to be evaluated.  This would be fine if it then enclosed the result in quotes, but since it doesn't I need to find either:
- A method of getting the EXECUTE to ignore a section of the string.
- A different technique for managing the insert trigger.

EDIT 2:
Here is how I solved this issue, and included the correct usage of the format() and USING features

CREATE OR REPLACE FUNCTION neurons_insert_function()
RETURNS TRIGGER AS $$
DECLARE
    sql_stmt varchar(200);
    r nada_neuron%ROWTYPE;
BEGIN
    EXECUTE format('INSERT INTO nada_neuron_' || NEW.experiment_id
    || ' (name, cell_type, geometry, keypoint, experiment_id, layer_id) '
    || ' SELECT $1, $2, $3, $4, $5, $6 ' )
    USING NEW.name, NEW.cell_type, NEW.geometry::text, NEW.keypoint::text, NEW.experiment_id, NEW.layer_id;
    RETURN NEW;
END;
$$
LANGUAGE plpgsql;

DROP TRIGGER IF EXISTS insert_neurons_trigger ON public.nada_neuron;
CREATE TRIGGER insert_neurons_trigger
    INSTEAD OF INSERT ON nada_neuron
    FOR EACH ROW EXECUTE PROCEDURE neurons_insert_function();


texts: How to add information to an attribute table?

texts:
I'm working with qgis 1.8.0 and did the following:

1.-Divided my shapefile into 830 features (small rectangles).
2.- Added to my shapefile several points containing weather information using the "add delimited text layer" tool.
3.- Did some weather data interpolations.
4.- Used the "point sampling tool" to collect the interpolated information for each and every one of the 830 features of the shapefile.

QUESTION: is there any way to add/import the information obtained with the point sampling tool to each of the 830 features of my original shapefile? I tried to modify the attribute table but couldn't find any way to paste information to a new column.. any ideas?

Thanks!


texts: Find COUNT of multiple tables that ST_INTERSECT a main table

texts:
I have 3 tables. 1 polygon main table polytable (400,000 polys with 1000 unique names in the polyname column), and 2 point tables pointtable1 (millions of points) and pointtable2 (1000s of points). I'd like to return the distinct list of features from polytable with the COUNT of how many points intersect from pointtable1 and also how many points intersect from pointtable2.
I want something like this:
polyname    pointtable1count   pointtable2count
red area    1                  2
blue area   0                  5
green area  9                  0
white area  0                  0

Notice how I'd like the polygon white area to still return even though there are no points that exist in it.
I've tried this query below but it seems to take too long. In a separate query when only using one of the point tables, I'm also having trouble returning polygons with 0 points within them (ie. white area would not be returned).
SELECT DISTINCT poly.polyname, COUNT(pnt1.*) as pointtable1count, COUNT(pnt2.*) as pointtable2count
FROM schema.polytable poly,
schema.pointtable1 pnt1,
schema.pointtable2 pnt2
WHERE ST_INTERSECTS(poly.geom, pnt1.geom) OR ST_INTERSECTS(poly.geom, pnt2.geom)
GROUP BY poly.polyname
ORDER BY poly.polyname ASC;


texts:
I like to know if there is a solution to the following problem. I currently have a simple ModelBuilder workflow in which includes a list of shapefiles being copied. Within the workflow, once copied, I  would like to batch delete fields from the copied files, but for some reason I am unable to perform such task. When I open up the delete field tool, I am only able to select the first row. The rest are all greyed out (disabled)

I have attached two images as reference.


texts: Multipolygon postgis layer exposed via geoserver is empty at certain zoom levels

texts:
I have the following issue:

There is a postgis table that contains a geometry column of type multipolygon (epsg is 4326, if it's relevant).

This layer is exposed via geoserver. When I preview this layer with openlayers, everyting is fine starting at scale of 1 : 35M.
However, when I zoom out, objects just dissapear. I had the style set to Point, so polygons being too small should not be the issue.

Currently I have a FeatureTypeStyle which uses PointSymbolizer with <MinScaleDenominator>10000000</MinScaleDenominator> and PolygonSymbolizer with <MaxScaleDenominator>10000000</MaxScaleDenominator>, and I've tried using only PointSymbolizer with the same result.

Geoserver version is 2.8.0.

What might be the problem here?


texts:
I need to get a reasonable Sentinel 2 Surface reflectance using Google Earth Engine, I tried but the values are between 10000 - 0

var ReflectTimeSeries = ui.Chart.image.seriesByRegion({
  imageCollection: temps2013,
  regions: westernRegions,
  reducer: ee.Reducer.mean(),
  scale: 1,
  xProperty: 'system:time_start',
  seriesProperty: 'label'
});

ReflectTimeSeries.setChartType('LineChart');
ReflectTimeSeries.setOptions({
  title: 'Sentinel 2 band values at three points near Potsdam City',
  hAxis: {
    title: 'Time Series'
  },
  vAxis: {
    title: 'Reflectance (intensity)',
    Divide : 10000
  },
  lineWidth: 1,
  pointSize: 4,
  series: {
    0: {color: COLOR.PARK},
    1: {color: COLOR.FARM},
    2: {color: COLOR.URBAN}
  }
});


texts: How to query exactly one particular layer with "Identify feature" tool in QGIS?

texts: Basic question: is there a way to unable attributes consulting for a given layer? I have several layers I need to visualize at the same time, but I would like to specify that I want to question one specific layer when I click the canvas using « Identify feature » option.


texts: Merge adjacent polygons and return all (adjacent and not adjacent) as Multipolygon?

texts:
I have a postgres database with administrative boundaries and their geometry.
Each of those boundaries have an ident number.

What I want to achieve:

First I want to select all rows which are starting with a particular zip code.
For example:

SELECT * FROM "post" WHERE "post"."ident" LIKE '101%'

All of those Geometries which are adjacent should be merged and those which are not should also be returned.

I have tried those statements:

This one I have from this answer:

Joining lots of small polygons to form larger polygon using PostGIS?

with t as (
select (st_dump(geom)).geom
from "post"
where "post"."ident" LIKE '593%'
) select ST_UNION(geom) from t;

But this statement returns only those polygons which are adjacent but others are not included

with t as (
select (st_dump(geom)).geom
from "post"
where "post"."ident" LIKE '593%'
) select ((ST_Dump(ST_UNION(geom)))) from t;

I get 2 rows, and it seems that this could be what I want, but I cannot get the result als ST_AsText.

Is there a statement which returns all polygons merged where necessary and return it ideally as geoJSON or maybe as text?


texts:
I want to project a raster from Sinusoidal to UTM. I used the following command in Python:

gdalwarp -overwrite -t_srs EPSG:32643 -of GTiff D:/IMAGE.tif D:/IMAGE2.tif

I'm getting the following error:

SyntaxError: invalid syntax

Anyone who can help?


texts: Openlayers 3 WMS custom projection wrong extent

texts:
I'm using OpenLayers 3 to view a WMS layer from a server in EPSG:3006, see the code:

Proj4js:

Proj4js.defs["EPSG:3006"] = "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs";

ol.js:

var map;
var projection = ol.proj.get('EPSG:3006');
function init(){
    map = new ol.Map({
        target:'map',
        renderer:'canvas',
        view: new ol.View({
            projection: projection,
            center:[617000,6728500],
            zoom:3
        })
    });

var topolayer = new ol.layer.Image({
    source: new ol.source.ImageWMS({
        preload: Infinity,
        url: 'URLgoesHere',
        params:{
            'LAYERS':"layernameGoesHere"
        }
    })
    });

map.addLayer(topolayer);

The result is a map centered over Sweden/Denmark, but it should result in a map centered over middle part of Sweden. I also added a button which calls a function to print the current extent of the ol.View (got help from a friend with this function, believe its from this post https://stackoverflow.com/questions/28166471/openlayer3-how-to-get-coordinates-of-viewport):

function extentCoords() {
var x=map.getView().calculateExtent(map.getSize());
alert(x);
}

For now it is in an alert() call just for experimenting and troubleshooting. When this function is called it returns values that are far away or totally wrong from the actual boundaries of the projection, for example when zoomed in over a small city in the mid part of Sweden it returns the following extent:

1905097,691, 8549623.590, 1910591.602, 8554759.203

The result should have X-values in the six figure range, rather than the seven figure range like now. The returned Y-values are in the seven figure range as expected, but not within the right values.

Also when changing the Center value, the map doesn't react to the new values.

The reason to why I need to extract these coordinates is because they are used to call a WCS through another function and download data within the same area/extent.

Is there something that I'm missing when using this custom projection, something needed to be transformed or maybe added in the ol.View part? Something wrong with the extent-function?

Sidenote: I'm kind of new to using OpenLayers, especially with WMS layers and custom projections. Also, I've got no control over the server as it is not mine.


texts: Plotting of non-spatial data with the same geometry

texts:
In order to create my tmap, I merged my spatial data (which includes State ZIP Codes, Cities and ZIP Codes Polygon Geometry) with my non-spatial data (Individual Personal Info, Individual's ZIP Code, and other variables). In total I have 1500 individuals so as you can imagine, I have more than one individual per ZIP Code.
Therefore, by the time I plot the tmap for a specific variable (e.g weight) I am not sure whether it plots only one individual per polygon or plots a consensus (e.g if most individuals with the same polygon geometry weight 132 lb, the polygon color is one that reflects 132 lb).
How can I check that and what can I do in order to represent all the individuals per polygon?
# Open libraries
knitr::opts_chunk$set(echo = TRUE, warning=FALSE)
suppressMessages(library(tidyverse))  # Modern data science workflow
suppressMessages(library(sf)) # Simple features for R
suppressMessages(library(dplyr))
suppressMessages(library(tmap))       # Thematic Maps
suppressMessages(library(tmaptools))  # Thematic Maps Tools
suppressMessages(library(RColorBrewer)) # ColorBrewer Palettes
suppressMessages(library(leaflet))    # Interactive web maps
suppressMessages(library(rgdal))      # Bindings for the Geospatial Data Abstraction Library
suppressMessages(library(rgeos))      # Interface to Geometry Engine - Open Source
suppressMessages(library(rio)) # for reading in data
suppressMessages(library(htmlwidgets))
suppressMessages(library(xlsx))  # trabajar con excel
suppressMessages(library(gifski))

# Read shapefile (includes ZIP Code and Geometry = Polygons)
cp_cdmx <- read_sf("~/CP_CdMx")

# Read file with districts and ZIP Code
deleg_cdmx <- read.csv("~/cp_deleg_cdmx_uniq.csv", header = FALSE)

# Exchange names from district file
colnames(deleg_cdmx) <- c("d_cp", "delegacion")
deleg_cdmx$d_cp <- as.character(deleg_cdmx$d_cp)

# Join geometry, ZIP code and district. This join is by ZIP code
cp_deleg_cdmx <- cp_cdmx %>%
  left_join(deleg_cdmx, by = 'd_cp')

# Import non-spatial data
basal_db <- rio::import("~/subset_basal.xlsx")

# Rename ZIP code column for common name with the other file
names(basal_db)[names(basal_db) == "cp_actual"] <- "d_cp"

# Process data
cp_deleg_cdmx$d_cp <- as.numeric(cp_deleg_cdmx$d_cp)
basal_db$d_cp <- as.numeric(basal_db$d_cp) # NAs are inserted where no ZIP code exist

# Merge spatial and non-spatial data
# The merge is by ZIP code
basal_map_cdmx <- merge(cp_deleg_cdmx, basal_db, by = "d_cp")

# Keep the data as sf class, so we will not lose the coordinate system
basal_cdmx <- st_as_sf(basal_map_cdmx)
st_crs(cp_deleg_cdmx)

# Plot Map
tmap_mode("plot")

# First layer = all CDMX ZIP codes
map <- tm_shape(cp_cdmx) +
  tm_borders(col = "#737373", alpha = 0.3) +

# Second layer
tm_shape(basal_cdmx) +
  tm_polygons(col = "actfisicatrabajo", title = "Post-Job Workout",
              palette = 'Blues',
              legend.hist = TRUE, legend.z = 0, legend.hist.z = 0.5, legend.shape.z = 0, border.alpha = 0.1) +
  tm_facets(by = "sexo", as.layers = FALSE, free.scales = TRUE, showNA = FALSE) +
  tm_legend(text.size=1, title.size=1.2, position = c("left","top"), bg.color = "white", bg.alpha= 0.1, frame= FALSE,
            height= 0.3, hist.width= 0.2, hist.height= 0.3, hist.bg.alpha= 0.5) +
  tm_scale_bar(color.dark = "gray60", position = c("right", "bottom")) +
  tm_view(colorNA = NULL, view.legend.position = c("right","top"), set.zoom.limits = c(5, 15)) +
  tm_layout(frame = "gray60", frame.lwd = 0, inner.margins = 0, outer.margins = 0, asp = NA, legend.hist.width = 0.15, legend.hist.height = 0.15,
            legend.height = 0.4, legend.width = 0.4, legend.outside.size = 0.2, panel.labels = c('Hombres','Mujeres'),
            panel.label.size = 1, panel.label.bg.color = "azure3")

All the files I use are in this Drive link


texts: Why does the output of gdal.RasterizeLayer differ from gdal_rasterize?

texts:
Problem

I've recently noticed that the output from the gdal.RasterizeLayer function differs from the result I'm getting using the gdal_rasterize utility.

I have geotiff of a SAR product that has been resampled to 500m and my aim is to create a corresponding land mask geotiff. The extent of the data is centered around Florence, Italy and covers the width of Italy with some sea to the east and west.

I'm using the GSHHS L1 vector dataset for the masking and I've written a Python function to perform the operation:

import os
import numpy as np
from osgeo import gdal, ogr, osr

# Source raster details
# src_rast_gt = (597878.6395908913, 500.0, 0.0, 4952884.831090534, 0.0, -500.0)
# src_rast_xsize = 560
# src_rast_ysize = 364
# src_dtype = gdal.GDT_Float32
# src_epsg = 32632

def create_land_mask(src_rast,land_vector='GSHHS_f_L1.shp'):
    # Get raster metadata
    rast_ds = gdal.Open(src_rast)
    if rast_ds is None:
        sys.exit('Failed to open src_file "{}"'.format(src_file))
    band = rast_ds.GetRasterBand(1)
    geotrans = rast_ds.GetGeoTransform()
    xsize, ysize = band.XSize, band.YSize

    mask_ds = ogr.Open(land_vector)
    if rast_ds is None:
        sys.exit('Failed to open land_vector "{}"'.format(land_vector))
    mask_lyr = mask_ds.GetLayer()

    # Create SpatialReference objects to hold raster and vector srs
    rast_srs = osr.SpatialReference(wkt=rast_ds.GetProjection())
    mask_srs = mask_lyr.GetSpatialRef()

    # Get raster corner points
    rast_tl = (geotrans[0] + 0*geotrans[1] + 0*geotrans[2],geotrans[3] + 0*geotrans[4] + 0*geotrans[5])
    rast_tr = (geotrans[0] + xsize*geotrans[1] + 0*geotrans[2],geotrans[3] + xsize*geotrans[4] + 0*geotrans[5])
    rast_br = (geotrans[0] + xsize*geotrans[1] + ysize*geotrans[2],geotrans[3] + xsize*geotrans[4] + ysize*geotrans[5])
    rast_bl = (geotrans[0] + 0*geotrans[1] + ysize*geotrans[2],geotrans[3] + 0*geotrans[4] + ysize*geotrans[5])
    coords = [rast_tl,rast_tr,rast_br,rast_bl,rast_tl]

    # Generate raster extent geometry
    rast_ext = ogr.Geometry(ogr.wkbPolygon)
    lin_ring = ogr.Geometry(ogr.wkbLinearRing)
    for c in coords:
        lin_ring.AddPoint_2D(*c)
    rast_ext.AddGeometry(lin_ring)
    rast_ext.AssignSpatialReference(rast_srs)

    if not rast_srs.IsSame(mask_srs):
        rast_ext.TransformTo(mask_srs)

    # Apply a spatial filter to the shapefile layer.
    mask_lyr.SetSpatialFilter(rast_ext)

    # Create raster to store mask
    drv = gdal.GetDriverByName('GTiff')
    out_path = os.path.splitext(src_rast)[0]+'_mask.tif'
    if os.path.exists(out_path):
        drv.Delete(out_path)
    out_rast = drv.Create(out_path,xsize,ysize,1,gdal.GDT_Byte)
    out_rast.SetGeoTransform(geotrans)
    out_rast.SetProjection(rast_srs.ExportToWkt())

    # Rasterize filtered layer into the mask tif
    gdal.RasterizeLayer(out_rast, [1], mask_lyr, burn_values=[1],options=['ALL_TOUCHED=TRUE'])

    # Invert with Numpy - option 'INVERT=TRUE' had no effect
    out_rast.FlushCache()
    out_band = out_rast.GetRasterBand(1)
    out_arr = out_band.ReadAsArray()
    out_arr = np.logical_not(out_arr).astype(np.uint8)
    out_band.WriteArray(out_arr)
    out_band.ComputeStatistics(False)

    out_rast = None
    mask_ds = None
    rast_ds = None

The corresponding gdal_rasterize command is below:

gdal_rasterize GSHHS_f_L1.shp gdal_rasterize_mask.tif -burn 1 -at -i -l GSHHS_f_L1 -te 10.202 42.997 13.764 44.723 -ts 364 560`

Results

The output of gdal_rasterize is as expected with 0 values over land and 1 over the sea:

However, the output of my function only masks out the coastlines:

I've run the same function on a sample of other SAR scenes over different locations and it works fine, successfully masking all pixels that fall on land. The issue seems to be with this location.

I've ran this comparison using GDAL 2.1.2 and 2.2.3 with the same results.

Question

Is there a problem with my code, or has anyone noticed similar behaviour?


texts: QGIS "Raster->Georeference" menu item not showing up

texts:
I am using QGIS 2.12.0 Lyon. I have installed a plug-in called "Georeferencer GDAL" and the version of that plug-in is 3.1.9 . [I am working on a Mac which is running OS X Yosemite 10.10.5.]

I am trying to follow tutorials on georeferencing. A critical step is "Click on Raster -> Georeference".

In my Raster menu, no item named "Georeference" appears. But the Plug-in Manager assures me that the plug-in "Georeferencer GDAL" is in fact installed. How can I get it to appear in the Raster menu, or in some other way that makes it accessible for my use?


texts: Apply new style to selected features in QGIS

texts:
Let's assume that I have selected some features in QGIS. Now all the features in that vector layer already has some style applied to it.
Now, I want to apply some new style to the selected features. Is this possible in QGIS?


texts:
What is a memory vector layer and how does it relate to the vector layer?

Background info:

I have a volunteer working with my organization, because of financial reasons we have him using QGIS.  He is relatively new to GIS so has some questions.  Being a ArcGIS user I can field most of his questions, as (it seems) the 2 programs have a lot of similarities.


texts: Georeferenced PDF is checked but doesn't georeference

texts: I have checked the georeference option for my exported PDF but when I import it into Avenza PDF Maps it says it is not georeferenced.  I have had this work for other maps but a new map I opened will not work.  My new map has an aerial and a polygon I drew and that is all.


texts: PostGIS - Why are my ST_Union results (stored in a new table) outputting as points and not polygons?

texts:
I've decided to move to PostGIS to merge polygons rather than spend hours doing in in QGIS. However, I'm hitting an issue where the results of an ST_Union query are being stored, or at least outputted as points and not merged polygons.

Here is my SQL query:

DROP TABLE opengeo.london_merged;
CREATE TABLE opengeo.london_merged (
  the_geom geometry(MultiPolygon),
  CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL)
);

INSERT INTO opengeo.london_merged
SELECT ST_Union(the_geom) AS the_geom
FROM opengeo.tq_building

This is the output from the new opengeo.london_merged table that I'm getting:

However, if I just run a straight this same SQL on the orignal opengeo.tq_building table as an SQL view in GeoServer (so not creating a new table, just outputting the results):

SELECT ST_Union(the_geom) AS the_geom
FROM opengeo.tq_building

I get the following, correct output:

There's obviously something going on with how the new table is created or how the data is stored / retrieved from it, I just don't yet know enough about PostGIS to put my finger on it.

How can I store the ST_Union output in a new table so that the output is the expected polygons?

Edit

I ran the following SQL view on the opengeo.london_merged table in GeoServer and it gives me the correct output, suggesting that the data is stored correctly but it being outputted / given to me wrong when I view the table normally. Perhaps a type or SRID issue? I don't know.

SELECT the_geom FROM opengeo.london_merged

I'm still keen to solve this.


texts: Contextily add_basemap uses wrong extent and zoom level

texts:
I'm trying to add a contextily basemap to a Matplotlib figure containing a GeoPandas data frame. When I just plot the data frame using df.plot the map extent is calculated correctly.
However, when I try adding a contextily basemap the map extent (and zoom level) is calculated wrongly and the following warning is shown:
UserWarning: The inferred zoom level of 27 is not valid for the current tile provider (valid zooms: 0 - 20).

I'm trying to execute the following code:
df = gpd.read_file('linz/StatBez_Linz_EPSG_4326.gml')
df = df.to_crs(epsg=3857)

fig = plt.figure(figsize=(16,9))
ax = plt.subplot()
ctx.add_basemap(ax = ax, source=ctx.providers.Stamen.Toner, crs=df.crs.to_string())
df.plot(color='none',edgecolor='green', ax = ax)

The output of df.tail() can be seen here:

The gml file is from data.gv.at


texts: How do I force the label to display above the lines

texts:
How do I force the label to stand out above the red lines to make the label more readable ?

Here is the sample of the image


texts:
I tried to shorten a given LineGeometry with some Python code.
layer = iface.activeLayer()
fs = layer.selectedFeatures()
for f in fs:
    fg = f.geometry()
    geomSingleType = QgsWkbTypes.isSingleType(fg.wkbType())
    if fg.type() == QgsWkbTypes.LineGeometry:
        if geomSingleType:
            x = fg.asPolyline()
            print("length: ", fg.length())
        else:
            x = fg.asMultiPolyline()
            print("MultiLine: ", x, "length: ", fg.length())
    newLine = fg.extendLine(0,(-2.5))
    print("length: ", newLine .length())

That returns a line with an equal length:
length:  11.995677929842055
length:  11.995677929842055

When I use a value > 0 the length is updated:
length:  11.995677929842055
length:  14.495677929927327

So I'm able to make a line longer, but I want to make it shorter.
Does anyone have any advice for me?
The code is part of a bigger project - so solutions with QGIS UI will not help me.


texts: How to reduce number of raster bands of a RGB picture?

texts:
For my Master Thesis in Ecology I want to calculate an NDVI Index.
The tricky part is, that I don't have sattelite pictures with bands and so on, but RGB pictures. The pictures where taken with a normal digital camera and the NIR pictures were taken by using a filter to block the light from 680 nm onwards.
The information of the red channel were obtained by extracting this channel out of a normal RGB picture, taken with the same digital camera as the NIR pictures. I stored this extracted channel as a TIF-file.

Now I loaded these two pictures in QGIS, render type "platted". Starting the raster calculator per picture 3 bands appear.

How can I reduce these raster channels in raster calculator to only one raster channel per picture?

Could I combine the information of the 3 bands to one somehow?

Could a simple addition of all three bands per picture help?


texts: How to set correct geo information on pleiades images with RPC model?

texts:
Now I have pleiades raster (img.tif) and a RPC model (rpc.XML). In the rpc.XML, there are the following lines which can be used to reset geo information.

  <RFM_Validity>
    <Direct_Model_Validity_Domain>
      <FIRST_ROW>1</FIRST_ROW>
      <FIRST_COL>1</FIRST_COL>
      <LAST_ROW>19676</LAST_ROW>
      <LAST_COL>18665</LAST_COL>
    </Direct_Model_Validity_Domain>
    <Inverse_Model_Validity_Domain>
      <FIRST_LON>7.393918390116374</FIRST_LON>
      <FIRST_LAT>46.89414037853461</FIRST_LAT>
      <LAST_LON>7.531577540047203</LAST_LON>
      <LAST_LAT>46.98982504658424</LAST_LAT>
    </Inverse_Model_Validity_Domain>
    <LONG_SCALE>0.06584113709872685</LONG_SCALE>
    <LONG_OFF>7.462641690331035</LONG_OFF>
    <LAT_SCALE>0.04690424904394419</LAT_SCALE>
    <LAT_OFF>46.94198271255942</LAT_OFF>
    <HEIGHT_SCALE>245</HEIGHT_SCALE>
    <HEIGHT_OFF>765</HEIGHT_OFF>
    <SAMP_SCALE>9332</SAMP_SCALE>
    <SAMP_OFF>9333</SAMP_OFF>
    <LINE_SCALE>9837.499999999996</LINE_SCALE>
    <LINE_OFF>9838.499999999996</LINE_OFF>
  </RFM_Validity>

And I use the following code to reset geo info of img.tif with the value below.
xmin = 7.393918390116374
ymin = 46.89414037853461
xmax = 7.531577540047203
ymax = 46.98982504658424
first_col = 1
last_col = 18665
first_row = 1
last_row = 19676

def set_geo_info(in_path, out_path, xmin, ymin, xmax, ymax,
                 first_row, first_col, last_row, last_col, *args):
  '''
  xmin, ymin, xmax, ymax are the longitude and latitude.
  '''
  print('Set Geo Information ...')
  xres = (xmax - xmin) / (last_col - first_col + 1)
  yres = (ymax - ymin) / (last_row - first_row + 1)

  geotransform = (xmin, xres, 0, ymax, 0, -yres)

  ds = gdal.Open(in_path)
  band = ds.GetRasterBand(1)
  arr = band.ReadAsArray()
  [y, x] = arr.shape

  driver = gdal.GetDriverByName("GTiff")
  outdata = driver.Create(out_path, x, y, 1, gdal.GDT_UInt16)
  outdata.SetGeoTransform(geotransform)
  outdata.GetRasterBand(1).WriteArray(arr)
  outdata.FlushCache()

  outdata = None
  arr = None
  band = None
  ds = None

  print("Done.")

The result is kind of correct, but not entirely.

As you can see, the red line is a shapefile (EPSG21781), while the railway in the result is not entirely overlap with shapefile. The distance of shifting is not consistent for this whole image. The displacement is larger near the edge of this image.

My questions:

Am I heading the right direction to reset geo information?
If 1. is correct, I guess the value which I fed into the gdal is wrong. How do I find the correct value from RPC model (rpc.XML)?
If 1. and 2. is wrong, how do I do it in a right way?


texts: Characterizing territory according to its relief using ArcGIS Desktop?

texts:
I want to characterize the territory according to its relief with ArcGis 10. I'm trying to define a wide territory by its relief, but I would like to do it in an automatic or semiautomatic form with the GIS.

I read on some webs that I can do it with the topographic position index, but I don't know how. Can somebody explain or give some references to learn how I can do this analysis?


texts: How to display the geojson polygon after 180th meridian in leaflet?

texts:
Like the vector grid plugin in the leaflet:
http://leaflet.github.io/Leaflet.VectorGrid/demo-geojson.html
see the above example at zoom level 1 or 2

Is there any way to display geojson, not as a tile layer?


texts:
I am using Google Earth Engine (JavaScript) to calculate brightness temperature by converting at-satellite radiance of Landsat 8. In literature the formula
T = K2 / (ln((K1/L)+1) -273.15 (T: brightness temperature; K1,K2: calibration constants 1 and 2; L: at-satellite radiance)
is used.
My code:
var bt_calc5 = function(image) {
  var bt_calc = image.expression(
    'K2 / (Math.log(K1 / L + 1) -273.15)', {
      'L' : image.select('TOAThermalRadiance')
    });
  return image.addBands(bt_calc);
};

It should add a brightness temperature band to all images of my collection, but results in Expression parse error which is why I think the syntax is not correct. I took the Math.log expression for logarithmus naturalis from the JavaScript documentation webpage. How do you formulate logarithms in Earth Engine?


texts:
I recently upgraded to QGIS 3.2 and followed the install instructions carefully.

Under the raster and vector menus I have only 1 tool for vector and only 3 for raster. No analysis, geometry, geoprocessing tools etc...

when I use the main search bar to look for any of the tools I would have used in QGIS 2 (polygonize, dissolve, clip) I get no results.

python and GDAL were install as per instructed in the QGIS 3.2 README

Am I missing something?


texts: clip one layer with the current map layers

texts:
I work in a little companie where the main goal of employees in GIS is to have a quick view of the issues on a study area.

It would be very useful to have a tool wich clip all the layers of the current map extent with the study area.

I made a quick search in ogr tools but I didn't find the appropriate one. Plus, I think that ogr can't work automatically with the current layers.

So I think the solution is Pyqgis but I have no skills in Python.

Is there an existing tool ? Any suggestions to find the simpliest way to do that ?

Thanks


texts: get valules for all fields from selected features

texts:
The following code works for retrieving values within a particular field(Name1), but I want to retrieve mulitple fields.  Am I to create multiple variables for each field?

IMxDocument pMxDocument = ArcMap.Application.Document as IMxDocument;
        IMap pMap = pMxDocument.FocusMap;
        ILayer pLayer = pMap.get_Layer(0);
        IFeatureLayer pFeatureLayer = (IFeatureLayer)pLayer;
        IFeatureClass pFeatureClass = pFeatureLayer.FeatureClass;

        ITable pTable = (ITable)pFeatureClass;
        ICursor pCursor = pTable.Search(null, false);
        int fldIndex = pFeatureClass.Fields.FindField("Name1");
        IRow pRow = pCursor.NextRow();
        while (pRow != null)
        {
            MessageBox.Show(pRow.get_Value(fldIndex).ToString());
            pRow = pCursor.NextRow();

        }

Edit:  What I'm trying to do is get the selected features of a feature class, get the values for fields like Name, address, city, etc and then pass those values onto a document.  I thought I would give some background to what I am doing since I may be on the wrong track to start off with.


texts:
I have a leaflet map and a kml file. I load the kml file using the plugin: leaflet-omnivore. Now I want to access the <folder> <name>(which is "RegionName"). I need it because I want to link the kml with an other json file which also has RegionName as a properties.name. But when output the loaded kml file via console.log(myKML) I can only access the <name> below <placemark>

The code I load the kml with:

var myKML = omnivore.kml('KML/Testdaten_KML.kml', null, L.mapbox.featureLayer());

The kml file is structured like:

<Folder>
  <name>RegionName</name>
  <Placemark>
    <name>Form 1</name>
    <Polygon>
      <outerBoundaryIs>
        <LinearRing>
          <tessellate>1</tessellate>
          <coordinates>
            11.82,52.60
            11.86,52.58
            11.90,52.61
            11.87,52.62
            11.82,52.60
          </coordinates>
        </LinearRing>
      </outerBoundaryIs>
    </Polygon>
  </Placemark>
</Folder>


texts: Seasonal analysis in Google Earth Engine

texts:
I want to perform a seasonal analysis of the MODIS NDVI dataset in GEE. I want to analyse for every year the months from May to September and thus, I need to remove the rest of the months from the dataset for every year.
I have no idea how to do this. I have found this post but this is just for one image for a year.

Any ideas?


texts: Transforming Lat/Long NASA points into raster for ArcMap?

texts:
I have this irradiation file and I would like it to upload it in ArcMap in order to overlay it with other layers to do clips, intersections, etc?

http://eosweb.larc.nasa.gov/sse/global/text/direct_radiation


texts: Calculating suitable area based on product of 2 binary layers

texts:
In QGIS I am attempting to calculate the suitable area based on two reclassified layers (Raster tools -> Reclassify by range) that have been overlaid with SAGA -> Raster calculus -> Product. Each layer is for one species, and I am trying to find the overlap area of suitable habitat. I am using binary values for suitability: 0 for unsuitable, 1 for suitable.
So, in the reclassification, each layer has only those two values. In the product, then, suitable habitat for both species (what I'm looking for) will have a product of 1, and all other values will have a product of 0.
I am looking for a way to calculate and be able to extract only the total area of overlap where the product is 1. How can I achieve that?


texts: QGIS docutils install (path,pythonpath,pythonhome)

texts:
New to QGIS and Python. Have QGIS installed and working great on Windows, want to install and use Docutils -- Python Documentation Utilities (presumably all simple and small .py programs).

I see that QGIS has a Python Console built in. Do I just use that? What is the command line to run install.py?

QGIS obviously installed Python. However it looks like the customary path variables were not set during the QIS install, such that Python can be run from a windows command (cmd.exe) window.

That seems to be the task here, so I proceed to try to enable Python to run independent of QGIS.

First I search for python.exe. I find it at:
C:\Program Files\QGIS Lyon\bin

So I set the windows path variable to include that.

After searching (and a lot of confusion as to where the Python distribution actually is, is it standard, or is it custom, questions as to how libraries, etc are laid out, and noting its all installed "underneath" C:\Program Files\QGIS Lyon) so it appears custom(?) Will I break that by using it for "generalized" Python programs?

I try the following environmental values:

PYTHONPATH
C:\Program Files\QGIS Lyon\apps\Python27\Lib

PYTHONHOME
C:\Program Files\QGIS Lyon\bin

Have no idea if those are even close to correct.
I then run the install for docutils

C:\Users\rick\Desktop\docutils-0.12>python install.py

Program starts, but then errors per below.

Seems like a "documentation utility" should be a snap to install. Plus there probably a straightforward way to run Python programs independent of QGIS? I also note an additional install of Python with Inkscape (a vector drawing program used in the QGIS training manual.) So already two installs of Python on the system already. (Correction, three installs. FWTools has yet another Python install.)

If a user wants to run Python programs (such as docutils and other small learning programs) should they download a fresh Python "stand alone" install? (looks like on this machine that would be the 4th install?)

errors

C:\Users\rick\Desktop\docutils-0.12>python install.py

This is a quick & dirty installation shortcut. It is equivalent to the
command::

python setup.py install

However, the shortcut lacks error checking and command-line option
processing.  If you need any kind of customization or help, please use
one of::

python setup.py install --help
python setup.py --help

Traceback (most recent call last):
  File "install.py", line 27, in
    dist.run_commands()
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dist.py", line 95
3, in run_commands
    self.run_command(cmd)
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dist.py", line 97
2, in run_command
    cmd_obj.run()
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\command\install.p
y", line 575, in run
    self.run_command(cmd_name)
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\cmd.py", line 326
, in run_command
    self.distribution.run_command(command)
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dist.py", line 97
2, in run_command
    cmd_obj.run()
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\command\install_l
ib.py", line 97, in run
    outfiles = self.install()
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\command\install_l
ib.py", line 115, in install
    outfiles = self.copy_tree(self.build_dir, self.install_dir)
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\cmd.py", line 377
, in copy_tree
    dry_run=self.dry_run)
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dir_util.py", lin
e 139, in copy_tree
    mkpath(dst, verbose=verbose)
  File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dir_util.py", lin
e 76, in mkpath
    "could not create '%s': %s" % (head, exc.args[-1]))
distutils.errors.DistutilsFileError: could not create 'C:\Program Files\QGIS Lyo
n\bin\Lib': Access is denied


texts: Downloading metro stations of Paris from OpenStreetMap?

texts:
I'm just getting acquainted with OSM (OpenStreetMap), so this might be a total newbie question, but I'm looking for OSM Point of Interest data for all metro stations in Paris. I'm sure this info should be available, as maps I've seen maps of Paris using OSM data (e.g., mapbox.com default maps) display these metro stations, but I'm not sure where to look.

After reviewing the answers provided and doing some digging myself, my conclusions are:

Using Overpass API seems to be the most user-friendly as far as I can tell (when you get the hang of the query language, that is). For example, the query below returns the actual routing (metro lines) plus the referenced metro stations in nice JSON format.

  http://overpass-api.de/api/interpreter?data=[out:json];(relation(48.7555,2.1406,48.9638,2.5571)[%22route%22=%22subway%22];%3E;);out;


texts: GPSBabel KML document has incorrect name

texts:
I'm using GPSBabel to convert GPX tracks in a directory to KML linestrings using the following:

@ECHO OFF
setlocal
for %%i in (C:\Users\...\*.gpx) do (
    gpsbabel.exe -t -i gpx -o kml,points=0 %%i %%~di%%~pikml\%%~ni.kml
)

This works great, except that all KML documents have the name "GPS device". I'm new to GPSBabel and after reading what appears to be the relevant documentation I'm not sure if it's possible to provide the output document's name. So far my best option seems to be post-processing with a find / replace operation but I'd like to avoid this.

How can I tell GPSBabel what to use as the KML document name?


texts: How can convert daily MODIS NDVI to monthly of them in Google Earth Engine?

texts:
The link of code:
https://code.earthengine.google.com/71d9f6ea42ee55b4e18d2adf22d9e678

As a result, I have to have 96 bands in the stack image so that the ID of each band returns the date associated with it, for example, the first band return 198801_198801 id . In fact, I want the monthly images of the winter in the stacked image.How to modify code?

code:

var monthly_map = function(collection, start, count, interval, units){

  var sequence = ee.List.sequence(0,ee.Number(count).subtract(1));
  var original_date = ee.Date(start);

  return ee.ImageCollection(sequence.map(function(i){

    var start_date = original_date.advance(ee.Number(interval).multiply(i),units);
    var end_date = original_date.advance(ee.Number(interval).multiply(ee.Number(i).add(1)),units);

    return collection.filterDate(start_date,end_date).mean()
    .set('system:time_start',start_date.millis())
    .set('system:time_end',end_date.millis());
    }))};

 var ndvii = ee.ImageCollection('NOAA/CDR/AVHRR/NDVI/V5')
.filterDate('1988-01-01','2020-01-01')
.filter(ee.Filter.calendarRange(1,3,'month'))
.filterBounds(table)
.map(function(img){
  //var id = img.id();
  var band = img.select('NDVI');
  var clip = band.clip(table);
  return clip;
  //.copyProperties(img,['system:time_start','system:time_end']);
});
print(ndvii);

var monthly_pro = monthly_map(ndvii,'1988-01-01',384, 30, 'days')
print(monthly_pro);
var monthlyband = monthly_pro.map(function(img1){
  var id = img1.id();
 var clippp = img1.clip(table);
  return clippp.rename(id)
})
.toBands();

print(monthlyband);

Export.image.toDrive({
  image: monthlyband,
  description: '1999',
  scale: 10000,
  crs: 'EPSG:4326',
  maxPixels:1e13,
  region: table});


texts:
I have the following problem: by default, you can simply change the opacity of layer by writing for example:

layer.setOpacity(0.3);

But when I add it to my map SelectFeature controls, my layer is put to a thing like OpenLayers.Vector.RootContainer. From this moment I can't set opacity of my layer.
I tried lots of ways but none of them work:

myVectorRootContainer.layers[0].setOpacity(0.4);
myVectorRootContainer.setOpacity(0.4);

Is there an option for using SelectFeature controls without RootContainers? Or a way to change opacity of a container?


texts: How to avoid Esri secure layer identity manager popup?

texts: I am new to Esri Identity manager, I am creating a secure layer service using identity manager, I have to avoid the secure layer login popup.


texts: Displaying graph preview in Python Toolbox?

texts:
I need to create a toolbox in python, which will be used for making graphs.

I need the parallel preview window during the graph creation. Half of the screen will be used for settings the parameters and the other half will be a preview that would change according to changes of parametres. Something like Create graph wizard in ArcGIS Desktop (picture below).

Does someone have python code to split the screen, or some advice, or web pages?

Picture - Create graph wizard


texts: What simplification should I apply to polygons at each zoom level

texts:
I'm sending geoJSON data to a web page and it's too big causing very slow download times.  I'm going to simplify on the server side.  My question is what tolerance I should use for my simplification at a given zoom level, is there a straight mapping or does it depend on the screen size?

My first guess is that I should work out how far one pixel is on my screen and simplify to half that distance?


texts: How to export a GeoTIFF to ArcGIS Online

texts:
I would like to add TIFs to ArcGIS Online as a layer using QGIS.
An example is the ValleyNAI layer here: https://tiles.arcgis.com/tiles/U0vgiXgwLpyQDfbW/arcgis/rest/services/ValleyNAI/MapServer
As far as I can tell, they were able to add TIFs to the tile layer. However, I cannot seem to figure out how to do so using QGIS.


texts: Grabbing Twitter locations with Python: Receiving ReadTimeOutError

texts:
I am running the following code with Python to grab Twitter locations for a specific bounding box:

import json
from tweepy import Stream
from tweepy import OAuthHandler
from tweepy.streaming import StreamListener

#Enter Twitter API Key information
consumer_key = 'ivgL3kTcnYxPtp8MLjuKuAGxq'
consumer_secret = 'n5SB7r11AwHQs0IcB3pEHpnGzdtNP9edufO5oEQc0zV8pKvlRE'
access_token = '946334425155035136-OLS2WzvtYq0qSagnhdnX00P4AcgONdw'
access_secret = 'JUDwrZQNCZSIb0B1J22vPA2SVckB4YYqwzosZhlBCU52n'

file = open("C:\Python27\Output2.csv", "w")
file.write("X,Y\n")

data_list = []
count = 0

class listener(StreamListener):

    def on_data(self, data):
        global count

        #How many tweets you want to find, could change to time based
        if count >= 0:
            json_data = json.loads(data)

            coords = json_data["coordinates"]
            if coords is not None:
               print coords["coordinates"]
               lon = coords["coordinates"][0]
               lat = coords["coordinates"][1]

               data_list.append(json_data)

               file.write(str(lon) + ",")
               file.write(str(lat) + "\n")

               count += 1
            return True
        else:
            file.close()
            return False

    def on_error(self, status):
        print status

auth = OAuthHandler(consumer_key, consumer_secret)
auth.set_access_token(access_token, access_secret)
twitterStream = Stream(auth, listener())
#What you want to search for here
twitterStream.filter(locations=[10.01,46.85,13.09,49.43])

This works well and I am getting coordinates. However, everytime after a while the program stops and I receive a read time out error which looks like this:

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    import UrbanTweets
  File "UrbanTweets.py", line 52, in <module>
    twitterStream.filter(locations=[10.01,46.85,13.09,49.43])
  File "C:\Python27\lib\site-packages\tweepy\streaming.py", line 445, in filter
    self._start(async)
  File "C:\Python27\lib\site-packages\tweepy\streaming.py", line 361, in _start
    self._run()
  File "C:\Python27\lib\site-packages\tweepy\streaming.py", line 294, in _run
    raise exception
ReadTimeoutError: HTTPSConnectionPool(host='stream.twitter.com', port=443): Read timed out.

Does anyone have an idea how to solve this issue?


texts:
I'm displaying coordinates in EPSG:900913 using ol.control.MousePosition() with OSM.
I'd like to convert these on the fly in another EPSG:2154.

I call in the header :

<script src="./proj4js-1-1-0/lib/proj4js-compressed.js" type="text/javascript"></script>

Then my html :

Proj4js.defs["EPSG:2154"] = "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs";

var map = new OpenLayers.Map({
        div: "map",
        projection: "EPSG:900913",
        controls:[
                new OpenLayers.Control.LayerSwitcher({'ascending':false}),
                new OpenLayers.Control.ScaleLine({
                    'geodesic': "true"
                }),
                 new OpenLayers.Control.Navigation(),
                 new OpenLayers.Control.KeyboardDefaults(),
                 new OpenLayers.Control.MousePosition(
                     {displayProjection: new OpenLayers.Projection("EPSG:2154")}
                 )
                ]
        });


texts: How to specify imprecise join criteria in ArcGIS for Desktop?

texts:
I need to join data from an Excel Spreadsheet to a Point feature class.

Namely, matching up Addresses from one piece of data to another, and then appending the rest of the data once the addresses match up.

The "Joins & Relates" tool worked fine, however, I was unable to specify how many characters should be considered a "Match" for the join.

I would like to be able to tell the tool that if 70% of all the characters match, or whatever, that should be considered a good match.

I apologize if my phrasing is off, but I hope I was able to get the point across.

I want to be able to specify what criteria a Join should look at.

I'm using ArcMap 10.2+.


texts: Is there a way to copy the ESRI symbol files for use in QGIS. if not id there a resource to get various svg files to use in QGIS? DIfferent linetypes as well?


texts: How to set the holeModifier on OpenLayers.Control.DrawFeature to always be active?

texts:
I would like to create a button "draw a hole" taht uses by default the holeModifier functionality but the user dont need to press the cttrlBtn.

The "standard" way I guess, would be to create a way of create holes in OpenLayers would be something like this:

var drawPolygon = new OpenLayers.Control.DrawFeature(ftthLayer,
                     OpenLayers.Handler.Polygon,{
                        handlerOptions: { holeModifier: "ctrlKey"} };

By this way the user needs to press the ctrlKey to draw a hole, or else a new polygon is created.
Is it possible to make the control DrawFeature always be in "holeModifier" mode?

I have tried setting the property drawingHole= true on the handler (OpenLayers.Handler.Polygon) that's connected to DrawFeature control , but that doesn't seem to help.

I know there are different ways to solve hole creation, but the ones I have seen havn't solved so that the user see you are drawing a hole. Like the way it's working when you use the holeModifier.

Happy for any input

Br
Mike


texts: Drawing longest polyline within polygon?

texts:
I am looking for a GIS tool/script (preference on QGIS but I can also use ArcGIS/Matlab) that allows creating a polyline layer from a polygon layer where a polyline is drawn for each polygon. The polyline is drawn being the main direction of the polygon and being within the polygon.

Below examples in image:

I do not have preference for the drawing method (vertex to vertex or edge to edge, skeletons), only the result is important:

This question can be closed as I found a suitable arcgis tool: Polyline to centerline

I would have preferred a Qgis tool but it is OK for now.

Results of the tool below on real polygon layer:


texts: Custom icons for Torque Categories in CartoDB

texts:
We’re developing a project with CartoDB and we want to show the elements clustering by zone and showing the most used one (Mode).

We're doing this with Torque Categories, and it works! :)

The questions is: how can I customize the icons of the markers?

I have tried with these options, and both make the icons disappear.

marker-file on #id[value=X]

#mapId[value=1] {
   marker-file: url(http://domain.com/Image001.png);
}
marker-file on #id[columnName=”value”]

#mapId[name="TEST"] {
   marker-file: url(http://domain.com/Image001.png);
}

I've read you cannot use marker-fill with marker-file, but I've removed all marker-fill's and still don't see the images.

How can it be done?


texts: Automatically update PostGIS layer in QGIS

texts:
Is there a way to update a PostGIS layer in QGIS (every 5 minutes for example) without needing to right click on the layer and going through the Update SQL layer dialogue as shown below?

The end goal is to have a live spatial layer that can be updated elsewhere by other users and for those updates to be reflected in the QGIS project as people work in that without them having to manually update the layer themselves.

FYI I am using QGIS 3.4.1


texts: How to create coordinate system (using Lat/Long) across IR & VIS satellite images

texts:
A recent project of mine has been to gather IR and VIS data from WeatherUnderground's online site for the purposes of creating fun and informational temporal GIFs. Now, I'd like to take it a step further and attempt to determine the cloud type at a particular location based on these two data sets (according, in part, to this paper). To do this, I will need a way to accurately compare pixel data across various regional data sets. For example, here are two images from the Central-South region:

Thankfully the state borders are clearly defined, so my first naive approach would be to manually create a coordinate system based on the edges of these states so that a cursor position would produce matching longitude/latitude coordinates, for both types of data sets. However, there are multiple regions I'd like to do this for (Central, Central-East/West/South) and realize this may not be the most efficient way to do so.

How can I efficiently and accurately compare the data from both of these image sets for the purposes of cloud classification?


texts: Resizing the map canvas to exactly the same size as a layer?

texts:
I have a problem where for the past few days I've been trying to generate binary black and white images of maps with labeled buildings from a GeoJSON file all programatically with Python.

I need to dimensions to match exactly with the original image (438 x 406). What I've done so far to generate the image below is...

# Imported the image first so the geojson matches with a correct coordinate system
image_layer = iface.addRasterLayer("path\\to\\image", "image_layer")

# Resize to the pixels I need my image as
iface.mapCanvas().resize(QSize(438,406))

# Add geojson object
vector_layer = iface.addVectorLayer("path\\to\\geojson", "geojson_layer", "ogr")

#Set background to black
iface.mapCanvas().setCanvasColor(QColor.fromRgb(0,0,0))

#Then I need to do some code to delete the image_layer

#Render the image
iface.mapCanvas().saveAsImage( IMAGE_DIR + FILENAME + ".png" )

Which then generates an image like this...(or will)

But the size is wrong...you can see it when I keep the image layer (see the extra black boundary)

I need it to be the exact same dimmensions as that original image, 438 x 406, but the image it generates is 434 x 402. Even hard coding it down to 434 x 402 then it is just 4 less than that at 430 x 398.

EDIT I tried doing...

iface.mapCanvas().setExtent(image_layer.extent())
iface.mapCanvas().refresh()

But it didn't work.

EDIT2

However, if I do iface.mapCanvas().setExtent(image_layer.extent()) outside of my script (ctrl c ctrl v into the qgis console) it resizes how I want it in QGIS but then saving it still isn't correct. (see below).

EDIT3 When I do the above manually (via the toolbar save as image) it works fine to generate the image below....I just need that done programatically!


texts: Citing ArcGIS Desktop when writing papers?

texts:
When writing papers, is there a standard style that should be used for citing ArcGIS, specifically ArcMap?

In particular, how do you reference a specific tool that you used in your analysis. I'm not referring to the citation style being used specifically in the paper.

Esri mention how to cite data produced from them, but I could not find any advice regarding citing the software.


texts: Point averaging alternative to GeoServer WPS heat map transformation

texts:
I have a point dataset in which each point has a value for a "rating" field. Points distributions don't follow any pattern and can overlap at identical locations. My requirement is to generate something like a heat map showing geographical variation of "rating", where high values are "hot" and low values are "cold".

My project already uses GeoServer's Web Processing Service and heatmap rendering transformation for a different requirement. Its existing heatmaps show the density of points that match given criteria, e.g. all points where "rating" < 1. This visualisation is meaningful because heat maps show density and only meaningful points are selected, however I can't use the same approach to show the variation across "rating". I can set the heat map process's weightAttr variable to "rating", but the resulting visualisation will not differentiate between areas with fewer points and areas with low "rating" values.

What I need is a transformation very similar to heat map but that instead of counting how many points are contained within each grid cell it calculates an average value across all values in each grid cell.

Is such a transformation possible within GeoServer / WPS, and
If not, is there documentation available on how to write and test my own rendering transformation?


texts: Togging on/off multiple layers - saving multiple layer preferences (i.e. ecological layers, utility layers, economical etc.)

texts:
I am conducting a routing review (60km), and there are hundreds of layers in the model. In order to review the route I am repeatedly toggling on/off layers inspecting the model for various factors (environmental, ecological, political etc.) If I have all the relevant layers turned on the model becomes over crowded and cluttered. Hence I am performing the same function of repeatidely turning on and off layers as I transverse the route. This is fine, however it is cumbersome/slow and prone to errors, as I could else forget to toggle a layer on during an inspection of a section of route.
Is there a way I can save a number of multi-layer searches (i.e. ecological layers, utility layers, economical etc.)
I regularly have the layers within a tree saved, however if the model is up dated (which occurs quite often when new info is obtained from the client / site), these preferences seem to be forgotten.


texts:
When i'm running my MODIS Download R script on my Windows, i'll get this error :

############################
Remote connection failed! Re-try: 49
Getting file from: LPDAAC
############################

After this, the script crashes because it tried to do gdalinfo on an inexisting image.
I'm just doing getHdf function on LPDAAC servers using my EarthExplorer Username and Password. I'll get this error on 3 different Windows PC, my personnal computer, university computer and a friend pc so on three different networks. My code runs perfectly on my Mac, do you have any idea ?

MODISoptions(gdalPath = gdalPath, MODISserverOrder = "LPDAAC")

#We have to log in into lpdaac servers using NASAs ID and password
lpdaacLogin(server = "LPDAAC")

image <- getHdf(product = "MOD09GA", begin = toString(currentDate), end = toString(currentDate), extent = c(longMin, longMax, latMin, latMax))

Here is my ModisOptions() :

'MRT_HOME' not set/found! MRT is NOT enabled! See:
'https://lpdaac.usgs.gov/tools/modis_reprojection_tool'

 STORAGE:
 _______________
 localArcPath : C:/Users/Romain/AppData/Local/Temp/RtmpGUSr8w/MODIS_ARC/
 outDirPath   :
 C:/Users/Romain/AppData/Local/Temp/RtmpGUSr8w/MODIS_ARC/PROCESSED/

 DOWNLOAD:
 _______________
 MODISserverOrder : LPDAAC
 dlmethod         : auto
 stubbornness     : high
 wait             : 0.5
 quiet            : FALSE

 PROCESSING:
 _______________
 GDAL           : GDAL 2.2.4, released 2018/03/19
 MRT            : Not available. Use 'MODIS:::checkTools('MRT')' for more
 information!
 pixelSize      : asIn
 outProj        : asIn
 resamplingType : NN
 dataFormat     : GTiff

Note : i choosed LPDAAC because LAADS FTP servers are down since the 20 april 2018.


texts: How to add north arrow to the project frame in QGIS

texts: I search for a way to add a north arrow to the project frame (and not in the composer frame- that i know how to do), in the blue mark place:


texts: QGIS: evis -- trying to use url rather than a path

texts:
QGIS 3.2 Mac OS 10.13.5

I have just added some photos to my project using importphoto to generate the layer and evis to view them.

I am now trying to get my photos independent of hard coded paths on my computer to make it much easier access the project from the 'Net.

To this end I have installed apache on the database server and have it serving the photos from an indexed directory and this works fine from the browser.

I first tried putting http://localhost:1812/gis_photos in the relative section and ticked the relative box in the path section.   From the docs I believe that this should take the filename from the path and append it to the url and then retrieve the resulting url.  I get no errors and no photo, just the metadata displayed.

I then tried editing the layer attributes to put the url in the path field with the same result.

Apache logs show no access attempts which suggests that the url is not being accessed at all.

Update:  I have logged an issue.

Further update:  This is a bug and it was present in 2.18 and therefore will have to wait before being fixed.


texts: Getting user input in ArcPy using GetParameterAsText()

texts:
I am writing a script to be included in a toolbox but every time I am getting an error when I try to include the parameter using GetParameterAsText().

import numpy
from arcpy import env
from arcpy.sa import *

input_raster= arcpy.GetParameterAsText(0)
inputFeature= arcpy.GetParameterAsText(1)
clipRaster = arcpy.GetParameterAsText(2)
floodedArea = arcpy.GetParameterAsText(3)
output_folder = arcpy.GetParameterAsText(4)

and the error I get is...

Traceback (most recent call last):
  File "C:\Python27\dealingwithraster.py", line 16, in <module>
    output_folder =arcpy.GetParameterAsText(4)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\__init__.py", line 655, in GetParameterAsText
    return gp.getParameterAsText(index)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\geoprocessing\_base.py", line 226, in getParameterAsText
    self._gp.GetParameterAsText(*gp_fixargs(args, True)))
RuntimeError: Object: Error in getting parameter as text

Failed to execute (dealingWithRaster).

Help needed because I can't seem to figure out the problem.


texts: How to select all features with maximum classification value using a processing algorithm?

texts:
I am building a Processing Model working with large vector Datasets with classification values. My vector file contains a lot of features (40.000) with one column containing classification values ranking from min (e.g. 1) to max (e.g. 15). Note that the min and max values of the classification can vary from Model-run to Model-run. I need to find a way to select all the features with the highest classification value to create a Buffer around them in a next step.

I already have a python script that writes the maximum value of a selected column in a number that than can be used in Processing:

##MaxFinder=name
##Layer=vector
##Fields=Field Layer
##value_KRK_max=output number 0

layer = processing.getObject( Layer )
idx = layer.fieldNameIndex( Fields )

KRK_max = layer.maximumValue(idx)

My idea was to define a function in the "Select by expression" tool, but I cannot make it work. Any suggestions?


texts: Reduce or heighten a DTM with overlaying polygons

texts:
I am using QGIS 3.01.

I have a DTM over a historic area. There has been done archaeological surveys there, which suggest that the ground was several meteres below todays surface many hundred years ago.

I was wondering if there is a tool or a method in QGIS to reduce/increase the height of certain areas in the raster based on values from polygons?

Like in the attached picture, with a attribute for meters that the raster calculator could do computations on or something?


texts: Python: rasterio how to cut out a roof with bigger buffor?


texts: Cannot open Spatialite db from OpenStreetMap in QGIS

texts:
I am new to QGIS, OpenStreetmap and Spatialite. I'm trying to get the (OSM) map data into a (SQLite) compact form I can eventually use for realtime reverse geocoding in a vehicle. Quantum is just for exploring the data, at this point.

I've downloaded the Saskatchewan highway layer from CloudMade (http://downloads.cloudmade.com/americas/northern_america/canada/saskatchewan) and converted it to Spatialite using

spatialite_osm_raw -o saskatchewan.highway.osm -d SK_Hi.sqlite

and

spatialite_osm_net -o saskatchewan.highway.osm -d SK_hi_net.sqlite -T Roads

Both processes run to completion, although there are lots of unresolved nodes while building the network. Neither map layer can be opened in QGis. In either case I am getting a "Failure exploring tables from: Z:/GIS/OpenStreetMap/SK/sk_hi_net.sqlite. no such column: type" when I try to connect to the database.

I am able to browse the db using the Quantum DB Manager, and I see that a number of the tables are empty. Specifically "geometry_columns_field_infos" all of the views_geometry tables and virts_geometry tables.

Is this a problem with the datasource, the toolchain, or my understanding (in increasing order of likelihood)?


texts: ArcGIS Arcade Expression: Data from a Service Query

texts:
This maybe not be possible but trying anyway. I have a Map Service with geometry features published. Two of the columns are 'Zone1' and 'Zone2' with values like 'ABC' or 'DEF'. The values of 'ABC' or 'DEF' can be present in both columns in the same row.

There is also a static 'lookup table' with a column named ZONECLASS and values like 'ABC' or 'DEF'. The static table is also published as layer in the along with the layer with geometry. So far so good.

What I am trying to do is that inside ArcGIS Online's Arcade Expression, query the static table's layer and get that data displayed inside ArcGIS Online. The query works in the browser as:

http://servicepath/MapServer/1/query?where=ZONECLASS=%ABC27&outFields=DOCNAME

But how to bring value of the DOCNAME field inside ArcGIS Online? I have built an expression as above url but that's just a string. I need to get the output.

Is that possible or do I need to seek an alternative solution?


texts:
I am writing out a shp, in ESRI shapefile driver. In R, the projection of this shp is

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
+y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I wrote this shp out as below-

writeOGR(x_projected, dsn="address path","name of shp to be stored", driver="ESRI Shapefile", overwrite_layer=TRUE)

However, on reading this projected, written shp back in, I see that the datum specified by +towgs84=0,0,0,0,0,0,0 does not exist. Also, when plotting this shp in ArcMap or QGIS, the datum is known confirming that the datum has not been written out.

Why is writeOGR() not writing out the datum i.e. why is the +towgd84=0,0,0,0,0,0,0 component missing?


texts: Plot Ranges of NDVI Values With Different Colors in Google Earth Engine

texts:
I just perform a NDVI classification, and I want to know if I can plot or display the values of the NDVI in different ranges in the same image. For example, from 0.0 to 0.3 is 'blue', from 0.3 to 0.5 'green', etc.

Here is part of my code:

//cargar coleccion de landsat
///////////////////////////////////////////////////////////////
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')

//geometry with the AOI
var roi = geometry;

//filtrar la imagen de acuerdo a la cobertura de nubes en el poligono
////////////////////////////////////////////////////////////////////
var l8f = l8.filterBounds(roi);

var withCloudiness = l8f.map(function(img_zone) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(img_zone).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: roi,
    scale: 30,
  });
  return img_zone.set(cloudiness);
});

var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));

var image = ee.Image(filteredCollection.filterBounds(roi)
    .filterDate('2017-03-30', '2017-07-30')
    .sort('CLOUD_COVER')
    .first());

    var red = image.select('B4');

    //infrarrojo cercano
    var NIR = image.select('B5')

    var nvdi = (NIR.subtract(red)).divide(NIR.add(red)).rename('NDVI');

var filtered_1 = nvdi.gt(0.3); // boolean (0,1) where nvdi > 0.3
var filtered_2 = nvdi.gt(0.5); // boolean (0,1) where nvdi > 0.5
var filtered_3 = nvdi.gt(0.7); // boolean (0,1) where nvdi > 0.7
var filtered_4 = nvdi.gt(0.8); // boolean (0,1) where nvdi > 0.8

var multiplied_1 = nvdi.multiply(filtered_1); // multiply to set <= 0.3 to 0
var multiplied_2 = nvdi.multiply(filtered_2); // multiply to set <= 0.5 to 0
var multiplied_3 = nvdi.multiply(filtered_3); // multiply to set <= 0.7 to 0
var multiplied_4 = nvdi.multiply(filtered_4); // multiply to set <= 0.8 to 0

Map.addLayer(multiplied,{min: -1, max: 1},'ndvi_rice_1'); // plot

I know that there are many ways of select range of values, i´m using less than [lt], greater than [gt].


texts: How to manage QGIS styles saved in Spatialite

texts:
At QGIS 2.14 the style files we have saved in our Spatialite database do not show the labels anymore.

Luckily I have the Layer Definition Files saved that still display the labels, but in trying to re-save (overwrite) the style to the Spatialite database (using the same name), I get the following:

While I'm able to save with a new name, how can I remove the old styles so I don't end up with too many old styles?

EDIT: Turns out the style saved for 2.14 doesn't show the labels either, so maybe this is a bigger issue...